1       SUBROUTINE DSTERF( N, D, E, INFO )
  2 *
  3 *  -- LAPACK routine (version 3.3.1) --
  4 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  5 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  6 *  -- April 2011                                                      --
  7 *
  8 *     .. Scalar Arguments ..
  9       INTEGER            INFO, N
 10 *     ..
 11 *     .. Array Arguments ..
 12       DOUBLE PRECISION   D( * ), E( * )
 13 *     ..
 14 *
 15 *  Purpose
 16 *  =======
 17 *
 18 *  DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
 19 *  using the Pal-Walker-Kahan variant of the QL or QR algorithm.
 20 *
 21 *  Arguments
 22 *  =========
 23 *
 24 *  N       (input) INTEGER
 25 *          The order of the matrix.  N >= 0.
 26 *
 27 *  D       (input/output) DOUBLE PRECISION array, dimension (N)
 28 *          On entry, the n diagonal elements of the tridiagonal matrix.
 29 *          On exit, if INFO = 0, the eigenvalues in ascending order.
 30 *
 31 *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
 32 *          On entry, the (n-1) subdiagonal elements of the tridiagonal
 33 *          matrix.
 34 *          On exit, E has been destroyed.
 35 *
 36 *  INFO    (output) INTEGER
 37 *          = 0:  successful exit
 38 *          < 0:  if INFO = -i, the i-th argument had an illegal value
 39 *          > 0:  the algorithm failed to find all of the eigenvalues in
 40 *                a total of 30*N iterations; if INFO = i, then i
 41 *                elements of E have not converged to zero.
 42 *
 43 *  =====================================================================
 44 *
 45 *     .. Parameters ..
 46       DOUBLE PRECISION   ZERO, ONE, TWO, THREE
 47       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
 48      $                   THREE = 3.0D0 )
 49       INTEGER            MAXIT
 50       PARAMETER          ( MAXIT = 30 )
 51 *     ..
 52 *     .. Local Scalars ..
 53       INTEGER            I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M,
 54      $                   NMAXIT
 55       DOUBLE PRECISION   ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
 56      $                   OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
 57      $                   SIGMA, SSFMAX, SSFMIN, RMAX
 58 *     ..
 59 *     .. External Functions ..
 60       DOUBLE PRECISION   DLAMCH, DLANST, DLAPY2
 61       EXTERNAL           DLAMCH, DLANST, DLAPY2
 62 *     ..
 63 *     .. External Subroutines ..
 64       EXTERNAL           DLAE2, DLASCL, DLASRT, XERBLA
 65 *     ..
 66 *     .. Intrinsic Functions ..
 67       INTRINSIC          ABSSIGNSQRT
 68 *     ..
 69 *     .. Executable Statements ..
 70 *
 71 *     Test the input parameters.
 72 *
 73       INFO = 0
 74 *
 75 *     Quick return if possible
 76 *
 77       IF( N.LT.0 ) THEN
 78          INFO = -1
 79          CALL XERBLA( 'DSTERF'-INFO )
 80          RETURN
 81       END IF
 82       IF( N.LE.1 )
 83      $   RETURN
 84 *
 85 *     Determine the unit roundoff for this environment.
 86 *
 87       EPS = DLAMCH( 'E' )
 88       EPS2 = EPS**2
 89       SAFMIN = DLAMCH( 'S' )
 90       SAFMAX = ONE / SAFMIN
 91       SSFMAX = SQRT( SAFMAX ) / THREE
 92       SSFMIN = SQRT( SAFMIN ) / EPS2
 93       RMAX = DLAMCH( 'O' )
 94 *
 95 *     Compute the eigenvalues of the tridiagonal matrix.
 96 *
 97       NMAXIT = N*MAXIT
 98       SIGMA = ZERO
 99       JTOT = 0
