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