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