100 *
101 *     Determine where the matrix splits and choose QL or QR iteration
102 *     for each block, according to whether top or bottom diagonal
103 *     element is smaller.
104 *
105       L1 = 1
106 *
107    10 CONTINUE
108       IF( L1.GT.N )
109      $   GO TO 170
110       IF( L1.GT.1 )
111      $   E( L1-1 ) = ZERO
112       DO 20 M = L1, N - 1
113          IFABS( E( M ) ).LE.SQRTABS( D( M ) ) )*SQRTABS( D( M+
114      $       1 ) ) ) )*EPS ) THEN
115             E( M ) = ZERO
116             GO TO 30
117          END IF
118    20 CONTINUE
119       M = N
120 *
121    30 CONTINUE
122       L = L1
123       LSV = L
124       LEND = M
125       LENDSV = LEND
126       L1 = M + 1
127       IF( LEND.EQ.L )
128      $   GO TO 10
129 *
130 *     Scale submatrix in rows and columns L to LEND
131 *
132       ANORM = DLANST( 'M', LEND-L+1, D( L ), E( L ) )
133       ISCALE = 0
134       IF( ANORM.EQ.ZERO )
135      $   GO TO 10      
136       IF( (ANORM.GT.SSFMAX) ) THEN
137          ISCALE = 1
138          CALL DLASCL( 'G'00, ANORM, SSFMAX, LEND-L+11, D( L ), N,
139      $                INFO )
140          CALL DLASCL( 'G'00, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
141      $                INFO )
142       ELSE IF( ANORM.LT.SSFMIN ) THEN
143          ISCALE = 2
144          CALL DLASCL( 'G'00, ANORM, SSFMIN, LEND-L+11, D( L ), N,
145      $                INFO )
146          CALL DLASCL( 'G'00, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
147      $                INFO )
148       END IF
149 *
150       DO 40 I = L, LEND - 1
151          E( I ) = E( I )**2
152    40 CONTINUE
153 *
154 *     Choose between QL and QR iteration
155 *
156       IFABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
157          LEND = LSV
158          L = LENDSV
159       END IF
160 *
161       IF( LEND.GE.L ) THEN
162 *
163 *        QL Iteration
164 *
165 *        Look for small subdiagonal element.
166 *
167    50    CONTINUE
168          IF( L.NE.LEND ) THEN
169             DO 60 M = L, LEND - 1
170                IFABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) )
171      $            GO TO 70
172    60       CONTINUE
173          END IF
174          M = LEND
175 *
176    70    CONTINUE
177          IF( M.LT.LEND )
178      $      E( M ) = ZERO
179          P = D( L )
180          IF( M.EQ.L )
181      $      GO TO 90
182 *
183 *        If remaining matrix is 2 by 2, use DLAE2 to compute its
184 *        eigenvalues.
185 *
186          IF( M.EQ.L+1 ) THEN
187             RTE = SQRT( E( L ) )
188             CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 )
189             D( L ) = RT1
190             D( L+1 ) = RT2
191             E( L ) = ZERO
192             L = L + 2
193             IF( L.LE.LEND )
194      $         GO TO 50
195             GO TO 150
196          END IF
197 *
198          IF( JTOT.EQ.NMAXIT )
199      $      GO TO 150
200          JTOT = JTOT + 1
201 *
202 *        Form shift.
203 *
204          RTE = SQRT( E( L ) )
205          SIGMA = ( D( L+1 )-P ) / ( TWO*RTE )
206          R = DLAPY2( SIGMA, ONE )
207          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
208 *
209          C = ONE
210          S = ZERO
211          GAMMA = D( M ) - SIGMA
212          P = GAMMA*GAMMA
213 *
214 *        Inner loop
215 *
216          DO 80 I = M - 1, L, -1
217             BB = E( I )
218             R = P + BB
219             IF( I.NE.M-1 )
220      $         E( I+1 ) = S*R
221             OLDC = C
222             C = P / R
223             S = BB / R
224             OLDGAM = GAMMA
225             ALPHA = D( I )
226             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
227             D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )
228             IF( C.NE.ZERO ) THEN
229                P = ( GAMMA*GAMMA ) / C
230             ELSE
231                P = OLDC*BB
232             END IF
233    80    CONTINUE
234 *
235          E( L ) = S*P
236          D( L ) = SIGMA + GAMMA
237          GO TO 50
238 *
239 *        Eigenvalue found.
240 *
241    90    CONTINUE
242          D( L ) = P
243 *
244          L = L + 1
245          IF( L.LE.LEND )
246      $      GO TO 50
247          GO TO 150
248 *
249       ELSE
250 *
251 *        QR Iteration
252 *
253 *        Look for small superdiagonal element.
254 *
255   100    CONTINUE
256          DO 110 M = L, LEND + 1-1
257             IFABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) )
258      $         GO TO 120
259   110    CONTINUE
260          M = LEND
261 *
262   120    CONTINUE
263          IF( M.GT.LEND )
264      $      E( M-1 ) = ZERO
265          P = D( L )
266          IF( M.EQ.L )
267      $      GO TO 140
268 *
269 *        If remaining matrix is 2 by 2, use DLAE2 to compute its
270 *        eigenvalues.
271 *
272          IF( M.EQ.L-1 ) THEN
273             RTE = SQRT( E( L-1 ) )
274             CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 )
275             D( L ) = RT1
276             D( L-1 ) = RT2
277             E( L-1 ) = ZERO
278             L = L - 2
279             IF( L.GE.LEND )
280      $         GO TO 100
281             GO TO 150
282          END IF
283 *
284          IF( JTOT.EQ.NMAXIT )
285      $      GO TO 150
286          JTOT = JTOT + 1
287 *
288 *        Form shift.
289 *
290          RTE = SQRT( E( L-1 ) )
291          SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )
292          R = DLAPY2( SIGMA, ONE )
293          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
294 *
295          C = ONE
296          S = ZERO
297          GAMMA = D( M ) - SIGMA
298          P = GAMMA*GAMMA
299 *
300 *        Inner loop
301 *
302          DO 130 I = M, L - 1
303             BB = E( I )
304             R = P + BB
305             IF( I.NE.M )
306      $         E( I-1 ) = S*R
307             OLDC = C
308             C = P / R
309             S = BB / R
310             OLDGAM = GAMMA
311             ALPHA = D( I+1 )
312             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
313             D( I ) = OLDGAM + ( ALPHA-GAMMA )
314             IF( C.NE.ZERO ) THEN
315                P = ( GAMMA*GAMMA ) / C
316             ELSE
317                P = OLDC*BB
318             END IF
319   130    CONTINUE
320 *
321          E( L-1 ) = S*P
322          D( L ) = SIGMA + GAMMA
323          GO TO 100
324 *
325 *        Eigenvalue found.
326 *
327   140    CONTINUE
328          D( L ) = P
329 *
330          L = L - 1
331          IF( L.GE.LEND )
332      $      GO TO 100
333          GO TO 150
334 *
335       END IF
336 *
337 *     Undo scaling if necessary
338 *
339   150 CONTINUE
340       IF( ISCALE.EQ.1 )
341      $   CALL DLASCL( 'G'00, SSFMAX, ANORM, LENDSV-LSV+11,
342      $                D( LSV ), N, INFO )
343       IF( ISCALE.EQ.2 )
344      $   CALL DLASCL( 'G'00, SSFMIN, ANORM, LENDSV-LSV+11,
345      $                D( LSV ), N, INFO )
346 *
347 *     Check for no convergence to an eigenvalue after a total
348 *     of N*MAXIT iterations.
349 *
350       IF( JTOT.LT.NMAXIT )
351      $   GO TO 10
352       DO 160 I = 1, N - 1
353          IF( E( I ).NE.ZERO )
354      $      INFO = INFO + 1
355   160 CONTINUE
356       GO TO 180
357 *
358 *     Sort eigenvalues in increasing order.
359 *
360   170 CONTINUE
361       CALL DLASRT( 'I', N, D, INFO )
362 *
363   180 CONTINUE
364       RETURN
365 *
366 *     End of DSTERF
367 *
368       END