001:       SUBROUTINE SSTERF( N, D, E, 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:       INTEGER            INFO, N
009: *     ..
010: *     .. Array Arguments ..
011:       REAL               D( * ), E( * )
012: *     ..
013: *
014: *  Purpose
015: *  =======
016: *
017: *  SSTERF computes all eigenvalues of a symmetric tridiagonal matrix
018: *  using the Pal-Walker-Kahan variant of the QL or QR algorithm.
019: *
020: *  Arguments
021: *  =========
022: *
023: *  N       (input) INTEGER
024: *          The order of the matrix.  N >= 0.
025: *
026: *  D       (input/output) REAL array, dimension (N)
027: *          On entry, the n diagonal elements of the tridiagonal matrix.
028: *          On exit, if INFO = 0, the eigenvalues in ascending order.
029: *
030: *  E       (input/output) REAL array, dimension (N-1)
031: *          On entry, the (n-1) subdiagonal elements of the tridiagonal
032: *          matrix.
033: *          On exit, E has been destroyed.
034: *
035: *  INFO    (output) INTEGER
036: *          = 0:  successful exit
037: *          < 0:  if INFO = -i, the i-th argument had an illegal value
038: *          > 0:  the algorithm failed to find all of the eigenvalues in
039: *                a total of 30*N iterations; if INFO = i, then i
040: *                elements of E have not converged to zero.
041: *
042: *  =====================================================================
043: *
044: *     .. Parameters ..
045:       REAL               ZERO, ONE, TWO, THREE
046:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
047:      $                   THREE = 3.0E0 )
048:       INTEGER            MAXIT
049:       PARAMETER          ( MAXIT = 30 )
050: *     ..
051: *     .. Local Scalars ..
052:       INTEGER            I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M,
053:      $                   NMAXIT
054:       REAL               ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
055:      $                   OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
056:      $                   SIGMA, SSFMAX, SSFMIN
057: *     ..
058: *     .. External Functions ..
059:       REAL               SLAMCH, SLANST, SLAPY2
060:       EXTERNAL           SLAMCH, SLANST, SLAPY2
061: *     ..
062: *     .. External Subroutines ..
063:       EXTERNAL           SLAE2, SLASCL, SLASRT, XERBLA
064: *     ..
065: *     .. Intrinsic Functions ..
066:       INTRINSIC          ABS, SIGN, SQRT
067: *     ..
068: *     .. Executable Statements ..
069: *
070: *     Test the input parameters.
071: *
072:       INFO = 0
073: *
074: *     Quick return if possible
075: *
076:       IF( N.LT.0 ) THEN
077:          INFO = -1
078:          CALL XERBLA( 'SSTERF', -INFO )
079:          RETURN
080:       END IF
081:       IF( N.LE.1 )
082:      $   RETURN
083: *
084: *     Determine the unit roundoff for this environment.
085: *
086:       EPS = SLAMCH( 'E' )
087:       EPS2 = EPS**2
088:       SAFMIN = SLAMCH( 'S' )
089:       SAFMAX = ONE / SAFMIN
090:       SSFMAX = SQRT( SAFMAX ) / THREE
091:       SSFMIN = SQRT( SAFMIN ) / EPS2
092: *
093: *     Compute the eigenvalues of the tridiagonal matrix.
094: *
095:       NMAXIT = N*MAXIT
096:       SIGMA = ZERO
097:       JTOT = 0
098: *
099: *     Determine where the matrix splits and choose QL or QR iteration
100: *     for each block, according to whether top or bottom diagonal
101: *     element is smaller.
102: *
103:       L1 = 1
104: *
105:    10 CONTINUE
106:       IF( L1.GT.N )
107:      $   GO TO 170
108:       IF( L1.GT.1 )
109:      $   E( L1-1 ) = ZERO
110:       DO 20 M = L1, N - 1
111:          IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*
112:      $       SQRT( ABS( D( M+1 ) ) ) )*EPS ) THEN
113:             E( M ) = ZERO
114:             GO TO 30
115:          END IF
116:    20 CONTINUE
117:       M = N
118: *
119:    30 CONTINUE
120:       L = L1
121:       LSV = L
122:       LEND = M
123:       LENDSV = LEND
124:       L1 = M + 1
125:       IF( LEND.EQ.L )
126:      $   GO TO 10
127: *
128: *     Scale submatrix in rows and columns L to LEND
129: *
130:       ANORM = SLANST( 'I', LEND-L+1, D( L ), E( L ) )
131:       ISCALE = 0
132:       IF( ANORM.GT.SSFMAX ) THEN
133:          ISCALE = 1
134:          CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
135:      $                INFO )
136:          CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
137:      $                INFO )
138:       ELSE IF( ANORM.LT.SSFMIN ) THEN
139:          ISCALE = 2
140:          CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
141:      $                INFO )
142:          CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
143:      $                INFO )
144:       END IF
145: *
146:       DO 40 I = L, LEND - 1
147:          E( I ) = E( I )**2
148:    40 CONTINUE
149: *
150: *     Choose between QL and QR iteration
151: *
152:       IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
153:          LEND = LSV
154:          L = LENDSV
155:       END IF
156: *
157:       IF( LEND.GE.L ) THEN
158: *
159: *        QL Iteration
160: *
161: *        Look for small subdiagonal element.
162: *
163:    50    CONTINUE
164:          IF( L.NE.LEND ) THEN
165:             DO 60 M = L, LEND - 1
166:                IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) )
167:      $            GO TO 70
168:    60       CONTINUE
169:          END IF
170:          M = LEND
171: *
172:    70    CONTINUE
173:          IF( M.LT.LEND )
174:      $      E( M ) = ZERO
175:          P = D( L )
176:          IF( M.EQ.L )
177:      $      GO TO 90
178: *
179: *        If remaining matrix is 2 by 2, use SLAE2 to compute its
180: *        eigenvalues.
181: *
182:          IF( M.EQ.L+1 ) THEN
183:             RTE = SQRT( E( L ) )
184:             CALL SLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 )
185:             D( L ) = RT1
186:             D( L+1 ) = RT2
187:             E( L ) = ZERO
188:             L = L + 2
189:             IF( L.LE.LEND )
190:      $         GO TO 50
191:             GO TO 150
192:          END IF
193: *
194:          IF( JTOT.EQ.NMAXIT )
195:      $      GO TO 150
196:          JTOT = JTOT + 1
197: *
198: *        Form shift.
199: *
200:          RTE = SQRT( E( L ) )
201:          SIGMA = ( D( L+1 )-P ) / ( TWO*RTE )
202:          R = SLAPY2( SIGMA, ONE )
203:          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
204: *
205:          C = ONE
206:          S = ZERO
207:          GAMMA = D( M ) - SIGMA
208:          P = GAMMA*GAMMA
209: *
210: *        Inner loop
211: *
212:          DO 80 I = M - 1, L, -1
213:             BB = E( I )
214:             R = P + BB
215:             IF( I.NE.M-1 )
216:      $         E( I+1 ) = S*R
217:             OLDC = C
218:             C = P / R
219:             S = BB / R
220:             OLDGAM = GAMMA
221:             ALPHA = D( I )
222:             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
223:             D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )
224:             IF( C.NE.ZERO ) THEN
225:                P = ( GAMMA*GAMMA ) / C
226:             ELSE
227:                P = OLDC*BB
228:             END IF
229:    80    CONTINUE
230: *
231:          E( L ) = S*P
232:          D( L ) = SIGMA + GAMMA
233:          GO TO 50
234: *
235: *        Eigenvalue found.
236: *
237:    90    CONTINUE
238:          D( L ) = P
239: *
240:          L = L + 1
241:          IF( L.LE.LEND )
242:      $      GO TO 50
243:          GO TO 150
244: *
245:       ELSE
246: *
247: *        QR Iteration
248: *
249: *        Look for small superdiagonal element.
250: *
251:   100    CONTINUE
252:          DO 110 M = L, LEND + 1, -1
253:             IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) )
254:      $         GO TO 120
255:   110    CONTINUE
256:          M = LEND
257: *
258:   120    CONTINUE
259:          IF( M.GT.LEND )
260:      $      E( M-1 ) = ZERO
261:          P = D( L )
262:          IF( M.EQ.L )
263:      $      GO TO 140
264: *
265: *        If remaining matrix is 2 by 2, use SLAE2 to compute its
266: *        eigenvalues.
267: *
268:          IF( M.EQ.L-1 ) THEN
269:             RTE = SQRT( E( L-1 ) )
270:             CALL SLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 )
271:             D( L ) = RT1
272:             D( L-1 ) = RT2
273:             E( L-1 ) = ZERO
274:             L = L - 2
275:             IF( L.GE.LEND )
276:      $         GO TO 100
277:             GO TO 150
278:          END IF
279: *
280:          IF( JTOT.EQ.NMAXIT )
281:      $      GO TO 150
282:          JTOT = JTOT + 1
283: *
284: *        Form shift.
285: *
286:          RTE = SQRT( E( L-1 ) )
287:          SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )
288:          R = SLAPY2( SIGMA, ONE )
289:          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
290: *
291:          C = ONE
292:          S = ZERO
293:          GAMMA = D( M ) - SIGMA
294:          P = GAMMA*GAMMA
295: *
296: *        Inner loop
297: *
298:          DO 130 I = M, L - 1
299:             BB = E( I )
300:             R = P + BB
301:             IF( I.NE.M )
302:      $         E( I-1 ) = S*R
303:             OLDC = C
304:             C = P / R
305:             S = BB / R
306:             OLDGAM = GAMMA
307:             ALPHA = D( I+1 )
308:             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
309:             D( I ) = OLDGAM + ( ALPHA-GAMMA )
310:             IF( C.NE.ZERO ) THEN
311:                P = ( GAMMA*GAMMA ) / C
312:             ELSE
313:                P = OLDC*BB
314:             END IF
315:   130    CONTINUE
316: *
317:          E( L-1 ) = S*P
318:          D( L ) = SIGMA + GAMMA
319:          GO TO 100
320: *
321: *        Eigenvalue found.
322: *
323:   140    CONTINUE
324:          D( L ) = P
325: *
326:          L = L - 1
327:          IF( L.GE.LEND )
328:      $      GO TO 100
329:          GO TO 150
330: *
331:       END IF
332: *
333: *     Undo scaling if necessary
334: *
335:   150 CONTINUE
336:       IF( ISCALE.EQ.1 )
337:      $   CALL SLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
338:      $                D( LSV ), N, INFO )
339:       IF( ISCALE.EQ.2 )
340:      $   CALL SLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
341:      $                D( LSV ), N, INFO )
342: *
343: *     Check for no convergence to an eigenvalue after a total
344: *     of N*MAXIT iterations.
345: *
346:       IF( JTOT.LT.NMAXIT )
347:      $   GO TO 10
348:       DO 160 I = 1, N - 1
349:          IF( E( I ).NE.ZERO )
350:      $      INFO = INFO + 1
351:   160 CONTINUE
352:       GO TO 180
353: *
354: *     Sort eigenvalues in increasing order.
355: *
356:   170 CONTINUE
357:       CALL SLASRT( 'I', N, D, INFO )
358: *
359:   180 CONTINUE
360:       RETURN
361: *
362: *     End of SSTERF
363: *
364:       END
365: