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