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