 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, 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 CHARACTER COMPZ 00010 INTEGER INFO, LDZ, N 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * DSTEQR computes all eigenvalues and, optionally, eigenvectors of a 00020 * symmetric tridiagonal matrix using the implicit QL or QR method. 00021 * The eigenvectors of a full or band symmetric matrix can also be found 00022 * if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to 00023 * tridiagonal form. 00024 * 00025 * Arguments 00026 * ========= 00027 * 00028 * COMPZ (input) CHARACTER*1 00029 * = 'N': Compute eigenvalues only. 00030 * = 'V': Compute eigenvalues and eigenvectors of the original 00031 * symmetric matrix. On entry, Z must contain the 00032 * orthogonal matrix used to reduce the original matrix 00033 * to tridiagonal form. 00034 * = 'I': Compute eigenvalues and eigenvectors of the 00035 * tridiagonal matrix. Z is initialized to the identity 00036 * matrix. 00037 * 00038 * N (input) INTEGER 00039 * The order of the matrix. N >= 0. 00040 * 00041 * D (input/output) DOUBLE PRECISION array, dimension (N) 00042 * On entry, the diagonal elements of the tridiagonal matrix. 00043 * On exit, if INFO = 0, the eigenvalues in ascending order. 00044 * 00045 * E (input/output) DOUBLE PRECISION array, dimension (N-1) 00046 * On entry, the (n-1) subdiagonal elements of the tridiagonal 00047 * matrix. 00048 * On exit, E has been destroyed. 00049 * 00050 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N) 00051 * On entry, if COMPZ = 'V', then Z contains the orthogonal 00052 * matrix used in the reduction to tridiagonal form. 00053 * On exit, if INFO = 0, then if COMPZ = 'V', Z contains the 00054 * orthonormal eigenvectors of the original symmetric matrix, 00055 * and if COMPZ = 'I', Z contains the orthonormal eigenvectors 00056 * of the symmetric tridiagonal matrix. 00057 * If COMPZ = 'N', then Z is not referenced. 00058 * 00059 * LDZ (input) INTEGER 00060 * The leading dimension of the array Z. LDZ >= 1, and if 00061 * eigenvectors are desired, then LDZ >= max(1,N). 00062 * 00063 * WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2)) 00064 * If COMPZ = 'N', then WORK is not referenced. 00065 * 00066 * INFO (output) INTEGER 00067 * = 0: successful exit 00068 * < 0: if INFO = -i, the i-th argument had an illegal value 00069 * > 0: the algorithm has failed to find all the eigenvalues in 00070 * a total of 30*N iterations; if INFO = i, then i 00071 * elements of E have not converged to zero; on exit, D 00072 * and E contain the elements of a symmetric tridiagonal 00073 * matrix which is orthogonally similar to the original 00074 * matrix. 00075 * 00076 * ===================================================================== 00077 * 00078 * .. Parameters .. 00079 DOUBLE PRECISION ZERO, ONE, TWO, THREE 00080 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 00081 $ THREE = 3.0D0 ) 00082 INTEGER MAXIT 00083 PARAMETER ( MAXIT = 30 ) 00084 * .. 00085 * .. Local Scalars .. 00086 INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND, 00087 $ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1, 00088 $ NM1, NMAXIT 00089 DOUBLE PRECISION ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2, 00090 $ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST 00091 * .. 00092 * .. External Functions .. 00093 LOGICAL LSAME 00094 DOUBLE PRECISION DLAMCH, DLANST, DLAPY2 00095 EXTERNAL LSAME, DLAMCH, DLANST, DLAPY2 00096 * .. 00097 * .. External Subroutines .. 00098 EXTERNAL DLAE2, DLAEV2, DLARTG, DLASCL, DLASET, DLASR, 00099 $ DLASRT, DSWAP, XERBLA 00100 * .. 00101 * .. Intrinsic Functions .. 00102 INTRINSIC ABS, MAX, SIGN, SQRT 00103 * .. 00104 * .. Executable Statements .. 00105 * 00106 * Test the input parameters. 00107 * 00108 INFO = 0 00109 * 00110 IF( LSAME( COMPZ, 'N' ) ) THEN 00111 ICOMPZ = 0 00112 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 00113 ICOMPZ = 1 00114 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 00115 ICOMPZ = 2 00116 ELSE 00117 ICOMPZ = -1 00118 END IF 00119 IF( ICOMPZ.LT.0 ) THEN 00120 INFO = -1 00121 ELSE IF( N.LT.0 ) THEN 00122 INFO = -2 00123 ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, 00124 $ N ) ) ) THEN 00125 INFO = -6 00126 END IF 00127 IF( INFO.NE.0 ) THEN 00128 CALL XERBLA( 'DSTEQR', -INFO ) 00129 RETURN 00130 END IF 00131 * 00132 * Quick return if possible 00133 * 00134 IF( N.EQ.0 ) 00135 $ RETURN 00136 * 00137 IF( N.EQ.1 ) THEN 00138 IF( ICOMPZ.EQ.2 ) 00139 $ Z( 1, 1 ) = ONE 00140 RETURN 00141 END IF 00142 * 00143 * Determine the unit roundoff and over/underflow thresholds. 00144 * 00145 EPS = DLAMCH( 'E' ) 00146 EPS2 = EPS**2 00147 SAFMIN = DLAMCH( 'S' ) 00148 SAFMAX = ONE / SAFMIN 00149 SSFMAX = SQRT( SAFMAX ) / THREE 00150 SSFMIN = SQRT( SAFMIN ) / EPS2 00151 * 00152 * Compute the eigenvalues and eigenvectors of the tridiagonal 00153 * matrix. 00154 * 00155 IF( ICOMPZ.EQ.2 ) 00156 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ ) 00157 * 00158 NMAXIT = N*MAXIT 00159 JTOT = 0 00160 * 00161 * Determine where the matrix splits and choose QL or QR iteration 00162 * for each block, according to whether top or bottom diagonal 00163 * element is smaller. 00164 * 00165 L1 = 1 00166 NM1 = N - 1 00167 * 00168 10 CONTINUE 00169 IF( L1.GT.N ) 00170 $ GO TO 160 00171 IF( L1.GT.1 ) 00172 $ E( L1-1 ) = ZERO 00173 IF( L1.LE.NM1 ) THEN 00174 DO 20 M = L1, NM1 00175 TST = ABS( E( M ) ) 00176 IF( TST.EQ.ZERO ) 00177 $ GO TO 30 00178 IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ 00179 $ 1 ) ) ) )*EPS ) THEN 00180 E( M ) = ZERO 00181 GO TO 30 00182 END IF 00183 20 CONTINUE 00184 END IF 00185 M = N 00186 * 00187 30 CONTINUE 00188 L = L1 00189 LSV = L 00190 LEND = M 00191 LENDSV = LEND 00192 L1 = M + 1 00193 IF( LEND.EQ.L ) 00194 $ GO TO 10 00195 * 00196 * Scale submatrix in rows and columns L to LEND 00197 * 00198 ANORM = DLANST( 'M', LEND-L+1, D( L ), E( L ) ) 00199 ISCALE = 0 00200 IF( ANORM.EQ.ZERO ) 00201 $ GO TO 10 00202 IF( ANORM.GT.SSFMAX ) THEN 00203 ISCALE = 1 00204 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, 00205 $ INFO ) 00206 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, 00207 $ INFO ) 00208 ELSE IF( ANORM.LT.SSFMIN ) THEN 00209 ISCALE = 2 00210 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, 00211 $ INFO ) 00212 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, 00213 $ INFO ) 00214 END IF 00215 * 00216 * Choose between QL and QR iteration 00217 * 00218 IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN 00219 LEND = LSV 00220 L = LENDSV 00221 END IF 00222 * 00223 IF( LEND.GT.L ) THEN 00224 * 00225 * QL Iteration 00226 * 00227 * Look for small subdiagonal element. 00228 * 00229 40 CONTINUE 00230 IF( L.NE.LEND ) THEN 00231 LENDM1 = LEND - 1 00232 DO 50 M = L, LENDM1 00233 TST = ABS( E( M ) )**2 00234 IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+ 00235 $ SAFMIN )GO TO 60 00236 50 CONTINUE 00237 END IF 00238 * 00239 M = LEND 00240 * 00241 60 CONTINUE 00242 IF( M.LT.LEND ) 00243 $ E( M ) = ZERO 00244 P = D( L ) 00245 IF( M.EQ.L ) 00246 $ GO TO 80 00247 * 00248 * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 00249 * to compute its eigensystem. 00250 * 00251 IF( M.EQ.L+1 ) THEN 00252 IF( ICOMPZ.GT.0 ) THEN 00253 CALL DLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S ) 00254 WORK( L ) = C 00255 WORK( N-1+L ) = S 00256 CALL DLASR( 'R', 'V', 'B', N, 2, WORK( L ), 00257 $ WORK( N-1+L ), Z( 1, L ), LDZ ) 00258 ELSE 00259 CALL DLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 ) 00260 END IF 00261 D( L ) = RT1 00262 D( L+1 ) = RT2 00263 E( L ) = ZERO 00264 L = L + 2 00265 IF( L.LE.LEND ) 00266 $ GO TO 40 00267 GO TO 140 00268 END IF 00269 * 00270 IF( JTOT.EQ.NMAXIT ) 00271 $ GO TO 140 00272 JTOT = JTOT + 1 00273 * 00274 * Form shift. 00275 * 00276 G = ( D( L+1 )-P ) / ( TWO*E( L ) ) 00277 R = DLAPY2( G, ONE ) 00278 G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) ) 00279 * 00280 S = ONE 00281 C = ONE 00282 P = ZERO 00283 * 00284 * Inner loop 00285 * 00286 MM1 = M - 1 00287 DO 70 I = MM1, L, -1 00288 F = S*E( I ) 00289 B = C*E( I ) 00290 CALL DLARTG( G, F, C, S, R ) 00291 IF( I.NE.M-1 ) 00292 $ E( I+1 ) = R 00293 G = D( I+1 ) - P 00294 R = ( D( I )-G )*S + TWO*C*B 00295 P = S*R 00296 D( I+1 ) = G + P 00297 G = C*R - B 00298 * 00299 * If eigenvectors are desired, then save rotations. 00300 * 00301 IF( ICOMPZ.GT.0 ) THEN 00302 WORK( I ) = C 00303 WORK( N-1+I ) = -S 00304 END IF 00305 * 00306 70 CONTINUE 00307 * 00308 * If eigenvectors are desired, then apply saved rotations. 00309 * 00310 IF( ICOMPZ.GT.0 ) THEN 00311 MM = M - L + 1 00312 CALL DLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ), 00313 $ Z( 1, L ), LDZ ) 00314 END IF 00315 * 00316 D( L ) = D( L ) - P 00317 E( L ) = G 00318 GO TO 40 00319 * 00320 * Eigenvalue found. 00321 * 00322 80 CONTINUE 00323 D( L ) = P 00324 * 00325 L = L + 1 00326 IF( L.LE.LEND ) 00327 $ GO TO 40 00328 GO TO 140 00329 * 00330 ELSE 00331 * 00332 * QR Iteration 00333 * 00334 * Look for small superdiagonal element. 00335 * 00336 90 CONTINUE 00337 IF( L.NE.LEND ) THEN 00338 LENDP1 = LEND + 1 00339 DO 100 M = L, LENDP1, -1 00340 TST = ABS( E( M-1 ) )**2 00341 IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+ 00342 $ SAFMIN )GO TO 110 00343 100 CONTINUE 00344 END IF 00345 * 00346 M = LEND 00347 * 00348 110 CONTINUE 00349 IF( M.GT.LEND ) 00350 $ E( M-1 ) = ZERO 00351 P = D( L ) 00352 IF( M.EQ.L ) 00353 $ GO TO 130 00354 * 00355 * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 00356 * to compute its eigensystem. 00357 * 00358 IF( M.EQ.L-1 ) THEN 00359 IF( ICOMPZ.GT.0 ) THEN 00360 CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S ) 00361 WORK( M ) = C 00362 WORK( N-1+M ) = S 00363 CALL DLASR( 'R', 'V', 'F', N, 2, WORK( M ), 00364 $ WORK( N-1+M ), Z( 1, L-1 ), LDZ ) 00365 ELSE 00366 CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 ) 00367 END IF 00368 D( L-1 ) = RT1 00369 D( L ) = RT2 00370 E( L-1 ) = ZERO 00371 L = L - 2 00372 IF( L.GE.LEND ) 00373 $ GO TO 90 00374 GO TO 140 00375 END IF 00376 * 00377 IF( JTOT.EQ.NMAXIT ) 00378 $ GO TO 140 00379 JTOT = JTOT + 1 00380 * 00381 * Form shift. 00382 * 00383 G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) ) 00384 R = DLAPY2( G, ONE ) 00385 G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) ) 00386 * 00387 S = ONE 00388 C = ONE 00389 P = ZERO 00390 * 00391 * Inner loop 00392 * 00393 LM1 = L - 1 00394 DO 120 I = M, LM1 00395 F = S*E( I ) 00396 B = C*E( I ) 00397 CALL DLARTG( G, F, C, S, R ) 00398 IF( I.NE.M ) 00399 $ E( I-1 ) = R 00400 G = D( I ) - P 00401 R = ( D( I+1 )-G )*S + TWO*C*B 00402 P = S*R 00403 D( I ) = G + P 00404 G = C*R - B 00405 * 00406 * If eigenvectors are desired, then save rotations. 00407 * 00408 IF( ICOMPZ.GT.0 ) THEN 00409 WORK( I ) = C 00410 WORK( N-1+I ) = S 00411 END IF 00412 * 00413 120 CONTINUE 00414 * 00415 * If eigenvectors are desired, then apply saved rotations. 00416 * 00417 IF( ICOMPZ.GT.0 ) THEN 00418 MM = L - M + 1 00419 CALL DLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ), 00420 $ Z( 1, M ), LDZ ) 00421 END IF 00422 * 00423 D( L ) = D( L ) - P 00424 E( LM1 ) = G 00425 GO TO 90 00426 * 00427 * Eigenvalue found. 00428 * 00429 130 CONTINUE 00430 D( L ) = P 00431 * 00432 L = L - 1 00433 IF( L.GE.LEND ) 00434 $ GO TO 90 00435 GO TO 140 00436 * 00437 END IF 00438 * 00439 * Undo scaling if necessary 00440 * 00441 140 CONTINUE 00442 IF( ISCALE.EQ.1 ) THEN 00443 CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, 00444 $ D( LSV ), N, INFO ) 00445 CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ), 00446 $ N, INFO ) 00447 ELSE IF( ISCALE.EQ.2 ) THEN 00448 CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, 00449 $ D( LSV ), N, INFO ) 00450 CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ), 00451 $ N, INFO ) 00452 END IF 00453 * 00454 * Check for no convergence to an eigenvalue after a total 00455 * of N*MAXIT iterations. 00456 * 00457 IF( JTOT.LT.NMAXIT ) 00458 $ GO TO 10 00459 DO 150 I = 1, N - 1 00460 IF( E( I ).NE.ZERO ) 00461 $ INFO = INFO + 1 00462 150 CONTINUE 00463 GO TO 190 00464 * 00465 * Order eigenvalues and eigenvectors. 00466 * 00467 160 CONTINUE 00468 IF( ICOMPZ.EQ.0 ) THEN 00469 * 00470 * Use Quick Sort 00471 * 00472 CALL DLASRT( 'I', N, D, INFO ) 00473 * 00474 ELSE 00475 * 00476 * Use Selection Sort to minimize swaps of eigenvectors 00477 * 00478 DO 180 II = 2, N 00479 I = II - 1 00480 K = I 00481 P = D( I ) 00482 DO 170 J = II, N 00483 IF( D( J ).LT.P ) THEN 00484 K = J 00485 P = D( J ) 00486 END IF 00487 170 CONTINUE 00488 IF( K.NE.I ) THEN 00489 D( K ) = D( I ) 00490 D( I ) = P 00491 CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 ) 00492 END IF 00493 180 CONTINUE 00494 END IF 00495 * 00496 190 CONTINUE 00497 RETURN 00498 * 00499 * End of DSTEQR 00500 * 00501 END
1.7.3 
