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