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