![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SSTERF 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SSTERF + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssterf.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssterf.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssterf.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SSTERF( N, D, E, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER INFO, N 00025 * .. 00026 * .. Array Arguments .. 00027 * REAL D( * ), E( * ) 00028 * .. 00029 * 00030 * 00031 *> \par Purpose: 00032 * ============= 00033 *> 00034 *> \verbatim 00035 *> 00036 *> SSTERF computes all eigenvalues of a symmetric tridiagonal matrix 00037 *> using the Pal-Walker-Kahan variant of the QL or QR algorithm. 00038 *> \endverbatim 00039 * 00040 * Arguments: 00041 * ========== 00042 * 00043 *> \param[in] N 00044 *> \verbatim 00045 *> N is INTEGER 00046 *> The order of the matrix. N >= 0. 00047 *> \endverbatim 00048 *> 00049 *> \param[in,out] D 00050 *> \verbatim 00051 *> D is REAL array, dimension (N) 00052 *> On entry, the n diagonal elements of the tridiagonal matrix. 00053 *> On exit, if INFO = 0, the eigenvalues in ascending order. 00054 *> \endverbatim 00055 *> 00056 *> \param[in,out] E 00057 *> \verbatim 00058 *> E is REAL array, dimension (N-1) 00059 *> On entry, the (n-1) subdiagonal elements of the tridiagonal 00060 *> matrix. 00061 *> On exit, E has been destroyed. 00062 *> \endverbatim 00063 *> 00064 *> \param[out] INFO 00065 *> \verbatim 00066 *> INFO is INTEGER 00067 *> = 0: successful exit 00068 *> < 0: if INFO = -i, the i-th argument had an illegal value 00069 *> > 0: the algorithm failed to find all of the eigenvalues in 00070 *> a total of 30*N iterations; if INFO = i, then i 00071 *> elements of E have not converged to zero. 00072 *> \endverbatim 00073 * 00074 * Authors: 00075 * ======== 00076 * 00077 *> \author Univ. of Tennessee 00078 *> \author Univ. of California Berkeley 00079 *> \author Univ. of Colorado Denver 00080 *> \author NAG Ltd. 00081 * 00082 *> \date November 2011 00083 * 00084 *> \ingroup auxOTHERcomputational 00085 * 00086 * ===================================================================== 00087 SUBROUTINE SSTERF( N, D, E, INFO ) 00088 * 00089 * -- LAPACK computational routine (version 3.4.0) -- 00090 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00091 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00092 * November 2011 00093 * 00094 * .. Scalar Arguments .. 00095 INTEGER INFO, N 00096 * .. 00097 * .. Array Arguments .. 00098 REAL D( * ), E( * ) 00099 * .. 00100 * 00101 * ===================================================================== 00102 * 00103 * .. Parameters .. 00104 REAL ZERO, ONE, TWO, THREE 00105 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, 00106 $ THREE = 3.0E0 ) 00107 INTEGER MAXIT 00108 PARAMETER ( MAXIT = 30 ) 00109 * .. 00110 * .. Local Scalars .. 00111 INTEGER I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M, 00112 $ NMAXIT 00113 REAL ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC, 00114 $ OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN, 00115 $ SIGMA, SSFMAX, SSFMIN 00116 * .. 00117 * .. External Functions .. 00118 REAL SLAMCH, SLANST, SLAPY2 00119 EXTERNAL SLAMCH, SLANST, SLAPY2 00120 * .. 00121 * .. External Subroutines .. 00122 EXTERNAL SLAE2, SLASCL, SLASRT, XERBLA 00123 * .. 00124 * .. Intrinsic Functions .. 00125 INTRINSIC ABS, SIGN, SQRT 00126 * .. 00127 * .. Executable Statements .. 00128 * 00129 * Test the input parameters. 00130 * 00131 INFO = 0 00132 * 00133 * Quick return if possible 00134 * 00135 IF( N.LT.0 ) THEN 00136 INFO = -1 00137 CALL XERBLA( 'SSTERF', -INFO ) 00138 RETURN 00139 END IF 00140 IF( N.LE.1 ) 00141 $ RETURN 00142 * 00143 * Determine the unit roundoff for this environment. 00144 * 00145 EPS = SLAMCH( 'E' ) 00146 EPS2 = EPS**2 00147 SAFMIN = SLAMCH( 'S' ) 00148 SAFMAX = ONE / SAFMIN 00149 SSFMAX = SQRT( SAFMAX ) / THREE 00150 SSFMIN = SQRT( SAFMIN ) / EPS2 00151 * 00152 * Compute the eigenvalues of the tridiagonal matrix. 00153 * 00154 NMAXIT = N*MAXIT 00155 SIGMA = ZERO 00156 JTOT = 0 00157 * 00158 * Determine where the matrix splits and choose QL or QR iteration 00159 * for each block, according to whether top or bottom diagonal 00160 * element is smaller. 00161 * 00162 L1 = 1 00163 * 00164 10 CONTINUE 00165 IF( L1.GT.N ) 00166 $ GO TO 170 00167 IF( L1.GT.1 ) 00168 $ E( L1-1 ) = ZERO 00169 DO 20 M = L1, N - 1 00170 IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )* 00171 $ SQRT( ABS( D( M+1 ) ) ) )*EPS ) THEN 00172 E( M ) = ZERO 00173 GO TO 30 00174 END IF 00175 20 CONTINUE 00176 M = N 00177 * 00178 30 CONTINUE 00179 L = L1 00180 LSV = L 00181 LEND = M 00182 LENDSV = LEND 00183 L1 = M + 1 00184 IF( LEND.EQ.L ) 00185 $ GO TO 10 00186 * 00187 * Scale submatrix in rows and columns L to LEND 00188 * 00189 ANORM = SLANST( 'M', LEND-L+1, D( L ), E( L ) ) 00190 ISCALE = 0 00191 IF( ANORM.EQ.ZERO ) 00192 $ GO TO 10 00193 IF( ANORM.GT.SSFMAX ) THEN 00194 ISCALE = 1 00195 CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, 00196 $ INFO ) 00197 CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, 00198 $ INFO ) 00199 ELSE IF( ANORM.LT.SSFMIN ) THEN 00200 ISCALE = 2 00201 CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, 00202 $ INFO ) 00203 CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, 00204 $ INFO ) 00205 END IF 00206 * 00207 DO 40 I = L, LEND - 1 00208 E( I ) = E( I )**2 00209 40 CONTINUE 00210 * 00211 * Choose between QL and QR iteration 00212 * 00213 IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN 00214 LEND = LSV 00215 L = LENDSV 00216 END IF 00217 * 00218 IF( LEND.GE.L ) THEN 00219 * 00220 * QL Iteration 00221 * 00222 * Look for small subdiagonal element. 00223 * 00224 50 CONTINUE 00225 IF( L.NE.LEND ) THEN 00226 DO 60 M = L, LEND - 1 00227 IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) ) 00228 $ GO TO 70 00229 60 CONTINUE 00230 END IF 00231 M = LEND 00232 * 00233 70 CONTINUE 00234 IF( M.LT.LEND ) 00235 $ E( M ) = ZERO 00236 P = D( L ) 00237 IF( M.EQ.L ) 00238 $ GO TO 90 00239 * 00240 * If remaining matrix is 2 by 2, use SLAE2 to compute its 00241 * eigenvalues. 00242 * 00243 IF( M.EQ.L+1 ) THEN 00244 RTE = SQRT( E( L ) ) 00245 CALL SLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 ) 00246 D( L ) = RT1 00247 D( L+1 ) = RT2 00248 E( L ) = ZERO 00249 L = L + 2 00250 IF( L.LE.LEND ) 00251 $ GO TO 50 00252 GO TO 150 00253 END IF 00254 * 00255 IF( JTOT.EQ.NMAXIT ) 00256 $ GO TO 150 00257 JTOT = JTOT + 1 00258 * 00259 * Form shift. 00260 * 00261 RTE = SQRT( E( L ) ) 00262 SIGMA = ( D( L+1 )-P ) / ( TWO*RTE ) 00263 R = SLAPY2( SIGMA, ONE ) 00264 SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) 00265 * 00266 C = ONE 00267 S = ZERO 00268 GAMMA = D( M ) - SIGMA 00269 P = GAMMA*GAMMA 00270 * 00271 * Inner loop 00272 * 00273 DO 80 I = M - 1, L, -1 00274 BB = E( I ) 00275 R = P + BB 00276 IF( I.NE.M-1 ) 00277 $ E( I+1 ) = S*R 00278 OLDC = C 00279 C = P / R 00280 S = BB / R 00281 OLDGAM = GAMMA 00282 ALPHA = D( I ) 00283 GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM 00284 D( I+1 ) = OLDGAM + ( ALPHA-GAMMA ) 00285 IF( C.NE.ZERO ) THEN 00286 P = ( GAMMA*GAMMA ) / C 00287 ELSE 00288 P = OLDC*BB 00289 END IF 00290 80 CONTINUE 00291 * 00292 E( L ) = S*P 00293 D( L ) = SIGMA + GAMMA 00294 GO TO 50 00295 * 00296 * Eigenvalue found. 00297 * 00298 90 CONTINUE 00299 D( L ) = P 00300 * 00301 L = L + 1 00302 IF( L.LE.LEND ) 00303 $ GO TO 50 00304 GO TO 150 00305 * 00306 ELSE 00307 * 00308 * QR Iteration 00309 * 00310 * Look for small superdiagonal element. 00311 * 00312 100 CONTINUE 00313 DO 110 M = L, LEND + 1, -1 00314 IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) ) 00315 $ GO TO 120 00316 110 CONTINUE 00317 M = LEND 00318 * 00319 120 CONTINUE 00320 IF( M.GT.LEND ) 00321 $ E( M-1 ) = ZERO 00322 P = D( L ) 00323 IF( M.EQ.L ) 00324 $ GO TO 140 00325 * 00326 * If remaining matrix is 2 by 2, use SLAE2 to compute its 00327 * eigenvalues. 00328 * 00329 IF( M.EQ.L-1 ) THEN 00330 RTE = SQRT( E( L-1 ) ) 00331 CALL SLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 ) 00332 D( L ) = RT1 00333 D( L-1 ) = RT2 00334 E( L-1 ) = ZERO 00335 L = L - 2 00336 IF( L.GE.LEND ) 00337 $ GO TO 100 00338 GO TO 150 00339 END IF 00340 * 00341 IF( JTOT.EQ.NMAXIT ) 00342 $ GO TO 150 00343 JTOT = JTOT + 1 00344 * 00345 * Form shift. 00346 * 00347 RTE = SQRT( E( L-1 ) ) 00348 SIGMA = ( D( L-1 )-P ) / ( TWO*RTE ) 00349 R = SLAPY2( SIGMA, ONE ) 00350 SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) 00351 * 00352 C = ONE 00353 S = ZERO 00354 GAMMA = D( M ) - SIGMA 00355 P = GAMMA*GAMMA 00356 * 00357 * Inner loop 00358 * 00359 DO 130 I = M, L - 1 00360 BB = E( I ) 00361 R = P + BB 00362 IF( I.NE.M ) 00363 $ E( I-1 ) = S*R 00364 OLDC = C 00365 C = P / R 00366 S = BB / R 00367 OLDGAM = GAMMA 00368 ALPHA = D( I+1 ) 00369 GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM 00370 D( I ) = OLDGAM + ( ALPHA-GAMMA ) 00371 IF( C.NE.ZERO ) THEN 00372 P = ( GAMMA*GAMMA ) / C 00373 ELSE 00374 P = OLDC*BB 00375 END IF 00376 130 CONTINUE 00377 * 00378 E( L-1 ) = S*P 00379 D( L ) = SIGMA + GAMMA 00380 GO TO 100 00381 * 00382 * Eigenvalue found. 00383 * 00384 140 CONTINUE 00385 D( L ) = P 00386 * 00387 L = L - 1 00388 IF( L.GE.LEND ) 00389 $ GO TO 100 00390 GO TO 150 00391 * 00392 END IF 00393 * 00394 * Undo scaling if necessary 00395 * 00396 150 CONTINUE 00397 IF( ISCALE.EQ.1 ) 00398 $ CALL SLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, 00399 $ D( LSV ), N, INFO ) 00400 IF( ISCALE.EQ.2 ) 00401 $ CALL SLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, 00402 $ D( LSV ), N, INFO ) 00403 * 00404 * Check for no convergence to an eigenvalue after a total 00405 * of N*MAXIT iterations. 00406 * 00407 IF( JTOT.LT.NMAXIT ) 00408 $ GO TO 10 00409 DO 160 I = 1, N - 1 00410 IF( E( I ).NE.ZERO ) 00411 $ INFO = INFO + 1 00412 160 CONTINUE 00413 GO TO 180 00414 * 00415 * Sort eigenvalues in increasing order. 00416 * 00417 170 CONTINUE 00418 CALL SLASRT( 'I', N, D, INFO ) 00419 * 00420 180 CONTINUE 00421 RETURN 00422 * 00423 * End of SSTERF 00424 * 00425 END