001:       SUBROUTINE CSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       CHARACTER          COMPZ
009:       INTEGER            INFO, LDZ, N
010: *     ..
011: *     .. Array Arguments ..
012:       REAL               D( * ), E( * ), WORK( * )
013:       COMPLEX            Z( LDZ, * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  CSTEQR computes all eigenvalues and, optionally, eigenvectors of a
020: *  symmetric tridiagonal matrix using the implicit QL or QR method.
021: *  The eigenvectors of a full or band complex Hermitian matrix can also
022: *  be found if CHETRD or CHPTRD or CHBTRD has been used to reduce this
023: *  matrix to tridiagonal form.
024: *
025: *  Arguments
026: *  =========
027: *
028: *  COMPZ   (input) CHARACTER*1
029: *          = 'N':  Compute eigenvalues only.
030: *          = 'V':  Compute eigenvalues and eigenvectors of the original
031: *                  Hermitian matrix.  On entry, Z must contain the
032: *                  unitary matrix used to reduce the original matrix
033: *                  to tridiagonal form.
034: *          = 'I':  Compute eigenvalues and eigenvectors of the
035: *                  tridiagonal matrix.  Z is initialized to the identity
036: *                  matrix.
037: *
038: *  N       (input) INTEGER
039: *          The order of the matrix.  N >= 0.
040: *
041: *  D       (input/output) REAL array, dimension (N)
042: *          On entry, the diagonal elements of the tridiagonal matrix.
043: *          On exit, if INFO = 0, the eigenvalues in ascending order.
044: *
045: *  E       (input/output) REAL array, dimension (N-1)
046: *          On entry, the (n-1) subdiagonal elements of the tridiagonal
047: *          matrix.
048: *          On exit, E has been destroyed.
049: *
050: *  Z       (input/output) COMPLEX array, dimension (LDZ, N)
051: *          On entry, if  COMPZ = 'V', then Z contains the unitary
052: *          matrix used in the reduction to tridiagonal form.
053: *          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
054: *          orthonormal eigenvectors of the original Hermitian matrix,
055: *          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
056: *          of the symmetric tridiagonal matrix.
057: *          If COMPZ = 'N', then Z is not referenced.
058: *
059: *  LDZ     (input) INTEGER
060: *          The leading dimension of the array Z.  LDZ >= 1, and if
061: *          eigenvectors are desired, then  LDZ >= max(1,N).
062: *
063: *  WORK    (workspace) REAL array, dimension (max(1,2*N-2))
064: *          If COMPZ = 'N', then WORK is not referenced.
065: *
066: *  INFO    (output) INTEGER
067: *          = 0:  successful exit
068: *          < 0:  if INFO = -i, the i-th argument had an illegal value
069: *          > 0:  the algorithm has failed to find all the eigenvalues in
070: *                a total of 30*N iterations; if INFO = i, then i
071: *                elements of E have not converged to zero; on exit, D
072: *                and E contain the elements of a symmetric tridiagonal
073: *                matrix which is unitarily similar to the original
074: *                matrix.
075: *
076: *  =====================================================================
077: *
078: *     .. Parameters ..
079:       REAL               ZERO, ONE, TWO, THREE
080:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
081:      $                   THREE = 3.0E0 )
082:       COMPLEX            CZERO, CONE
083:       PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ),
084:      $                   CONE = ( 1.0E0, 0.0E0 ) )
085:       INTEGER            MAXIT
086:       PARAMETER          ( MAXIT = 30 )
087: *     ..
088: *     .. Local Scalars ..
089:       INTEGER            I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
090:      $                   LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
091:      $                   NM1, NMAXIT
092:       REAL               ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
093:      $                   S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
094: *     ..
095: *     .. External Functions ..
096:       LOGICAL            LSAME
097:       REAL               SLAMCH, SLANST, SLAPY2
098:       EXTERNAL           LSAME, SLAMCH, SLANST, SLAPY2
099: *     ..
100: *     .. External Subroutines ..
101:       EXTERNAL           CLASET, CLASR, CSWAP, SLAE2, SLAEV2, SLARTG,
102:      $                   SLASCL, SLASRT, XERBLA
103: *     ..
104: *     .. Intrinsic Functions ..
105:       INTRINSIC          ABS, MAX, SIGN, SQRT
106: *     ..
107: *     .. Executable Statements ..
108: *
109: *     Test the input parameters.
110: *
111:       INFO = 0
112: *
113:       IF( LSAME( COMPZ, 'N' ) ) THEN
114:          ICOMPZ = 0
115:       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
116:          ICOMPZ = 1
117:       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
118:          ICOMPZ = 2
119:       ELSE
120:          ICOMPZ = -1
121:       END IF
122:       IF( ICOMPZ.LT.0 ) THEN
123:          INFO = -1
124:       ELSE IF( N.LT.0 ) THEN
125:          INFO = -2
126:       ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
127:      $         N ) ) ) THEN
128:          INFO = -6
129:       END IF
130:       IF( INFO.NE.0 ) THEN
131:          CALL XERBLA( 'CSTEQR', -INFO )
132:          RETURN
133:       END IF
134: *
135: *     Quick return if possible
136: *
137:       IF( N.EQ.0 )
138:      $   RETURN
139: *
140:       IF( N.EQ.1 ) THEN
141:          IF( ICOMPZ.EQ.2 )
142:      $      Z( 1, 1 ) = CONE
143:          RETURN
144:       END IF
145: *
146: *     Determine the unit roundoff and over/underflow thresholds.
147: *
148:       EPS = SLAMCH( 'E' )
149:       EPS2 = EPS**2
150:       SAFMIN = SLAMCH( 'S' )
151:       SAFMAX = ONE / SAFMIN
152:       SSFMAX = SQRT( SAFMAX ) / THREE
153:       SSFMIN = SQRT( SAFMIN ) / EPS2
154: *
155: *     Compute the eigenvalues and eigenvectors of the tridiagonal
156: *     matrix.
157: *
158:       IF( ICOMPZ.EQ.2 )
159:      $   CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
160: *
161:       NMAXIT = N*MAXIT
162:       JTOT = 0
163: *
164: *     Determine where the matrix splits and choose QL or QR iteration
165: *     for each block, according to whether top or bottom diagonal
166: *     element is smaller.
167: *
168:       L1 = 1
169:       NM1 = N - 1
170: *
171:    10 CONTINUE
172:       IF( L1.GT.N )
173:      $   GO TO 160
174:       IF( L1.GT.1 )
175:      $   E( L1-1 ) = ZERO
176:       IF( L1.LE.NM1 ) THEN
177:          DO 20 M = L1, NM1
178:             TST = ABS( E( M ) )
179:             IF( TST.EQ.ZERO )
180:      $         GO TO 30
181:             IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
182:      $          1 ) ) ) )*EPS ) THEN
183:                E( M ) = ZERO
184:                GO TO 30
185:             END IF
186:    20    CONTINUE
187:       END IF
188:       M = N
189: *
190:    30 CONTINUE
191:       L = L1
192:       LSV = L
193:       LEND = M
194:       LENDSV = LEND
195:       L1 = M + 1
196:       IF( LEND.EQ.L )
197:      $   GO TO 10
198: *
199: *     Scale submatrix in rows and columns L to LEND
200: *
201:       ANORM = SLANST( 'I', LEND-L+1, D( L ), E( L ) )
202:       ISCALE = 0
203:       IF( ANORM.EQ.ZERO )
204:      $   GO TO 10
205:       IF( ANORM.GT.SSFMAX ) THEN
206:          ISCALE = 1
207:          CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
208:      $                INFO )
209:          CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
210:      $                INFO )
211:       ELSE IF( ANORM.LT.SSFMIN ) THEN
212:          ISCALE = 2
213:          CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
214:      $                INFO )
215:          CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
216:      $                INFO )
217:       END IF
218: *
219: *     Choose between QL and QR iteration
220: *
221:       IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
222:          LEND = LSV
223:          L = LENDSV
224:       END IF
225: *
226:       IF( LEND.GT.L ) THEN
227: *
228: *        QL Iteration
229: *
230: *        Look for small subdiagonal element.
231: *
232:    40    CONTINUE
233:          IF( L.NE.LEND ) THEN
234:             LENDM1 = LEND - 1
235:             DO 50 M = L, LENDM1
236:                TST = ABS( E( M ) )**2
237:                IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+
238:      $             SAFMIN )GO TO 60
239:    50       CONTINUE
240:          END IF
241: *
242:          M = LEND
243: *
244:    60    CONTINUE
245:          IF( M.LT.LEND )
246:      $      E( M ) = ZERO
247:          P = D( L )
248:          IF( M.EQ.L )
249:      $      GO TO 80
250: *
251: *        If remaining matrix is 2-by-2, use SLAE2 or SLAEV2
252: *        to compute its eigensystem.
253: *
254:          IF( M.EQ.L+1 ) THEN
255:             IF( ICOMPZ.GT.0 ) THEN
256:                CALL SLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S )
257:                WORK( L ) = C
258:                WORK( N-1+L ) = S
259:                CALL CLASR( 'R', 'V', 'B', N, 2, WORK( L ),
260:      $                     WORK( N-1+L ), Z( 1, L ), LDZ )
261:             ELSE
262:                CALL SLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 )
263:             END IF
264:             D( L ) = RT1
265:             D( L+1 ) = RT2
266:             E( L ) = ZERO
267:             L = L + 2
268:             IF( L.LE.LEND )
269:      $         GO TO 40
270:             GO TO 140
271:          END IF
272: *
273:          IF( JTOT.EQ.NMAXIT )
274:      $      GO TO 140
275:          JTOT = JTOT + 1
276: *
277: *        Form shift.
278: *
279:          G = ( D( L+1 )-P ) / ( TWO*E( L ) )
280:          R = SLAPY2( G, ONE )
281:          G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) )
282: *
283:          S = ONE
284:          C = ONE
285:          P = ZERO
286: *
287: *        Inner loop
288: *
289:          MM1 = M - 1
290:          DO 70 I = MM1, L, -1
291:             F = S*E( I )
292:             B = C*E( I )
293:             CALL SLARTG( G, F, C, S, R )
294:             IF( I.NE.M-1 )
295:      $         E( I+1 ) = R
296:             G = D( I+1 ) - P
297:             R = ( D( I )-G )*S + TWO*C*B
298:             P = S*R
299:             D( I+1 ) = G + P
300:             G = C*R - B
301: *
302: *           If eigenvectors are desired, then save rotations.
303: *
304:             IF( ICOMPZ.GT.0 ) THEN
305:                WORK( I ) = C
306:                WORK( N-1+I ) = -S
307:             END IF
308: *
309:    70    CONTINUE
310: *
311: *        If eigenvectors are desired, then apply saved rotations.
312: *
313:          IF( ICOMPZ.GT.0 ) THEN
314:             MM = M - L + 1
315:             CALL CLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ),
316:      $                  Z( 1, L ), LDZ )
317:          END IF
318: *
319:          D( L ) = D( L ) - P
320:          E( L ) = G
321:          GO TO 40
322: *
323: *        Eigenvalue found.
324: *
325:    80    CONTINUE
326:          D( L ) = P
327: *
328:          L = L + 1
329:          IF( L.LE.LEND )
330:      $      GO TO 40
331:          GO TO 140
332: *
333:       ELSE
334: *
335: *        QR Iteration
336: *
337: *        Look for small superdiagonal element.
338: *
339:    90    CONTINUE
340:          IF( L.NE.LEND ) THEN
341:             LENDP1 = LEND + 1
342:             DO 100 M = L, LENDP1, -1
343:                TST = ABS( E( M-1 ) )**2
344:                IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+
345:      $             SAFMIN )GO TO 110
346:   100       CONTINUE
347:          END IF
348: *
349:          M = LEND
350: *
351:   110    CONTINUE
352:          IF( M.GT.LEND )
353:      $      E( M-1 ) = ZERO
354:          P = D( L )
355:          IF( M.EQ.L )
356:      $      GO TO 130
357: *
358: *        If remaining matrix is 2-by-2, use SLAE2 or SLAEV2
359: *        to compute its eigensystem.
360: *
361:          IF( M.EQ.L-1 ) THEN
362:             IF( ICOMPZ.GT.0 ) THEN
363:                CALL SLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S )
364:                WORK( M ) = C
365:                WORK( N-1+M ) = S
366:                CALL CLASR( 'R', 'V', 'F', N, 2, WORK( M ),
367:      $                     WORK( N-1+M ), Z( 1, L-1 ), LDZ )
368:             ELSE
369:                CALL SLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 )
370:             END IF
371:             D( L-1 ) = RT1
372:             D( L ) = RT2
373:             E( L-1 ) = ZERO
374:             L = L - 2
375:             IF( L.GE.LEND )
376:      $         GO TO 90
377:             GO TO 140
378:          END IF
379: *
380:          IF( JTOT.EQ.NMAXIT )
381:      $      GO TO 140
382:          JTOT = JTOT + 1
383: *
384: *        Form shift.
385: *
386:          G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
387:          R = SLAPY2( G, ONE )
388:          G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
389: *
390:          S = ONE
391:          C = ONE
392:          P = ZERO
393: *
394: *        Inner loop
395: *
396:          LM1 = L - 1
397:          DO 120 I = M, LM1
398:             F = S*E( I )
399:             B = C*E( I )
400:             CALL SLARTG( G, F, C, S, R )
401:             IF( I.NE.M )
402:      $         E( I-1 ) = R
403:             G = D( I ) - P
404:             R = ( D( I+1 )-G )*S + TWO*C*B
405:             P = S*R
406:             D( I ) = G + P
407:             G = C*R - B
408: *
409: *           If eigenvectors are desired, then save rotations.
410: *
411:             IF( ICOMPZ.GT.0 ) THEN
412:                WORK( I ) = C
413:                WORK( N-1+I ) = S
414:             END IF
415: *
416:   120    CONTINUE
417: *
418: *        If eigenvectors are desired, then apply saved rotations.
419: *
420:          IF( ICOMPZ.GT.0 ) THEN
421:             MM = L - M + 1
422:             CALL CLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ),
423:      $                  Z( 1, M ), LDZ )
424:          END IF
425: *
426:          D( L ) = D( L ) - P
427:          E( LM1 ) = G
428:          GO TO 90
429: *
430: *        Eigenvalue found.
431: *
432:   130    CONTINUE
433:          D( L ) = P
434: *
435:          L = L - 1
436:          IF( L.GE.LEND )
437:      $      GO TO 90
438:          GO TO 140
439: *
440:       END IF
441: *
442: *     Undo scaling if necessary
443: *
444:   140 CONTINUE
445:       IF( ISCALE.EQ.1 ) THEN
446:          CALL SLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
447:      $                D( LSV ), N, INFO )
448:          CALL SLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ),
449:      $                N, INFO )
450:       ELSE IF( ISCALE.EQ.2 ) THEN
451:          CALL SLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
452:      $                D( LSV ), N, INFO )
453:          CALL SLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ),
454:      $                N, INFO )
455:       END IF
456: *
457: *     Check for no convergence to an eigenvalue after a total
458: *     of N*MAXIT iterations.
459: *
460:       IF( JTOT.EQ.NMAXIT ) THEN
461:          DO 150 I = 1, N - 1
462:             IF( E( I ).NE.ZERO )
463:      $         INFO = INFO + 1
464:   150    CONTINUE
465:          RETURN
466:       END IF
467:       GO TO 10
468: *
469: *     Order eigenvalues and eigenvectors.
470: *
471:   160 CONTINUE
472:       IF( ICOMPZ.EQ.0 ) THEN
473: *
474: *        Use Quick Sort
475: *
476:          CALL SLASRT( 'I', N, D, INFO )
477: *
478:       ELSE
479: *
480: *        Use Selection Sort to minimize swaps of eigenvectors
481: *
482:          DO 180 II = 2, N
483:             I = II - 1
484:             K = I
485:             P = D( I )
486:             DO 170 J = II, N
487:                IF( D( J ).LT.P ) THEN
488:                   K = J
489:                   P = D( J )
490:                END IF
491:   170       CONTINUE
492:             IF( K.NE.I ) THEN
493:                D( K ) = D( I )
494:                D( I ) = P
495:                CALL CSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
496:             END IF
497:   180    CONTINUE
498:       END IF
499:       RETURN
500: *
501: *     End of CSTEQR
502: *
503:       END
504: