|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DSTERF( N, D, E, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 INTEGER INFO, N 00010 * .. 00011 * .. Array Arguments .. 00012 DOUBLE PRECISION D( * ), E( * ) 00013 * .. 00014 * 00015 * Purpose 00016 * ======= 00017 * 00018 * DSTERF computes all eigenvalues of a symmetric tridiagonal matrix 00019 * using the Pal-Walker-Kahan variant of the QL or QR algorithm. 00020 * 00021 * Arguments 00022 * ========= 00023 * 00024 * N (input) INTEGER 00025 * The order of the matrix. N >= 0. 00026 * 00027 * D (input/output) DOUBLE PRECISION array, dimension (N) 00028 * On entry, the n diagonal elements of the tridiagonal matrix. 00029 * On exit, if INFO = 0, the eigenvalues in ascending order. 00030 * 00031 * E (input/output) DOUBLE PRECISION array, dimension (N-1) 00032 * On entry, the (n-1) subdiagonal elements of the tridiagonal 00033 * matrix. 00034 * On exit, E has been destroyed. 00035 * 00036 * INFO (output) INTEGER 00037 * = 0: successful exit 00038 * < 0: if INFO = -i, the i-th argument had an illegal value 00039 * > 0: the algorithm failed to find all of the eigenvalues in 00040 * a total of 30*N iterations; if INFO = i, then i 00041 * elements of E have not converged to zero. 00042 * 00043 * ===================================================================== 00044 * 00045 * .. Parameters .. 00046 DOUBLE PRECISION ZERO, ONE, TWO, THREE 00047 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 00048 $ THREE = 3.0D0 ) 00049 INTEGER MAXIT 00050 PARAMETER ( MAXIT = 30 ) 00051 * .. 00052 * .. Local Scalars .. 00053 INTEGER I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M, 00054 $ NMAXIT 00055 DOUBLE PRECISION ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC, 00056 $ OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN, 00057 $ SIGMA, SSFMAX, SSFMIN 00058 * .. 00059 * .. External Functions .. 00060 DOUBLE PRECISION DLAMCH, DLANST, DLAPY2 00061 EXTERNAL DLAMCH, DLANST, DLAPY2 00062 * .. 00063 * .. External Subroutines .. 00064 EXTERNAL DLAE2, DLASCL, DLASRT, XERBLA 00065 * .. 00066 * .. Intrinsic Functions .. 00067 INTRINSIC ABS, SIGN, SQRT 00068 * .. 00069 * .. Executable Statements .. 00070 * 00071 * Test the input parameters. 00072 * 00073 INFO = 0 00074 * 00075 * Quick return if possible 00076 * 00077 IF( N.LT.0 ) THEN 00078 INFO = -1 00079 CALL XERBLA( 'DSTERF', -INFO ) 00080 RETURN 00081 END IF 00082 IF( N.LE.1 ) 00083 $ RETURN 00084 * 00085 * Determine the unit roundoff for this environment. 00086 * 00087 EPS = DLAMCH( 'E' ) 00088 EPS2 = EPS**2 00089 SAFMIN = DLAMCH( 'S' ) 00090 SAFMAX = ONE / SAFMIN 00091 SSFMAX = SQRT( SAFMAX ) / THREE 00092 SSFMIN = SQRT( SAFMIN ) / EPS2 00093 * 00094 * Compute the eigenvalues of the tridiagonal matrix. 00095 * 00096 NMAXIT = N*MAXIT 00097 SIGMA = ZERO 00098 JTOT = 0 00099 * 00100 * Determine where the matrix splits and choose QL or QR iteration 00101 * for each block, according to whether top or bottom diagonal 00102 * element is smaller. 00103 * 00104 L1 = 1 00105 * 00106 10 CONTINUE 00107 IF( L1.GT.N ) 00108 $ GO TO 170 00109 IF( L1.GT.1 ) 00110 $ E( L1-1 ) = ZERO 00111 DO 20 M = L1, N - 1 00112 IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ 00113 $ 1 ) ) ) )*EPS ) THEN 00114 E( M ) = ZERO 00115 GO TO 30 00116 END IF 00117 20 CONTINUE 00118 M = N 00119 * 00120 30 CONTINUE 00121 L = L1 00122 LSV = L 00123 LEND = M 00124 LENDSV = LEND 00125 L1 = M + 1 00126 IF( LEND.EQ.L ) 00127 $ GO TO 10 00128 * 00129 * Scale submatrix in rows and columns L to LEND 00130 * 00131 ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) ) 00132 ISCALE = 0 00133 IF( ANORM.GT.SSFMAX ) THEN 00134 ISCALE = 1 00135 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, 00136 $ INFO ) 00137 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, 00138 $ INFO ) 00139 ELSE IF( ANORM.LT.SSFMIN ) THEN 00140 ISCALE = 2 00141 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, 00142 $ INFO ) 00143 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, 00144 $ INFO ) 00145 END IF 00146 * 00147 DO 40 I = L, LEND - 1 00148 E( I ) = E( I )**2 00149 40 CONTINUE 00150 * 00151 * Choose between QL and QR iteration 00152 * 00153 IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN 00154 LEND = LSV 00155 L = LENDSV 00156 END IF 00157 * 00158 IF( LEND.GE.L ) THEN 00159 * 00160 * QL Iteration 00161 * 00162 * Look for small subdiagonal element. 00163 * 00164 50 CONTINUE 00165 IF( L.NE.LEND ) THEN 00166 DO 60 M = L, LEND - 1 00167 IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) ) 00168 $ GO TO 70 00169 60 CONTINUE 00170 END IF 00171 M = LEND 00172 * 00173 70 CONTINUE 00174 IF( M.LT.LEND ) 00175 $ E( M ) = ZERO 00176 P = D( L ) 00177 IF( M.EQ.L ) 00178 $ GO TO 90 00179 * 00180 * If remaining matrix is 2 by 2, use DLAE2 to compute its 00181 * eigenvalues. 00182 * 00183 IF( M.EQ.L+1 ) THEN 00184 RTE = SQRT( E( L ) ) 00185 CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 ) 00186 D( L ) = RT1 00187 D( L+1 ) = RT2 00188 E( L ) = ZERO 00189 L = L + 2 00190 IF( L.LE.LEND ) 00191 $ GO TO 50 00192 GO TO 150 00193 END IF 00194 * 00195 IF( JTOT.EQ.NMAXIT ) 00196 $ GO TO 150 00197 JTOT = JTOT + 1 00198 * 00199 * Form shift. 00200 * 00201 RTE = SQRT( E( L ) ) 00202 SIGMA = ( D( L+1 )-P ) / ( TWO*RTE ) 00203 R = DLAPY2( SIGMA, ONE ) 00204 SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) 00205 * 00206 C = ONE 00207 S = ZERO 00208 GAMMA = D( M ) - SIGMA 00209 P = GAMMA*GAMMA 00210 * 00211 * Inner loop 00212 * 00213 DO 80 I = M - 1, L, -1 00214 BB = E( I ) 00215 R = P + BB 00216 IF( I.NE.M-1 ) 00217 $ E( I+1 ) = S*R 00218 OLDC = C 00219 C = P / R 00220 S = BB / R 00221 OLDGAM = GAMMA 00222 ALPHA = D( I ) 00223 GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM 00224 D( I+1 ) = OLDGAM + ( ALPHA-GAMMA ) 00225 IF( C.NE.ZERO ) THEN 00226 P = ( GAMMA*GAMMA ) / C 00227 ELSE 00228 P = OLDC*BB 00229 END IF 00230 80 CONTINUE 00231 * 00232 E( L ) = S*P 00233 D( L ) = SIGMA + GAMMA 00234 GO TO 50 00235 * 00236 * Eigenvalue found. 00237 * 00238 90 CONTINUE 00239 D( L ) = P 00240 * 00241 L = L + 1 00242 IF( L.LE.LEND ) 00243 $ GO TO 50 00244 GO TO 150 00245 * 00246 ELSE 00247 * 00248 * QR Iteration 00249 * 00250 * Look for small superdiagonal element. 00251 * 00252 100 CONTINUE 00253 DO 110 M = L, LEND + 1, -1 00254 IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) ) 00255 $ GO TO 120 00256 110 CONTINUE 00257 M = LEND 00258 * 00259 120 CONTINUE 00260 IF( M.GT.LEND ) 00261 $ E( M-1 ) = ZERO 00262 P = D( L ) 00263 IF( M.EQ.L ) 00264 $ GO TO 140 00265 * 00266 * If remaining matrix is 2 by 2, use DLAE2 to compute its 00267 * eigenvalues. 00268 * 00269 IF( M.EQ.L-1 ) THEN 00270 RTE = SQRT( E( L-1 ) ) 00271 CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 ) 00272 D( L ) = RT1 00273 D( L-1 ) = RT2 00274 E( L-1 ) = ZERO 00275 L = L - 2 00276 IF( L.GE.LEND ) 00277 $ GO TO 100 00278 GO TO 150 00279 END IF 00280 * 00281 IF( JTOT.EQ.NMAXIT ) 00282 $ GO TO 150 00283 JTOT = JTOT + 1 00284 * 00285 * Form shift. 00286 * 00287 RTE = SQRT( E( L-1 ) ) 00288 SIGMA = ( D( L-1 )-P ) / ( TWO*RTE ) 00289 R = DLAPY2( SIGMA, ONE ) 00290 SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) 00291 * 00292 C = ONE 00293 S = ZERO 00294 GAMMA = D( M ) - SIGMA 00295 P = GAMMA*GAMMA 00296 * 00297 * Inner loop 00298 * 00299 DO 130 I = M, L - 1 00300 BB = E( I ) 00301 R = P + BB 00302 IF( I.NE.M ) 00303 $ E( I-1 ) = S*R 00304 OLDC = C 00305 C = P / R 00306 S = BB / R 00307 OLDGAM = GAMMA 00308 ALPHA = D( I+1 ) 00309 GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM 00310 D( I ) = OLDGAM + ( ALPHA-GAMMA ) 00311 IF( C.NE.ZERO ) THEN 00312 P = ( GAMMA*GAMMA ) / C 00313 ELSE 00314 P = OLDC*BB 00315 END IF 00316 130 CONTINUE 00317 * 00318 E( L-1 ) = S*P 00319 D( L ) = SIGMA + GAMMA 00320 GO TO 100 00321 * 00322 * Eigenvalue found. 00323 * 00324 140 CONTINUE 00325 D( L ) = P 00326 * 00327 L = L - 1 00328 IF( L.GE.LEND ) 00329 $ GO TO 100 00330 GO TO 150 00331 * 00332 END IF 00333 * 00334 * Undo scaling if necessary 00335 * 00336 150 CONTINUE 00337 IF( ISCALE.EQ.1 ) 00338 $ CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, 00339 $ D( LSV ), N, INFO ) 00340 IF( ISCALE.EQ.2 ) 00341 $ CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, 00342 $ D( LSV ), N, INFO ) 00343 * 00344 * Check for no convergence to an eigenvalue after a total 00345 * of N*MAXIT iterations. 00346 * 00347 IF( JTOT.LT.NMAXIT ) 00348 $ GO TO 10 00349 DO 160 I = 1, N - 1 00350 IF( E( I ).NE.ZERO ) 00351 $ INFO = INFO + 1 00352 160 CONTINUE 00353 GO TO 180 00354 * 00355 * Sort eigenvalues in increasing order. 00356 * 00357 170 CONTINUE 00358 CALL DLASRT( 'I', N, D, INFO ) 00359 * 00360 180 CONTINUE 00361 RETURN 00362 * 00363 * End of DSTERF 00364 * 00365 END
1.7.3 
