github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/internal/testdata/dsterftest/dsterf.f (about) 1 *> \brief \b DSTERF 2 * 3 * =========== DOCUMENTATION =========== 4 * 5 * Online html documentation available at 6 * http://www.netlib.org/lapack/explore-html/ 7 * 8 *> \htmlonly 9 *> Download DSTERF + dependencies 10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsterf.f"> 11 *> [TGZ]</a> 12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsterf.f"> 13 *> [ZIP]</a> 14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsterf.f"> 15 *> [TXT]</a> 16 *> \endhtmlonly 17 * 18 * Definition: 19 * =========== 20 * 21 * SUBROUTINE DSTERF( N, D, E, INFO ) 22 * 23 * .. Scalar Arguments .. 24 * INTEGER INFO, N 25 * .. 26 * .. Array Arguments .. 27 * DOUBLE PRECISION D( * ), E( * ) 28 * .. 29 * 30 * 31 *> \par Purpose: 32 * ============= 33 *> 34 *> \verbatim 35 *> 36 *> DSTERF computes all eigenvalues of a symmetric tridiagonal matrix 37 *> using the Pal-Walker-Kahan variant of the QL or QR algorithm. 38 *> \endverbatim 39 * 40 * Arguments: 41 * ========== 42 * 43 *> \param[in] N 44 *> \verbatim 45 *> N is INTEGER 46 *> The order of the matrix. N >= 0. 47 *> \endverbatim 48 *> 49 *> \param[in,out] D 50 *> \verbatim 51 *> D is DOUBLE PRECISION array, dimension (N) 52 *> On entry, the n diagonal elements of the tridiagonal matrix. 53 *> On exit, if INFO = 0, the eigenvalues in ascending order. 54 *> \endverbatim 55 *> 56 *> \param[in,out] E 57 *> \verbatim 58 *> E is DOUBLE PRECISION array, dimension (N-1) 59 *> On entry, the (n-1) subdiagonal elements of the tridiagonal 60 *> matrix. 61 *> On exit, E has been destroyed. 62 *> \endverbatim 63 *> 64 *> \param[out] INFO 65 *> \verbatim 66 *> INFO is INTEGER 67 *> = 0: successful exit 68 *> < 0: if INFO = -i, the i-th argument had an illegal value 69 *> > 0: the algorithm failed to find all of the eigenvalues in 70 *> a total of 30*N iterations; if INFO = i, then i 71 *> elements of E have not converged to zero. 72 *> \endverbatim 73 * 74 * Authors: 75 * ======== 76 * 77 *> \author Univ. of Tennessee 78 *> \author Univ. of California Berkeley 79 *> \author Univ. of Colorado Denver 80 *> \author NAG Ltd. 81 * 82 *> \date November 2011 83 * 84 *> \ingroup auxOTHERcomputational 85 * 86 * ===================================================================== 87 SUBROUTINE DSTERF( N, D, E, INFO ) 88 * 89 * -- LAPACK computational routine (version 3.4.0) -- 90 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 91 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 92 * November 2011 93 * 94 * .. Scalar Arguments .. 95 INTEGER INFO, N 96 * .. 97 * .. Array Arguments .. 98 DOUBLE PRECISION D( * ), E( * ) 99 * .. 100 * 101 * ===================================================================== 102 * 103 * .. Parameters .. 104 DOUBLE PRECISION ZERO, ONE, TWO, THREE 105 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 106 $ THREE = 3.0D0 ) 107 INTEGER MAXIT 108 PARAMETER ( MAXIT = 30 ) 109 * .. 110 * .. Local Scalars .. 111 INTEGER I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M, 112 $ NMAXIT 113 DOUBLE PRECISION ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC, 114 $ OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN, 115 $ SIGMA, SSFMAX, SSFMIN, RMAX 116 * .. 117 * .. External Functions .. 118 DOUBLE PRECISION DLAMCH, DLANST, DLAPY2 119 EXTERNAL DLAMCH, DLANST, DLAPY2 120 * .. 121 * .. External Subroutines .. 122 EXTERNAL DLAE2, DLASCL, DLASRT, XERBLA 123 * .. 124 * .. Intrinsic Functions .. 125 INTRINSIC ABS, SIGN, SQRT 126 * .. 127 * .. Executable Statements .. 128 * 129 * Test the input parameters. 130 * 131 INFO = 0 132 * 133 * Quick return if possible 134 * 135 IF( N.LT.0 ) THEN 136 INFO = -1 137 CALL XERBLA( 'DSTERF', -INFO ) 138 RETURN 139 END IF 140 IF( N.LE.1 ) 141 $ RETURN 142 * 143 * Determine the unit roundoff for this environment. 144 * 145 EPS = DLAMCH( 'E' ) 146 EPS2 = EPS**2 147 SAFMIN = DLAMCH( 'S' ) 148 SAFMAX = ONE / SAFMIN 149 SSFMAX = SQRT( SAFMAX ) / THREE 150 SSFMIN = SQRT( SAFMIN ) / EPS2 151 RMAX = DLAMCH( 'O' ) 152 * 153 * Compute the eigenvalues of the tridiagonal matrix. 154 * 155 NMAXIT = N*MAXIT 156 SIGMA = ZERO 157 JTOT = 0 158 * 159 * Determine where the matrix splits and choose QL or QR iteration 160 * for each block, according to whether top or bottom diagonal 161 * element is smaller. 162 * 163 L1 = 1 164 * 165 10 CONTINUE 166 print *, "l1 = ", l1 167 IF( L1.GT.N ) THEN 168 print *, "going to 170" 169 GO TO 170 170 end if 171 IF( L1.GT.1 ) 172 $ E( L1-1 ) = ZERO 173 DO 20 M = L1, N - 1 174 IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ 175 $ 1 ) ) ) )*EPS ) THEN 176 E( M ) = ZERO 177 GO TO 30 178 END IF 179 20 CONTINUE 180 M = N 181 * 182 30 CONTINUE 183 print *, "30, d" 184 print *, d(1:n) 185 L = L1 186 LSV = L 187 LEND = M 188 LENDSV = LEND 189 L1 = M + 1 190 IF( LEND.EQ.L ) 191 $ GO TO 10 192 * 193 * Scale submatrix in rows and columns L to LEND 194 * 195 ANORM = DLANST( 'M', LEND-L+1, D( L ), E( L ) ) 196 ISCALE = 0 197 IF( ANORM.EQ.ZERO ) 198 $ GO TO 10 199 IF( (ANORM.GT.SSFMAX) ) THEN 200 ISCALE = 1 201 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, 202 $ INFO ) 203 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, 204 $ INFO ) 205 ELSE IF( ANORM.LT.SSFMIN ) THEN 206 ISCALE = 2 207 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, 208 $ INFO ) 209 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, 210 $ INFO ) 211 END IF 212 * 213 DO 40 I = L, LEND - 1 214 E( I ) = E( I )**2 215 40 CONTINUE 216 * 217 * Choose between QL and QR iteration 218 * 219 IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN 220 LEND = LSV 221 L = LENDSV 222 END IF 223 * 224 IF( LEND.GE.L ) THEN 225 print *, "ql, d" 226 print *, d(1:n) 227 * 228 * QL Iteration 229 * 230 * Look for small subdiagonal element. 231 * 232 50 CONTINUE 233 IF( L.NE.LEND ) THEN 234 DO 60 M = L, LEND - 1 235 IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) ) 236 $ GO TO 70 237 60 CONTINUE 238 END IF 239 M = LEND 240 * 241 70 CONTINUE 242 IF( M.LT.LEND ) 243 $ E( M ) = ZERO 244 P = D( L ) 245 IF( M.EQ.L ) 246 $ GO TO 90 247 * 248 * If remaining matrix is 2 by 2, use DLAE2 to compute its 249 * eigenvalues. 250 * 251 IF( M.EQ.L+1 ) THEN 252 RTE = SQRT( E( L ) ) 253 CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 ) 254 D( L ) = RT1 255 D( L+1 ) = RT2 256 E( L ) = ZERO 257 L = L + 2 258 IF( L.LE.LEND ) 259 $ GO TO 50 260 GO TO 150 261 END IF 262 * 263 IF( JTOT.EQ.NMAXIT ) 264 $ GO TO 150 265 JTOT = JTOT + 1 266 * 267 * Form shift. 268 * 269 RTE = SQRT( E( L ) ) 270 SIGMA = ( D( L+1 )-P ) / ( TWO*RTE ) 271 R = DLAPY2( SIGMA, ONE ) 272 SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) 273 * 274 C = ONE 275 S = ZERO 276 GAMMA = D( M ) - SIGMA 277 P = GAMMA*GAMMA 278 * 279 * Inner loop 280 * 281 print *, "inner loop d before" 282 print *, d(1:n) 283 284 DO 80 I = M - 1, L, -1 285 print *, "inner loop" 286 print *, "ei", e(i) 287 BB = E( I ) 288 R = P + BB 289 print *, "bb,p,r" 290 print *, bb,p,r 291 IF( I.NE.M-1 ) THEN 292 print *, s,r 293 E( I+1 ) = S*R 294 end if 295 OLDC = C 296 C = P / R 297 S = BB / R 298 OLDGAM = GAMMA 299 print *, "di", d(i) 300 ALPHA = D( I ) 301 GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM 302 print *,"og, a, ga", OLDGAM, ALPHA, GAMMA 303 D( I+1 ) = OLDGAM + ( ALPHA-GAMMA ) 304 IF( C.NE.ZERO ) THEN 305 P = ( GAMMA*GAMMA ) / C 306 ELSE 307 P = OLDC*BB 308 END IF 309 print *, "p, gamma = ", p,GAMMA 310 80 CONTINUE 311 * 312 E( L ) = S*P 313 D( L ) = SIGMA + GAMMA 314 315 print *, "inner loop d after" 316 print *, d(1:n) 317 GO TO 50 318 * 319 * Eigenvalue found. 320 * 321 90 CONTINUE 322 D( L ) = P 323 * 324 L = L + 1 325 IF( L.LE.LEND ) 326 $ GO TO 50 327 GO TO 150 328 * 329 ELSE 330 * 331 * QR Iteration 332 * 333 * Look for small superdiagonal element. 334 * 335 100 CONTINUE 336 DO 110 M = L, LEND + 1, -1 337 IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) ) 338 $ GO TO 120 339 110 CONTINUE 340 M = LEND 341 * 342 120 CONTINUE 343 IF( M.GT.LEND ) 344 $ E( M-1 ) = ZERO 345 P = D( L ) 346 IF( M.EQ.L ) 347 $ GO TO 140 348 * 349 * If remaining matrix is 2 by 2, use DLAE2 to compute its 350 * eigenvalues. 351 * 352 IF( M.EQ.L-1 ) THEN 353 RTE = SQRT( E( L-1 ) ) 354 CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 ) 355 D( L ) = RT1 356 D( L-1 ) = RT2 357 E( L-1 ) = ZERO 358 L = L - 2 359 IF( L.GE.LEND ) 360 $ GO TO 100 361 GO TO 150 362 END IF 363 * 364 IF( JTOT.EQ.NMAXIT ) 365 $ GO TO 150 366 JTOT = JTOT + 1 367 * 368 * Form shift. 369 * 370 RTE = SQRT( E( L-1 ) ) 371 SIGMA = ( D( L-1 )-P ) / ( TWO*RTE ) 372 R = DLAPY2( SIGMA, ONE ) 373 SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) 374 * 375 C = ONE 376 S = ZERO 377 GAMMA = D( M ) - SIGMA 378 P = GAMMA*GAMMA 379 * 380 * Inner loop 381 * 382 DO 130 I = M, L - 1 383 BB = E( I ) 384 R = P + BB 385 IF( I.NE.M ) 386 $ E( I-1 ) = S*R 387 OLDC = C 388 C = P / R 389 S = BB / R 390 OLDGAM = GAMMA 391 ALPHA = D( I+1 ) 392 GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM 393 D( I ) = OLDGAM + ( ALPHA-GAMMA ) 394 IF( C.NE.ZERO ) THEN 395 P = ( GAMMA*GAMMA ) / C 396 ELSE 397 P = OLDC*BB 398 END IF 399 130 CONTINUE 400 * 401 E( L-1 ) = S*P 402 D( L ) = SIGMA + GAMMA 403 GO TO 100 404 * 405 * Eigenvalue found. 406 * 407 140 CONTINUE 408 D( L ) = P 409 * 410 L = L - 1 411 IF( L.GE.LEND ) 412 $ GO TO 100 413 GO TO 150 414 * 415 END IF 416 * 417 * Undo scaling if necessary 418 * 419 150 CONTINUE 420 IF( ISCALE.EQ.1 ) 421 $ CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, 422 $ D( LSV ), N, INFO ) 423 IF( ISCALE.EQ.2 ) 424 $ CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, 425 $ D( LSV ), N, INFO ) 426 * 427 * Check for no convergence to an eigenvalue after a total 428 * of N*MAXIT iterations. 429 * 430 IF( JTOT.LT.NMAXIT ) 431 $ GO TO 10 432 DO 160 I = 1, N - 1 433 IF( E( I ).NE.ZERO ) 434 $ INFO = INFO + 1 435 160 CONTINUE 436 GO TO 180 437 * 438 * Sort eigenvalues in increasing order. 439 * 440 170 CONTINUE 441 CALL DLASRT( 'I', N, D, INFO ) 442 * 443 180 CONTINUE 444 RETURN 445 * 446 * End of DSTERF 447 * 448 END