Annotation of rpl/lapack/lapack/dsterf.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DSTERF( N, D, E, INFO )
! 2: *
! 3: * -- LAPACK routine (version 3.2) --
! 4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 6: * November 2006
! 7: *
! 8: * .. Scalar Arguments ..
! 9: INTEGER INFO, N
! 10: * ..
! 11: * .. Array Arguments ..
! 12: DOUBLE PRECISION D( * ), E( * )
! 13: * ..
! 14: *
! 15: * Purpose
! 16: * =======
! 17: *
! 18: * DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
! 19: * using the Pal-Walker-Kahan variant of the QL or QR algorithm.
! 20: *
! 21: * Arguments
! 22: * =========
! 23: *
! 24: * N (input) INTEGER
! 25: * The order of the matrix. N >= 0.
! 26: *
! 27: * D (input/output) DOUBLE PRECISION array, dimension (N)
! 28: * On entry, the n diagonal elements of the tridiagonal matrix.
! 29: * On exit, if INFO = 0, the eigenvalues in ascending order.
! 30: *
! 31: * E (input/output) DOUBLE PRECISION array, dimension (N-1)
! 32: * On entry, the (n-1) subdiagonal elements of the tridiagonal
! 33: * matrix.
! 34: * On exit, E has been destroyed.
! 35: *
! 36: * INFO (output) INTEGER
! 37: * = 0: successful exit
! 38: * < 0: if INFO = -i, the i-th argument had an illegal value
! 39: * > 0: the algorithm failed to find all of the eigenvalues in
! 40: * a total of 30*N iterations; if INFO = i, then i
! 41: * elements of E have not converged to zero.
! 42: *
! 43: * =====================================================================
! 44: *
! 45: * .. Parameters ..
! 46: DOUBLE PRECISION ZERO, ONE, TWO, THREE
! 47: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
! 48: $ THREE = 3.0D0 )
! 49: INTEGER MAXIT
! 50: PARAMETER ( MAXIT = 30 )
! 51: * ..
! 52: * .. Local Scalars ..
! 53: INTEGER I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M,
! 54: $ NMAXIT
! 55: DOUBLE PRECISION ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
! 56: $ OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
! 57: $ SIGMA, SSFMAX, SSFMIN
! 58: * ..
! 59: * .. External Functions ..
! 60: DOUBLE PRECISION DLAMCH, DLANST, DLAPY2
! 61: EXTERNAL DLAMCH, DLANST, DLAPY2
! 62: * ..
! 63: * .. External Subroutines ..
! 64: EXTERNAL DLAE2, DLASCL, DLASRT, XERBLA
! 65: * ..
! 66: * .. Intrinsic Functions ..
! 67: INTRINSIC ABS, SIGN, SQRT
! 68: * ..
! 69: * .. Executable Statements ..
! 70: *
! 71: * Test the input parameters.
! 72: *
! 73: INFO = 0
! 74: *
! 75: * Quick return if possible
! 76: *
! 77: IF( N.LT.0 ) THEN
! 78: INFO = -1
! 79: CALL XERBLA( 'DSTERF', -INFO )
! 80: RETURN
! 81: END IF
! 82: IF( N.LE.1 )
! 83: $ RETURN
! 84: *
! 85: * Determine the unit roundoff for this environment.
! 86: *
! 87: EPS = DLAMCH( 'E' )
! 88: EPS2 = EPS**2
! 89: SAFMIN = DLAMCH( 'S' )
! 90: SAFMAX = ONE / SAFMIN
! 91: SSFMAX = SQRT( SAFMAX ) / THREE
! 92: SSFMIN = SQRT( SAFMIN ) / EPS2
! 93: *
! 94: * Compute the eigenvalues of the tridiagonal matrix.
! 95: *
! 96: NMAXIT = N*MAXIT
! 97: SIGMA = ZERO
! 98: JTOT = 0
! 99: *
! 100: * Determine where the matrix splits and choose QL or QR iteration
! 101: * for each block, according to whether top or bottom diagonal
! 102: * element is smaller.
! 103: *
! 104: L1 = 1
! 105: *
! 106: 10 CONTINUE
! 107: IF( L1.GT.N )
! 108: $ GO TO 170
! 109: IF( L1.GT.1 )
! 110: $ E( L1-1 ) = ZERO
! 111: DO 20 M = L1, N - 1
! 112: IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
! 113: $ 1 ) ) ) )*EPS ) THEN
! 114: E( M ) = ZERO
! 115: GO TO 30
! 116: END IF
! 117: 20 CONTINUE
! 118: M = N
! 119: *
! 120: 30 CONTINUE
! 121: L = L1
! 122: LSV = L
! 123: LEND = M
! 124: LENDSV = LEND
! 125: L1 = M + 1
! 126: IF( LEND.EQ.L )
! 127: $ GO TO 10
! 128: *
! 129: * Scale submatrix in rows and columns L to LEND
! 130: *
! 131: ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) )
! 132: ISCALE = 0
! 133: IF( ANORM.GT.SSFMAX ) THEN
! 134: ISCALE = 1
! 135: CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
! 136: $ INFO )
! 137: CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
! 138: $ INFO )
! 139: ELSE IF( ANORM.LT.SSFMIN ) THEN
! 140: ISCALE = 2
! 141: CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
! 142: $ INFO )
! 143: CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
! 144: $ INFO )
! 145: END IF
! 146: *
! 147: DO 40 I = L, LEND - 1
! 148: E( I ) = E( I )**2
! 149: 40 CONTINUE
! 150: *
! 151: * Choose between QL and QR iteration
! 152: *
! 153: IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
! 154: LEND = LSV
! 155: L = LENDSV
! 156: END IF
! 157: *
! 158: IF( LEND.GE.L ) THEN
! 159: *
! 160: * QL Iteration
! 161: *
! 162: * Look for small subdiagonal element.
! 163: *
! 164: 50 CONTINUE
! 165: IF( L.NE.LEND ) THEN
! 166: DO 60 M = L, LEND - 1
! 167: IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) )
! 168: $ GO TO 70
! 169: 60 CONTINUE
! 170: END IF
! 171: M = LEND
! 172: *
! 173: 70 CONTINUE
! 174: IF( M.LT.LEND )
! 175: $ E( M ) = ZERO
! 176: P = D( L )
! 177: IF( M.EQ.L )
! 178: $ GO TO 90
! 179: *
! 180: * If remaining matrix is 2 by 2, use DLAE2 to compute its
! 181: * eigenvalues.
! 182: *
! 183: IF( M.EQ.L+1 ) THEN
! 184: RTE = SQRT( E( L ) )
! 185: CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 )
! 186: D( L ) = RT1
! 187: D( L+1 ) = RT2
! 188: E( L ) = ZERO
! 189: L = L + 2
! 190: IF( L.LE.LEND )
! 191: $ GO TO 50
! 192: GO TO 150
! 193: END IF
! 194: *
! 195: IF( JTOT.EQ.NMAXIT )
! 196: $ GO TO 150
! 197: JTOT = JTOT + 1
! 198: *
! 199: * Form shift.
! 200: *
! 201: RTE = SQRT( E( L ) )
! 202: SIGMA = ( D( L+1 )-P ) / ( TWO*RTE )
! 203: R = DLAPY2( SIGMA, ONE )
! 204: SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
! 205: *
! 206: C = ONE
! 207: S = ZERO
! 208: GAMMA = D( M ) - SIGMA
! 209: P = GAMMA*GAMMA
! 210: *
! 211: * Inner loop
! 212: *
! 213: DO 80 I = M - 1, L, -1
! 214: BB = E( I )
! 215: R = P + BB
! 216: IF( I.NE.M-1 )
! 217: $ E( I+1 ) = S*R
! 218: OLDC = C
! 219: C = P / R
! 220: S = BB / R
! 221: OLDGAM = GAMMA
! 222: ALPHA = D( I )
! 223: GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
! 224: D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )
! 225: IF( C.NE.ZERO ) THEN
! 226: P = ( GAMMA*GAMMA ) / C
! 227: ELSE
! 228: P = OLDC*BB
! 229: END IF
! 230: 80 CONTINUE
! 231: *
! 232: E( L ) = S*P
! 233: D( L ) = SIGMA + GAMMA
! 234: GO TO 50
! 235: *
! 236: * Eigenvalue found.
! 237: *
! 238: 90 CONTINUE
! 239: D( L ) = P
! 240: *
! 241: L = L + 1
! 242: IF( L.LE.LEND )
! 243: $ GO TO 50
! 244: GO TO 150
! 245: *
! 246: ELSE
! 247: *
! 248: * QR Iteration
! 249: *
! 250: * Look for small superdiagonal element.
! 251: *
! 252: 100 CONTINUE
! 253: DO 110 M = L, LEND + 1, -1
! 254: IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) )
! 255: $ GO TO 120
! 256: 110 CONTINUE
! 257: M = LEND
! 258: *
! 259: 120 CONTINUE
! 260: IF( M.GT.LEND )
! 261: $ E( M-1 ) = ZERO
! 262: P = D( L )
! 263: IF( M.EQ.L )
! 264: $ GO TO 140
! 265: *
! 266: * If remaining matrix is 2 by 2, use DLAE2 to compute its
! 267: * eigenvalues.
! 268: *
! 269: IF( M.EQ.L-1 ) THEN
! 270: RTE = SQRT( E( L-1 ) )
! 271: CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 )
! 272: D( L ) = RT1
! 273: D( L-1 ) = RT2
! 274: E( L-1 ) = ZERO
! 275: L = L - 2
! 276: IF( L.GE.LEND )
! 277: $ GO TO 100
! 278: GO TO 150
! 279: END IF
! 280: *
! 281: IF( JTOT.EQ.NMAXIT )
! 282: $ GO TO 150
! 283: JTOT = JTOT + 1
! 284: *
! 285: * Form shift.
! 286: *
! 287: RTE = SQRT( E( L-1 ) )
! 288: SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )
! 289: R = DLAPY2( SIGMA, ONE )
! 290: SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
! 291: *
! 292: C = ONE
! 293: S = ZERO
! 294: GAMMA = D( M ) - SIGMA
! 295: P = GAMMA*GAMMA
! 296: *
! 297: * Inner loop
! 298: *
! 299: DO 130 I = M, L - 1
! 300: BB = E( I )
! 301: R = P + BB
! 302: IF( I.NE.M )
! 303: $ E( I-1 ) = S*R
! 304: OLDC = C
! 305: C = P / R
! 306: S = BB / R
! 307: OLDGAM = GAMMA
! 308: ALPHA = D( I+1 )
! 309: GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
! 310: D( I ) = OLDGAM + ( ALPHA-GAMMA )
! 311: IF( C.NE.ZERO ) THEN
! 312: P = ( GAMMA*GAMMA ) / C
! 313: ELSE
! 314: P = OLDC*BB
! 315: END IF
! 316: 130 CONTINUE
! 317: *
! 318: E( L-1 ) = S*P
! 319: D( L ) = SIGMA + GAMMA
! 320: GO TO 100
! 321: *
! 322: * Eigenvalue found.
! 323: *
! 324: 140 CONTINUE
! 325: D( L ) = P
! 326: *
! 327: L = L - 1
! 328: IF( L.GE.LEND )
! 329: $ GO TO 100
! 330: GO TO 150
! 331: *
! 332: END IF
! 333: *
! 334: * Undo scaling if necessary
! 335: *
! 336: 150 CONTINUE
! 337: IF( ISCALE.EQ.1 )
! 338: $ CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
! 339: $ D( LSV ), N, INFO )
! 340: IF( ISCALE.EQ.2 )
! 341: $ CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
! 342: $ D( LSV ), N, INFO )
! 343: *
! 344: * Check for no convergence to an eigenvalue after a total
! 345: * of N*MAXIT iterations.
! 346: *
! 347: IF( JTOT.LT.NMAXIT )
! 348: $ GO TO 10
! 349: DO 160 I = 1, N - 1
! 350: IF( E( I ).NE.ZERO )
! 351: $ INFO = INFO + 1
! 352: 160 CONTINUE
! 353: GO TO 180
! 354: *
! 355: * Sort eigenvalues in increasing order.
! 356: *
! 357: 170 CONTINUE
! 358: CALL DLASRT( 'I', N, D, INFO )
! 359: *
! 360: 180 CONTINUE
! 361: RETURN
! 362: *
! 363: * End of DSTERF
! 364: *
! 365: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsterf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack