Annotation of rpl/lapack/lapack/dlasq2.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLASQ2( N, Z, INFO )
! 2: *
! 3: * -- LAPACK routine (version 3.2) --
! 4: *
! 5: * -- Contributed by Osni Marques of the Lawrence Berkeley National --
! 6: * -- Laboratory and Beresford Parlett of the Univ. of California at --
! 7: * -- Berkeley --
! 8: * -- November 2008 --
! 9: *
! 10: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 11: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 12: *
! 13: * .. Scalar Arguments ..
! 14: INTEGER INFO, N
! 15: * ..
! 16: * .. Array Arguments ..
! 17: DOUBLE PRECISION Z( * )
! 18: * ..
! 19: *
! 20: * Purpose
! 21: * =======
! 22: *
! 23: * DLASQ2 computes all the eigenvalues of the symmetric positive
! 24: * definite tridiagonal matrix associated with the qd array Z to high
! 25: * relative accuracy are computed to high relative accuracy, in the
! 26: * absence of denormalization, underflow and overflow.
! 27: *
! 28: * To see the relation of Z to the tridiagonal matrix, let L be a
! 29: * unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
! 30: * let U be an upper bidiagonal matrix with 1's above and diagonal
! 31: * Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
! 32: * symmetric tridiagonal to which it is similar.
! 33: *
! 34: * Note : DLASQ2 defines a logical variable, IEEE, which is true
! 35: * on machines which follow ieee-754 floating-point standard in their
! 36: * handling of infinities and NaNs, and false otherwise. This variable
! 37: * is passed to DLASQ3.
! 38: *
! 39: * Arguments
! 40: * =========
! 41: *
! 42: * N (input) INTEGER
! 43: * The number of rows and columns in the matrix. N >= 0.
! 44: *
! 45: * Z (input/output) DOUBLE PRECISION array, dimension ( 4*N )
! 46: * On entry Z holds the qd array. On exit, entries 1 to N hold
! 47: * the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
! 48: * trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
! 49: * N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
! 50: * holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
! 51: * shifts that failed.
! 52: *
! 53: * INFO (output) INTEGER
! 54: * = 0: successful exit
! 55: * < 0: if the i-th argument is a scalar and had an illegal
! 56: * value, then INFO = -i, if the i-th argument is an
! 57: * array and the j-entry had an illegal value, then
! 58: * INFO = -(i*100+j)
! 59: * > 0: the algorithm failed
! 60: * = 1, a split was marked by a positive value in E
! 61: * = 2, current block of Z not diagonalized after 30*N
! 62: * iterations (in inner while loop)
! 63: * = 3, termination criterion of outer while loop not met
! 64: * (program created more than N unreduced blocks)
! 65: *
! 66: * Further Details
! 67: * ===============
! 68: * Local Variables: I0:N0 defines a current unreduced segment of Z.
! 69: * The shifts are accumulated in SIGMA. Iteration count is in ITER.
! 70: * Ping-pong is controlled by PP (alternates between 0 and 1).
! 71: *
! 72: * =====================================================================
! 73: *
! 74: * .. Parameters ..
! 75: DOUBLE PRECISION CBIAS
! 76: PARAMETER ( CBIAS = 1.50D0 )
! 77: DOUBLE PRECISION ZERO, HALF, ONE, TWO, FOUR, HUNDRD
! 78: PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0,
! 79: $ TWO = 2.0D0, FOUR = 4.0D0, HUNDRD = 100.0D0 )
! 80: * ..
! 81: * .. Local Scalars ..
! 82: LOGICAL IEEE
! 83: INTEGER I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K,
! 84: $ KMIN, N0, NBIG, NDIV, NFAIL, PP, SPLT, TTYPE
! 85: DOUBLE PRECISION D, DEE, DEEMIN, DESIG, DMIN, DMIN1, DMIN2, DN,
! 86: $ DN1, DN2, E, EMAX, EMIN, EPS, G, OLDEMN, QMAX,
! 87: $ QMIN, S, SAFMIN, SIGMA, T, TAU, TEMP, TOL,
! 88: $ TOL2, TRACE, ZMAX
! 89: * ..
! 90: * .. External Subroutines ..
! 91: EXTERNAL DLASQ3, DLASRT, XERBLA
! 92: * ..
! 93: * .. External Functions ..
! 94: INTEGER ILAENV
! 95: DOUBLE PRECISION DLAMCH
! 96: EXTERNAL DLAMCH, ILAENV
! 97: * ..
! 98: * .. Intrinsic Functions ..
! 99: INTRINSIC ABS, DBLE, MAX, MIN, SQRT
! 100: * ..
! 101: * .. Executable Statements ..
! 102: *
! 103: * Test the input arguments.
! 104: * (in case DLASQ2 is not called by DLASQ1)
! 105: *
! 106: INFO = 0
! 107: EPS = DLAMCH( 'Precision' )
! 108: SAFMIN = DLAMCH( 'Safe minimum' )
! 109: TOL = EPS*HUNDRD
! 110: TOL2 = TOL**2
! 111: *
! 112: IF( N.LT.0 ) THEN
! 113: INFO = -1
! 114: CALL XERBLA( 'DLASQ2', 1 )
! 115: RETURN
! 116: ELSE IF( N.EQ.0 ) THEN
! 117: RETURN
! 118: ELSE IF( N.EQ.1 ) THEN
! 119: *
! 120: * 1-by-1 case.
! 121: *
! 122: IF( Z( 1 ).LT.ZERO ) THEN
! 123: INFO = -201
! 124: CALL XERBLA( 'DLASQ2', 2 )
! 125: END IF
! 126: RETURN
! 127: ELSE IF( N.EQ.2 ) THEN
! 128: *
! 129: * 2-by-2 case.
! 130: *
! 131: IF( Z( 2 ).LT.ZERO .OR. Z( 3 ).LT.ZERO ) THEN
! 132: INFO = -2
! 133: CALL XERBLA( 'DLASQ2', 2 )
! 134: RETURN
! 135: ELSE IF( Z( 3 ).GT.Z( 1 ) ) THEN
! 136: D = Z( 3 )
! 137: Z( 3 ) = Z( 1 )
! 138: Z( 1 ) = D
! 139: END IF
! 140: Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 )
! 141: IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN
! 142: T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) )
! 143: S = Z( 3 )*( Z( 2 ) / T )
! 144: IF( S.LE.T ) THEN
! 145: S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) )
! 146: ELSE
! 147: S = Z( 3 )*( Z( 2 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
! 148: END IF
! 149: T = Z( 1 ) + ( S+Z( 2 ) )
! 150: Z( 3 ) = Z( 3 )*( Z( 1 ) / T )
! 151: Z( 1 ) = T
! 152: END IF
! 153: Z( 2 ) = Z( 3 )
! 154: Z( 6 ) = Z( 2 ) + Z( 1 )
! 155: RETURN
! 156: END IF
! 157: *
! 158: * Check for negative data and compute sums of q's and e's.
! 159: *
! 160: Z( 2*N ) = ZERO
! 161: EMIN = Z( 2 )
! 162: QMAX = ZERO
! 163: ZMAX = ZERO
! 164: D = ZERO
! 165: E = ZERO
! 166: *
! 167: DO 10 K = 1, 2*( N-1 ), 2
! 168: IF( Z( K ).LT.ZERO ) THEN
! 169: INFO = -( 200+K )
! 170: CALL XERBLA( 'DLASQ2', 2 )
! 171: RETURN
! 172: ELSE IF( Z( K+1 ).LT.ZERO ) THEN
! 173: INFO = -( 200+K+1 )
! 174: CALL XERBLA( 'DLASQ2', 2 )
! 175: RETURN
! 176: END IF
! 177: D = D + Z( K )
! 178: E = E + Z( K+1 )
! 179: QMAX = MAX( QMAX, Z( K ) )
! 180: EMIN = MIN( EMIN, Z( K+1 ) )
! 181: ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) )
! 182: 10 CONTINUE
! 183: IF( Z( 2*N-1 ).LT.ZERO ) THEN
! 184: INFO = -( 200+2*N-1 )
! 185: CALL XERBLA( 'DLASQ2', 2 )
! 186: RETURN
! 187: END IF
! 188: D = D + Z( 2*N-1 )
! 189: QMAX = MAX( QMAX, Z( 2*N-1 ) )
! 190: ZMAX = MAX( QMAX, ZMAX )
! 191: *
! 192: * Check for diagonality.
! 193: *
! 194: IF( E.EQ.ZERO ) THEN
! 195: DO 20 K = 2, N
! 196: Z( K ) = Z( 2*K-1 )
! 197: 20 CONTINUE
! 198: CALL DLASRT( 'D', N, Z, IINFO )
! 199: Z( 2*N-1 ) = D
! 200: RETURN
! 201: END IF
! 202: *
! 203: TRACE = D + E
! 204: *
! 205: * Check for zero data.
! 206: *
! 207: IF( TRACE.EQ.ZERO ) THEN
! 208: Z( 2*N-1 ) = ZERO
! 209: RETURN
! 210: END IF
! 211: *
! 212: * Check whether the machine is IEEE conformable.
! 213: *
! 214: IEEE = ILAENV( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND.
! 215: $ ILAENV( 11, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1
! 216: *
! 217: * Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
! 218: *
! 219: DO 30 K = 2*N, 2, -2
! 220: Z( 2*K ) = ZERO
! 221: Z( 2*K-1 ) = Z( K )
! 222: Z( 2*K-2 ) = ZERO
! 223: Z( 2*K-3 ) = Z( K-1 )
! 224: 30 CONTINUE
! 225: *
! 226: I0 = 1
! 227: N0 = N
! 228: *
! 229: * Reverse the qd-array, if warranted.
! 230: *
! 231: IF( CBIAS*Z( 4*I0-3 ).LT.Z( 4*N0-3 ) ) THEN
! 232: IPN4 = 4*( I0+N0 )
! 233: DO 40 I4 = 4*I0, 2*( I0+N0-1 ), 4
! 234: TEMP = Z( I4-3 )
! 235: Z( I4-3 ) = Z( IPN4-I4-3 )
! 236: Z( IPN4-I4-3 ) = TEMP
! 237: TEMP = Z( I4-1 )
! 238: Z( I4-1 ) = Z( IPN4-I4-5 )
! 239: Z( IPN4-I4-5 ) = TEMP
! 240: 40 CONTINUE
! 241: END IF
! 242: *
! 243: * Initial split checking via dqd and Li's test.
! 244: *
! 245: PP = 0
! 246: *
! 247: DO 80 K = 1, 2
! 248: *
! 249: D = Z( 4*N0+PP-3 )
! 250: DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4
! 251: IF( Z( I4-1 ).LE.TOL2*D ) THEN
! 252: Z( I4-1 ) = -ZERO
! 253: D = Z( I4-3 )
! 254: ELSE
! 255: D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) )
! 256: END IF
! 257: 50 CONTINUE
! 258: *
! 259: * dqd maps Z to ZZ plus Li's test.
! 260: *
! 261: EMIN = Z( 4*I0+PP+1 )
! 262: D = Z( 4*I0+PP-3 )
! 263: DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4
! 264: Z( I4-2*PP-2 ) = D + Z( I4-1 )
! 265: IF( Z( I4-1 ).LE.TOL2*D ) THEN
! 266: Z( I4-1 ) = -ZERO
! 267: Z( I4-2*PP-2 ) = D
! 268: Z( I4-2*PP ) = ZERO
! 269: D = Z( I4+1 )
! 270: ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND.
! 271: $ SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN
! 272: TEMP = Z( I4+1 ) / Z( I4-2*PP-2 )
! 273: Z( I4-2*PP ) = Z( I4-1 )*TEMP
! 274: D = D*TEMP
! 275: ELSE
! 276: Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) )
! 277: D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) )
! 278: END IF
! 279: EMIN = MIN( EMIN, Z( I4-2*PP ) )
! 280: 60 CONTINUE
! 281: Z( 4*N0-PP-2 ) = D
! 282: *
! 283: * Now find qmax.
! 284: *
! 285: QMAX = Z( 4*I0-PP-2 )
! 286: DO 70 I4 = 4*I0 - PP + 2, 4*N0 - PP - 2, 4
! 287: QMAX = MAX( QMAX, Z( I4 ) )
! 288: 70 CONTINUE
! 289: *
! 290: * Prepare for the next iteration on K.
! 291: *
! 292: PP = 1 - PP
! 293: 80 CONTINUE
! 294: *
! 295: * Initialise variables to pass to DLASQ3.
! 296: *
! 297: TTYPE = 0
! 298: DMIN1 = ZERO
! 299: DMIN2 = ZERO
! 300: DN = ZERO
! 301: DN1 = ZERO
! 302: DN2 = ZERO
! 303: G = ZERO
! 304: TAU = ZERO
! 305: *
! 306: ITER = 2
! 307: NFAIL = 0
! 308: NDIV = 2*( N0-I0 )
! 309: *
! 310: DO 160 IWHILA = 1, N + 1
! 311: IF( N0.LT.1 )
! 312: $ GO TO 170
! 313: *
! 314: * While array unfinished do
! 315: *
! 316: * E(N0) holds the value of SIGMA when submatrix in I0:N0
! 317: * splits from the rest of the array, but is negated.
! 318: *
! 319: DESIG = ZERO
! 320: IF( N0.EQ.N ) THEN
! 321: SIGMA = ZERO
! 322: ELSE
! 323: SIGMA = -Z( 4*N0-1 )
! 324: END IF
! 325: IF( SIGMA.LT.ZERO ) THEN
! 326: INFO = 1
! 327: RETURN
! 328: END IF
! 329: *
! 330: * Find last unreduced submatrix's top index I0, find QMAX and
! 331: * EMIN. Find Gershgorin-type bound if Q's much greater than E's.
! 332: *
! 333: EMAX = ZERO
! 334: IF( N0.GT.I0 ) THEN
! 335: EMIN = ABS( Z( 4*N0-5 ) )
! 336: ELSE
! 337: EMIN = ZERO
! 338: END IF
! 339: QMIN = Z( 4*N0-3 )
! 340: QMAX = QMIN
! 341: DO 90 I4 = 4*N0, 8, -4
! 342: IF( Z( I4-5 ).LE.ZERO )
! 343: $ GO TO 100
! 344: IF( QMIN.GE.FOUR*EMAX ) THEN
! 345: QMIN = MIN( QMIN, Z( I4-3 ) )
! 346: EMAX = MAX( EMAX, Z( I4-5 ) )
! 347: END IF
! 348: QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) )
! 349: EMIN = MIN( EMIN, Z( I4-5 ) )
! 350: 90 CONTINUE
! 351: I4 = 4
! 352: *
! 353: 100 CONTINUE
! 354: I0 = I4 / 4
! 355: PP = 0
! 356: *
! 357: IF( N0-I0.GT.1 ) THEN
! 358: DEE = Z( 4*I0-3 )
! 359: DEEMIN = DEE
! 360: KMIN = I0
! 361: DO 110 I4 = 4*I0+1, 4*N0-3, 4
! 362: DEE = Z( I4 )*( DEE /( DEE+Z( I4-2 ) ) )
! 363: IF( DEE.LE.DEEMIN ) THEN
! 364: DEEMIN = DEE
! 365: KMIN = ( I4+3 )/4
! 366: END IF
! 367: 110 CONTINUE
! 368: IF( (KMIN-I0)*2.LT.N0-KMIN .AND.
! 369: $ DEEMIN.LE.HALF*Z(4*N0-3) ) THEN
! 370: IPN4 = 4*( I0+N0 )
! 371: PP = 2
! 372: DO 120 I4 = 4*I0, 2*( I0+N0-1 ), 4
! 373: TEMP = Z( I4-3 )
! 374: Z( I4-3 ) = Z( IPN4-I4-3 )
! 375: Z( IPN4-I4-3 ) = TEMP
! 376: TEMP = Z( I4-2 )
! 377: Z( I4-2 ) = Z( IPN4-I4-2 )
! 378: Z( IPN4-I4-2 ) = TEMP
! 379: TEMP = Z( I4-1 )
! 380: Z( I4-1 ) = Z( IPN4-I4-5 )
! 381: Z( IPN4-I4-5 ) = TEMP
! 382: TEMP = Z( I4 )
! 383: Z( I4 ) = Z( IPN4-I4-4 )
! 384: Z( IPN4-I4-4 ) = TEMP
! 385: 120 CONTINUE
! 386: END IF
! 387: END IF
! 388: *
! 389: * Put -(initial shift) into DMIN.
! 390: *
! 391: DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) )
! 392: *
! 393: * Now I0:N0 is unreduced.
! 394: * PP = 0 for ping, PP = 1 for pong.
! 395: * PP = 2 indicates that flipping was applied to the Z array and
! 396: * and that the tests for deflation upon entry in DLASQ3
! 397: * should not be performed.
! 398: *
! 399: NBIG = 30*( N0-I0+1 )
! 400: DO 140 IWHILB = 1, NBIG
! 401: IF( I0.GT.N0 )
! 402: $ GO TO 150
! 403: *
! 404: * While submatrix unfinished take a good dqds step.
! 405: *
! 406: CALL DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
! 407: $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
! 408: $ DN2, G, TAU )
! 409: *
! 410: PP = 1 - PP
! 411: *
! 412: * When EMIN is very small check for splits.
! 413: *
! 414: IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN
! 415: IF( Z( 4*N0 ).LE.TOL2*QMAX .OR.
! 416: $ Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN
! 417: SPLT = I0 - 1
! 418: QMAX = Z( 4*I0-3 )
! 419: EMIN = Z( 4*I0-1 )
! 420: OLDEMN = Z( 4*I0 )
! 421: DO 130 I4 = 4*I0, 4*( N0-3 ), 4
! 422: IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR.
! 423: $ Z( I4-1 ).LE.TOL2*SIGMA ) THEN
! 424: Z( I4-1 ) = -SIGMA
! 425: SPLT = I4 / 4
! 426: QMAX = ZERO
! 427: EMIN = Z( I4+3 )
! 428: OLDEMN = Z( I4+4 )
! 429: ELSE
! 430: QMAX = MAX( QMAX, Z( I4+1 ) )
! 431: EMIN = MIN( EMIN, Z( I4-1 ) )
! 432: OLDEMN = MIN( OLDEMN, Z( I4 ) )
! 433: END IF
! 434: 130 CONTINUE
! 435: Z( 4*N0-1 ) = EMIN
! 436: Z( 4*N0 ) = OLDEMN
! 437: I0 = SPLT + 1
! 438: END IF
! 439: END IF
! 440: *
! 441: 140 CONTINUE
! 442: *
! 443: INFO = 2
! 444: RETURN
! 445: *
! 446: * end IWHILB
! 447: *
! 448: 150 CONTINUE
! 449: *
! 450: 160 CONTINUE
! 451: *
! 452: INFO = 3
! 453: RETURN
! 454: *
! 455: * end IWHILA
! 456: *
! 457: 170 CONTINUE
! 458: *
! 459: * Move q's to the front.
! 460: *
! 461: DO 180 K = 2, N
! 462: Z( K ) = Z( 4*K-3 )
! 463: 180 CONTINUE
! 464: *
! 465: * Sort and compute sum of eigenvalues.
! 466: *
! 467: CALL DLASRT( 'D', N, Z, IINFO )
! 468: *
! 469: E = ZERO
! 470: DO 190 K = N, 1, -1
! 471: E = E + Z( K )
! 472: 190 CONTINUE
! 473: *
! 474: * Store trace, sum(eigenvalues) and information on performance.
! 475: *
! 476: Z( 2*N+1 ) = TRACE
! 477: Z( 2*N+2 ) = E
! 478: Z( 2*N+3 ) = DBLE( ITER )
! 479: Z( 2*N+4 ) = DBLE( NDIV ) / DBLE( N**2 )
! 480: Z( 2*N+5 ) = HUNDRD*NFAIL / DBLE( ITER )
! 481: RETURN
! 482: *
! 483: * End of DLASQ2
! 484: *
! 485: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasq2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack