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