 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DSTERF( N, D, E, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.3.1) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * -- April 2011 -- 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, RMAX 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 RMAX = DLAMCH( 'O' ) 00094 * 00095 * Compute the eigenvalues of the tridiagonal matrix. 00096 * 00097 NMAXIT = N*MAXIT 00098 SIGMA = ZERO 00099 JTOT = 0 00100 * 00101 * Determine where the matrix splits and choose QL or QR iteration 00102 * for each block, according to whether top or bottom diagonal 00103 * element is smaller. 00104 * 00105 L1 = 1 00106 * 00107 10 CONTINUE 00108 IF( L1.GT.N ) 00109 $ GO TO 170 00110 IF( L1.GT.1 ) 00111 $ E( L1-1 ) = ZERO 00112 DO 20 M = L1, N - 1 00113 IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ 00114 $ 1 ) ) ) )*EPS ) THEN 00115 E( M ) = ZERO 00116 GO TO 30 00117 END IF 00118 20 CONTINUE 00119 M = N 00120 * 00121 30 CONTINUE 00122 L = L1 00123 LSV = L 00124 LEND = M 00125 LENDSV = LEND 00126 L1 = M + 1 00127 IF( LEND.EQ.L ) 00128 $ GO TO 10 00129 * 00130 * Scale submatrix in rows and columns L to LEND 00131 * 00132 ANORM = DLANST( 'M', LEND-L+1, D( L ), E( L ) ) 00133 ISCALE = 0 00134 IF( ANORM.EQ.ZERO ) 00135 $ GO TO 10 00136 IF( (ANORM.GT.SSFMAX) ) THEN 00137 ISCALE = 1 00138 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, 00139 $ INFO ) 00140 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, 00141 $ INFO ) 00142 ELSE IF( ANORM.LT.SSFMIN ) THEN 00143 ISCALE = 2 00144 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, 00145 $ INFO ) 00146 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, 00147 $ INFO ) 00148 END IF 00149 * 00150 DO 40 I = L, LEND - 1 00151 E( I ) = E( I )**2 00152 40 CONTINUE 00153 * 00154 * Choose between QL and QR iteration 00155 * 00156 IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN 00157 LEND = LSV 00158 L = LENDSV 00159 END IF 00160 * 00161 IF( LEND.GE.L ) THEN 00162 * 00163 * QL Iteration 00164 * 00165 * Look for small subdiagonal element. 00166 * 00167 50 CONTINUE 00168 IF( L.NE.LEND ) THEN 00169 DO 60 M = L, LEND - 1 00170 IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) ) 00171 $ GO TO 70 00172 60 CONTINUE 00173 END IF 00174 M = LEND 00175 * 00176 70 CONTINUE 00177 IF( M.LT.LEND ) 00178 $ E( M ) = ZERO 00179 P = D( L ) 00180 IF( M.EQ.L ) 00181 $ GO TO 90 00182 * 00183 * If remaining matrix is 2 by 2, use DLAE2 to compute its 00184 * eigenvalues. 00185 * 00186 IF( M.EQ.L+1 ) THEN 00187 RTE = SQRT( E( L ) ) 00188 CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 ) 00189 D( L ) = RT1 00190 D( L+1 ) = RT2 00191 E( L ) = ZERO 00192 L = L + 2 00193 IF( L.LE.LEND ) 00194 $ GO TO 50 00195 GO TO 150 00196 END IF 00197 * 00198 IF( JTOT.EQ.NMAXIT ) 00199 $ GO TO 150 00200 JTOT = JTOT + 1 00201 * 00202 * Form shift. 00203 * 00204 RTE = SQRT( E( L ) ) 00205 SIGMA = ( D( L+1 )-P ) / ( TWO*RTE ) 00206 R = DLAPY2( SIGMA, ONE ) 00207 SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) 00208 * 00209 C = ONE 00210 S = ZERO 00211 GAMMA = D( M ) - SIGMA 00212 P = GAMMA*GAMMA 00213 * 00214 * Inner loop 00215 * 00216 DO 80 I = M - 1, L, -1 00217 BB = E( I ) 00218 R = P + BB 00219 IF( I.NE.M-1 ) 00220 $ E( I+1 ) = S*R 00221 OLDC = C 00222 C = P / R 00223 S = BB / R 00224 OLDGAM = GAMMA 00225 ALPHA = D( I ) 00226 GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM 00227 D( I+1 ) = OLDGAM + ( ALPHA-GAMMA ) 00228 IF( C.NE.ZERO ) THEN 00229 P = ( GAMMA*GAMMA ) / C 00230 ELSE 00231 P = OLDC*BB 00232 END IF 00233 80 CONTINUE 00234 * 00235 E( L ) = S*P 00236 D( L ) = SIGMA + GAMMA 00237 GO TO 50 00238 * 00239 * Eigenvalue found. 00240 * 00241 90 CONTINUE 00242 D( L ) = P 00243 * 00244 L = L + 1 00245 IF( L.LE.LEND ) 00246 $ GO TO 50 00247 GO TO 150 00248 * 00249 ELSE 00250 * 00251 * QR Iteration 00252 * 00253 * Look for small superdiagonal element. 00254 * 00255 100 CONTINUE 00256 DO 110 M = L, LEND + 1, -1 00257 IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) ) 00258 $ GO TO 120 00259 110 CONTINUE 00260 M = LEND 00261 * 00262 120 CONTINUE 00263 IF( M.GT.LEND ) 00264 $ E( M-1 ) = ZERO 00265 P = D( L ) 00266 IF( M.EQ.L ) 00267 $ GO TO 140 00268 * 00269 * If remaining matrix is 2 by 2, use DLAE2 to compute its 00270 * eigenvalues. 00271 * 00272 IF( M.EQ.L-1 ) THEN 00273 RTE = SQRT( E( L-1 ) ) 00274 CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 ) 00275 D( L ) = RT1 00276 D( L-1 ) = RT2 00277 E( L-1 ) = ZERO 00278 L = L - 2 00279 IF( L.GE.LEND ) 00280 $ GO TO 100 00281 GO TO 150 00282 END IF 00283 * 00284 IF( JTOT.EQ.NMAXIT ) 00285 $ GO TO 150 00286 JTOT = JTOT + 1 00287 * 00288 * Form shift. 00289 * 00290 RTE = SQRT( E( L-1 ) ) 00291 SIGMA = ( D( L-1 )-P ) / ( TWO*RTE ) 00292 R = DLAPY2( SIGMA, ONE ) 00293 SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) 00294 * 00295 C = ONE 00296 S = ZERO 00297 GAMMA = D( M ) - SIGMA 00298 P = GAMMA*GAMMA 00299 * 00300 * Inner loop 00301 * 00302 DO 130 I = M, L - 1 00303 BB = E( I ) 00304 R = P + BB 00305 IF( I.NE.M ) 00306 $ E( I-1 ) = S*R 00307 OLDC = C 00308 C = P / R 00309 S = BB / R 00310 OLDGAM = GAMMA 00311 ALPHA = D( I+1 ) 00312 GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM 00313 D( I ) = OLDGAM + ( ALPHA-GAMMA ) 00314 IF( C.NE.ZERO ) THEN 00315 P = ( GAMMA*GAMMA ) / C 00316 ELSE 00317 P = OLDC*BB 00318 END IF 00319 130 CONTINUE 00320 * 00321 E( L-1 ) = S*P 00322 D( L ) = SIGMA + GAMMA 00323 GO TO 100 00324 * 00325 * Eigenvalue found. 00326 * 00327 140 CONTINUE 00328 D( L ) = P 00329 * 00330 L = L - 1 00331 IF( L.GE.LEND ) 00332 $ GO TO 100 00333 GO TO 150 00334 * 00335 END IF 00336 * 00337 * Undo scaling if necessary 00338 * 00339 150 CONTINUE 00340 IF( ISCALE.EQ.1 ) 00341 $ CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, 00342 $ D( LSV ), N, INFO ) 00343 IF( ISCALE.EQ.2 ) 00344 $ CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, 00345 $ D( LSV ), N, INFO ) 00346 * 00347 * Check for no convergence to an eigenvalue after a total 00348 * of N*MAXIT iterations. 00349 * 00350 IF( JTOT.LT.NMAXIT ) 00351 $ GO TO 10 00352 DO 160 I = 1, N - 1 00353 IF( E( I ).NE.ZERO ) 00354 $ INFO = INFO + 1 00355 160 CONTINUE 00356 GO TO 180 00357 * 00358 * Sort eigenvalues in increasing order. 00359 * 00360 170 CONTINUE 00361 CALL DLASRT( 'I', N, D, INFO ) 00362 * 00363 180 CONTINUE 00364 RETURN 00365 * 00366 * End of DSTERF 00367 * 00368 END
1.7.3 
