![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DSTERF 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DSTERF + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsterf.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsterf.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsterf.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DSTERF( N, D, E, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER INFO, N 00025 * .. 00026 * .. Array Arguments .. 00027 * DOUBLE PRECISION D( * ), E( * ) 00028 * .. 00029 * 00030 * 00031 *> \par Purpose: 00032 * ============= 00033 *> 00034 *> \verbatim 00035 *> 00036 *> DSTERF 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DSTERF( 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 DOUBLE PRECISION D( * ), E( * ) 00099 * .. 00100 * 00101 * ===================================================================== 00102 * 00103 * .. Parameters .. 00104 DOUBLE PRECISION ZERO, ONE, TWO, THREE 00105 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 00106 $ THREE = 3.0D0 ) 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 DOUBLE PRECISION ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC, 00114 $ OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN, 00115 $ SIGMA, SSFMAX, SSFMIN, RMAX 00116 * .. 00117 * .. External Functions .. 00118 DOUBLE PRECISION DLAMCH, DLANST, DLAPY2 00119 EXTERNAL DLAMCH, DLANST, DLAPY2 00120 * .. 00121 * .. External Subroutines .. 00122 EXTERNAL DLAE2, DLASCL, DLASRT, 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( 'DSTERF', -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 = DLAMCH( 'E' ) 00146 EPS2 = EPS**2 00147 SAFMIN = DLAMCH( 'S' ) 00148 SAFMAX = ONE / SAFMIN 00149 SSFMAX = SQRT( SAFMAX ) / THREE 00150 SSFMIN = SQRT( SAFMIN ) / EPS2 00151 RMAX = DLAMCH( 'O' ) 00152 * 00153 * Compute the eigenvalues of the tridiagonal matrix. 00154 * 00155 NMAXIT = N*MAXIT 00156 SIGMA = ZERO 00157 JTOT = 0 00158 * 00159 * Determine where the matrix splits and choose QL or QR iteration 00160 * for each block, according to whether top or bottom diagonal 00161 * element is smaller. 00162 * 00163 L1 = 1 00164 * 00165 10 CONTINUE 00166 IF( L1.GT.N ) 00167 $ GO TO 170 00168 IF( L1.GT.1 ) 00169 $ E( L1-1 ) = ZERO 00170 DO 20 M = L1, N - 1 00171 IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ 00172 $ 1 ) ) ) )*EPS ) THEN 00173 E( M ) = ZERO 00174 GO TO 30 00175 END IF 00176 20 CONTINUE 00177 M = N 00178 * 00179 30 CONTINUE 00180 L = L1 00181 LSV = L 00182 LEND = M 00183 LENDSV = LEND 00184 L1 = M + 1 00185 IF( LEND.EQ.L ) 00186 $ GO TO 10 00187 * 00188 * Scale submatrix in rows and columns L to LEND 00189 * 00190 ANORM = DLANST( 'M', LEND-L+1, D( L ), E( L ) ) 00191 ISCALE = 0 00192 IF( ANORM.EQ.ZERO ) 00193 $ GO TO 10 00194 IF( (ANORM.GT.SSFMAX) ) THEN 00195 ISCALE = 1 00196 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, 00197 $ INFO ) 00198 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, 00199 $ INFO ) 00200 ELSE IF( ANORM.LT.SSFMIN ) THEN 00201 ISCALE = 2 00202 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, 00203 $ INFO ) 00204 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, 00205 $ INFO ) 00206 END IF 00207 * 00208 DO 40 I = L, LEND - 1 00209 E( I ) = E( I )**2 00210 40 CONTINUE 00211 * 00212 * Choose between QL and QR iteration 00213 * 00214 IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN 00215 LEND = LSV 00216 L = LENDSV 00217 END IF 00218 * 00219 IF( LEND.GE.L ) THEN 00220 * 00221 * QL Iteration 00222 * 00223 * Look for small subdiagonal element. 00224 * 00225 50 CONTINUE 00226 IF( L.NE.LEND ) THEN 00227 DO 60 M = L, LEND - 1 00228 IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) ) 00229 $ GO TO 70 00230 60 CONTINUE 00231 END IF 00232 M = LEND 00233 * 00234 70 CONTINUE 00235 IF( M.LT.LEND ) 00236 $ E( M ) = ZERO 00237 P = D( L ) 00238 IF( M.EQ.L ) 00239 $ GO TO 90 00240 * 00241 * If remaining matrix is 2 by 2, use DLAE2 to compute its 00242 * eigenvalues. 00243 * 00244 IF( M.EQ.L+1 ) THEN 00245 RTE = SQRT( E( L ) ) 00246 CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 ) 00247 D( L ) = RT1 00248 D( L+1 ) = RT2 00249 E( L ) = ZERO 00250 L = L + 2 00251 IF( L.LE.LEND ) 00252 $ GO TO 50 00253 GO TO 150 00254 END IF 00255 * 00256 IF( JTOT.EQ.NMAXIT ) 00257 $ GO TO 150 00258 JTOT = JTOT + 1 00259 * 00260 * Form shift. 00261 * 00262 RTE = SQRT( E( L ) ) 00263 SIGMA = ( D( L+1 )-P ) / ( TWO*RTE ) 00264 R = DLAPY2( SIGMA, ONE ) 00265 SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) 00266 * 00267 C = ONE 00268 S = ZERO 00269 GAMMA = D( M ) - SIGMA 00270 P = GAMMA*GAMMA 00271 * 00272 * Inner loop 00273 * 00274 DO 80 I = M - 1, L, -1 00275 BB = E( I ) 00276 R = P + BB 00277 IF( I.NE.M-1 ) 00278 $ E( I+1 ) = S*R 00279 OLDC = C 00280 C = P / R 00281 S = BB / R 00282 OLDGAM = GAMMA 00283 ALPHA = D( I ) 00284 GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM 00285 D( I+1 ) = OLDGAM + ( ALPHA-GAMMA ) 00286 IF( C.NE.ZERO ) THEN 00287 P = ( GAMMA*GAMMA ) / C 00288 ELSE 00289 P = OLDC*BB 00290 END IF 00291 80 CONTINUE 00292 * 00293 E( L ) = S*P 00294 D( L ) = SIGMA + GAMMA 00295 GO TO 50 00296 * 00297 * Eigenvalue found. 00298 * 00299 90 CONTINUE 00300 D( L ) = P 00301 * 00302 L = L + 1 00303 IF( L.LE.LEND ) 00304 $ GO TO 50 00305 GO TO 150 00306 * 00307 ELSE 00308 * 00309 * QR Iteration 00310 * 00311 * Look for small superdiagonal element. 00312 * 00313 100 CONTINUE 00314 DO 110 M = L, LEND + 1, -1 00315 IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) ) 00316 $ GO TO 120 00317 110 CONTINUE 00318 M = LEND 00319 * 00320 120 CONTINUE 00321 IF( M.GT.LEND ) 00322 $ E( M-1 ) = ZERO 00323 P = D( L ) 00324 IF( M.EQ.L ) 00325 $ GO TO 140 00326 * 00327 * If remaining matrix is 2 by 2, use DLAE2 to compute its 00328 * eigenvalues. 00329 * 00330 IF( M.EQ.L-1 ) THEN 00331 RTE = SQRT( E( L-1 ) ) 00332 CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 ) 00333 D( L ) = RT1 00334 D( L-1 ) = RT2 00335 E( L-1 ) = ZERO 00336 L = L - 2 00337 IF( L.GE.LEND ) 00338 $ GO TO 100 00339 GO TO 150 00340 END IF 00341 * 00342 IF( JTOT.EQ.NMAXIT ) 00343 $ GO TO 150 00344 JTOT = JTOT + 1 00345 * 00346 * Form shift. 00347 * 00348 RTE = SQRT( E( L-1 ) ) 00349 SIGMA = ( D( L-1 )-P ) / ( TWO*RTE ) 00350 R = DLAPY2( SIGMA, ONE ) 00351 SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) 00352 * 00353 C = ONE 00354 S = ZERO 00355 GAMMA = D( M ) - SIGMA 00356 P = GAMMA*GAMMA 00357 * 00358 * Inner loop 00359 * 00360 DO 130 I = M, L - 1 00361 BB = E( I ) 00362 R = P + BB 00363 IF( I.NE.M ) 00364 $ E( I-1 ) = S*R 00365 OLDC = C 00366 C = P / R 00367 S = BB / R 00368 OLDGAM = GAMMA 00369 ALPHA = D( I+1 ) 00370 GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM 00371 D( I ) = OLDGAM + ( ALPHA-GAMMA ) 00372 IF( C.NE.ZERO ) THEN 00373 P = ( GAMMA*GAMMA ) / C 00374 ELSE 00375 P = OLDC*BB 00376 END IF 00377 130 CONTINUE 00378 * 00379 E( L-1 ) = S*P 00380 D( L ) = SIGMA + GAMMA 00381 GO TO 100 00382 * 00383 * Eigenvalue found. 00384 * 00385 140 CONTINUE 00386 D( L ) = P 00387 * 00388 L = L - 1 00389 IF( L.GE.LEND ) 00390 $ GO TO 100 00391 GO TO 150 00392 * 00393 END IF 00394 * 00395 * Undo scaling if necessary 00396 * 00397 150 CONTINUE 00398 IF( ISCALE.EQ.1 ) 00399 $ CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, 00400 $ D( LSV ), N, INFO ) 00401 IF( ISCALE.EQ.2 ) 00402 $ CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, 00403 $ D( LSV ), N, INFO ) 00404 * 00405 * Check for no convergence to an eigenvalue after a total 00406 * of N*MAXIT iterations. 00407 * 00408 IF( JTOT.LT.NMAXIT ) 00409 $ GO TO 10 00410 DO 160 I = 1, N - 1 00411 IF( E( I ).NE.ZERO ) 00412 $ INFO = INFO + 1 00413 160 CONTINUE 00414 GO TO 180 00415 * 00416 * Sort eigenvalues in increasing order. 00417 * 00418 170 CONTINUE 00419 CALL DLASRT( 'I', N, D, INFO ) 00420 * 00421 180 CONTINUE 00422 RETURN 00423 * 00424 * End of DSTERF 00425 * 00426 END