Annotation of rpl/lapack/lapack/dlasq4.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
! 2: $ DN1, DN2, TAU, TTYPE, G )
! 3: *
! 4: * -- LAPACK routine (version 3.2) --
! 5: *
! 6: * -- Contributed by Osni Marques of the Lawrence Berkeley National --
! 7: * -- Laboratory and Beresford Parlett of the Univ. of California at --
! 8: * -- Berkeley --
! 9: * -- November 2008 --
! 10: *
! 11: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 12: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 13: *
! 14: * .. Scalar Arguments ..
! 15: INTEGER I0, N0, N0IN, PP, TTYPE
! 16: DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
! 17: * ..
! 18: * .. Array Arguments ..
! 19: DOUBLE PRECISION Z( * )
! 20: * ..
! 21: *
! 22: * Purpose
! 23: * =======
! 24: *
! 25: * DLASQ4 computes an approximation TAU to the smallest eigenvalue
! 26: * using values of d from the previous transform.
! 27: *
! 28: * I0 (input) INTEGER
! 29: * First index.
! 30: *
! 31: * N0 (input) INTEGER
! 32: * Last index.
! 33: *
! 34: * Z (input) DOUBLE PRECISION array, dimension ( 4*N )
! 35: * Z holds the qd array.
! 36: *
! 37: * PP (input) INTEGER
! 38: * PP=0 for ping, PP=1 for pong.
! 39: *
! 40: * NOIN (input) INTEGER
! 41: * The value of N0 at start of EIGTEST.
! 42: *
! 43: * DMIN (input) DOUBLE PRECISION
! 44: * Minimum value of d.
! 45: *
! 46: * DMIN1 (input) DOUBLE PRECISION
! 47: * Minimum value of d, excluding D( N0 ).
! 48: *
! 49: * DMIN2 (input) DOUBLE PRECISION
! 50: * Minimum value of d, excluding D( N0 ) and D( N0-1 ).
! 51: *
! 52: * DN (input) DOUBLE PRECISION
! 53: * d(N)
! 54: *
! 55: * DN1 (input) DOUBLE PRECISION
! 56: * d(N-1)
! 57: *
! 58: * DN2 (input) DOUBLE PRECISION
! 59: * d(N-2)
! 60: *
! 61: * TAU (output) DOUBLE PRECISION
! 62: * This is the shift.
! 63: *
! 64: * TTYPE (output) INTEGER
! 65: * Shift type.
! 66: *
! 67: * G (input/output) REAL
! 68: * G is passed as an argument in order to save its value between
! 69: * calls to DLASQ4.
! 70: *
! 71: * Further Details
! 72: * ===============
! 73: * CNST1 = 9/16
! 74: *
! 75: * =====================================================================
! 76: *
! 77: * .. Parameters ..
! 78: DOUBLE PRECISION CNST1, CNST2, CNST3
! 79: PARAMETER ( CNST1 = 0.5630D0, CNST2 = 1.010D0,
! 80: $ CNST3 = 1.050D0 )
! 81: DOUBLE PRECISION QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD
! 82: PARAMETER ( QURTR = 0.250D0, THIRD = 0.3330D0,
! 83: $ HALF = 0.50D0, ZERO = 0.0D0, ONE = 1.0D0,
! 84: $ TWO = 2.0D0, HUNDRD = 100.0D0 )
! 85: * ..
! 86: * .. Local Scalars ..
! 87: INTEGER I4, NN, NP
! 88: DOUBLE PRECISION A2, B1, B2, GAM, GAP1, GAP2, S
! 89: * ..
! 90: * .. Intrinsic Functions ..
! 91: INTRINSIC MAX, MIN, SQRT
! 92: * ..
! 93: * .. Executable Statements ..
! 94: *
! 95: * A negative DMIN forces the shift to take that absolute value
! 96: * TTYPE records the type of shift.
! 97: *
! 98: IF( DMIN.LE.ZERO ) THEN
! 99: TAU = -DMIN
! 100: TTYPE = -1
! 101: RETURN
! 102: END IF
! 103: *
! 104: NN = 4*N0 + PP
! 105: IF( N0IN.EQ.N0 ) THEN
! 106: *
! 107: * No eigenvalues deflated.
! 108: *
! 109: IF( DMIN.EQ.DN .OR. DMIN.EQ.DN1 ) THEN
! 110: *
! 111: B1 = SQRT( Z( NN-3 ) )*SQRT( Z( NN-5 ) )
! 112: B2 = SQRT( Z( NN-7 ) )*SQRT( Z( NN-9 ) )
! 113: A2 = Z( NN-7 ) + Z( NN-5 )
! 114: *
! 115: * Cases 2 and 3.
! 116: *
! 117: IF( DMIN.EQ.DN .AND. DMIN1.EQ.DN1 ) THEN
! 118: GAP2 = DMIN2 - A2 - DMIN2*QURTR
! 119: IF( GAP2.GT.ZERO .AND. GAP2.GT.B2 ) THEN
! 120: GAP1 = A2 - DN - ( B2 / GAP2 )*B2
! 121: ELSE
! 122: GAP1 = A2 - DN - ( B1+B2 )
! 123: END IF
! 124: IF( GAP1.GT.ZERO .AND. GAP1.GT.B1 ) THEN
! 125: S = MAX( DN-( B1 / GAP1 )*B1, HALF*DMIN )
! 126: TTYPE = -2
! 127: ELSE
! 128: S = ZERO
! 129: IF( DN.GT.B1 )
! 130: $ S = DN - B1
! 131: IF( A2.GT.( B1+B2 ) )
! 132: $ S = MIN( S, A2-( B1+B2 ) )
! 133: S = MAX( S, THIRD*DMIN )
! 134: TTYPE = -3
! 135: END IF
! 136: ELSE
! 137: *
! 138: * Case 4.
! 139: *
! 140: TTYPE = -4
! 141: S = QURTR*DMIN
! 142: IF( DMIN.EQ.DN ) THEN
! 143: GAM = DN
! 144: A2 = ZERO
! 145: IF( Z( NN-5 ) .GT. Z( NN-7 ) )
! 146: $ RETURN
! 147: B2 = Z( NN-5 ) / Z( NN-7 )
! 148: NP = NN - 9
! 149: ELSE
! 150: NP = NN - 2*PP
! 151: B2 = Z( NP-2 )
! 152: GAM = DN1
! 153: IF( Z( NP-4 ) .GT. Z( NP-2 ) )
! 154: $ RETURN
! 155: A2 = Z( NP-4 ) / Z( NP-2 )
! 156: IF( Z( NN-9 ) .GT. Z( NN-11 ) )
! 157: $ RETURN
! 158: B2 = Z( NN-9 ) / Z( NN-11 )
! 159: NP = NN - 13
! 160: END IF
! 161: *
! 162: * Approximate contribution to norm squared from I < NN-1.
! 163: *
! 164: A2 = A2 + B2
! 165: DO 10 I4 = NP, 4*I0 - 1 + PP, -4
! 166: IF( B2.EQ.ZERO )
! 167: $ GO TO 20
! 168: B1 = B2
! 169: IF( Z( I4 ) .GT. Z( I4-2 ) )
! 170: $ RETURN
! 171: B2 = B2*( Z( I4 ) / Z( I4-2 ) )
! 172: A2 = A2 + B2
! 173: IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
! 174: $ GO TO 20
! 175: 10 CONTINUE
! 176: 20 CONTINUE
! 177: A2 = CNST3*A2
! 178: *
! 179: * Rayleigh quotient residual bound.
! 180: *
! 181: IF( A2.LT.CNST1 )
! 182: $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
! 183: END IF
! 184: ELSE IF( DMIN.EQ.DN2 ) THEN
! 185: *
! 186: * Case 5.
! 187: *
! 188: TTYPE = -5
! 189: S = QURTR*DMIN
! 190: *
! 191: * Compute contribution to norm squared from I > NN-2.
! 192: *
! 193: NP = NN - 2*PP
! 194: B1 = Z( NP-2 )
! 195: B2 = Z( NP-6 )
! 196: GAM = DN2
! 197: IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 )
! 198: $ RETURN
! 199: A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 )
! 200: *
! 201: * Approximate contribution to norm squared from I < NN-2.
! 202: *
! 203: IF( N0-I0.GT.2 ) THEN
! 204: B2 = Z( NN-13 ) / Z( NN-15 )
! 205: A2 = A2 + B2
! 206: DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4
! 207: IF( B2.EQ.ZERO )
! 208: $ GO TO 40
! 209: B1 = B2
! 210: IF( Z( I4 ) .GT. Z( I4-2 ) )
! 211: $ RETURN
! 212: B2 = B2*( Z( I4 ) / Z( I4-2 ) )
! 213: A2 = A2 + B2
! 214: IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
! 215: $ GO TO 40
! 216: 30 CONTINUE
! 217: 40 CONTINUE
! 218: A2 = CNST3*A2
! 219: END IF
! 220: *
! 221: IF( A2.LT.CNST1 )
! 222: $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
! 223: ELSE
! 224: *
! 225: * Case 6, no information to guide us.
! 226: *
! 227: IF( TTYPE.EQ.-6 ) THEN
! 228: G = G + THIRD*( ONE-G )
! 229: ELSE IF( TTYPE.EQ.-18 ) THEN
! 230: G = QURTR*THIRD
! 231: ELSE
! 232: G = QURTR
! 233: END IF
! 234: S = G*DMIN
! 235: TTYPE = -6
! 236: END IF
! 237: *
! 238: ELSE IF( N0IN.EQ.( N0+1 ) ) THEN
! 239: *
! 240: * One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
! 241: *
! 242: IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN
! 243: *
! 244: * Cases 7 and 8.
! 245: *
! 246: TTYPE = -7
! 247: S = THIRD*DMIN1
! 248: IF( Z( NN-5 ).GT.Z( NN-7 ) )
! 249: $ RETURN
! 250: B1 = Z( NN-5 ) / Z( NN-7 )
! 251: B2 = B1
! 252: IF( B2.EQ.ZERO )
! 253: $ GO TO 60
! 254: DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
! 255: A2 = B1
! 256: IF( Z( I4 ).GT.Z( I4-2 ) )
! 257: $ RETURN
! 258: B1 = B1*( Z( I4 ) / Z( I4-2 ) )
! 259: B2 = B2 + B1
! 260: IF( HUNDRD*MAX( B1, A2 ).LT.B2 )
! 261: $ GO TO 60
! 262: 50 CONTINUE
! 263: 60 CONTINUE
! 264: B2 = SQRT( CNST3*B2 )
! 265: A2 = DMIN1 / ( ONE+B2**2 )
! 266: GAP2 = HALF*DMIN2 - A2
! 267: IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
! 268: S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
! 269: ELSE
! 270: S = MAX( S, A2*( ONE-CNST2*B2 ) )
! 271: TTYPE = -8
! 272: END IF
! 273: ELSE
! 274: *
! 275: * Case 9.
! 276: *
! 277: S = QURTR*DMIN1
! 278: IF( DMIN1.EQ.DN1 )
! 279: $ S = HALF*DMIN1
! 280: TTYPE = -9
! 281: END IF
! 282: *
! 283: ELSE IF( N0IN.EQ.( N0+2 ) ) THEN
! 284: *
! 285: * Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
! 286: *
! 287: * Cases 10 and 11.
! 288: *
! 289: IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN
! 290: TTYPE = -10
! 291: S = THIRD*DMIN2
! 292: IF( Z( NN-5 ).GT.Z( NN-7 ) )
! 293: $ RETURN
! 294: B1 = Z( NN-5 ) / Z( NN-7 )
! 295: B2 = B1
! 296: IF( B2.EQ.ZERO )
! 297: $ GO TO 80
! 298: DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
! 299: IF( Z( I4 ).GT.Z( I4-2 ) )
! 300: $ RETURN
! 301: B1 = B1*( Z( I4 ) / Z( I4-2 ) )
! 302: B2 = B2 + B1
! 303: IF( HUNDRD*B1.LT.B2 )
! 304: $ GO TO 80
! 305: 70 CONTINUE
! 306: 80 CONTINUE
! 307: B2 = SQRT( CNST3*B2 )
! 308: A2 = DMIN2 / ( ONE+B2**2 )
! 309: GAP2 = Z( NN-7 ) + Z( NN-9 ) -
! 310: $ SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2
! 311: IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
! 312: S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
! 313: ELSE
! 314: S = MAX( S, A2*( ONE-CNST2*B2 ) )
! 315: END IF
! 316: ELSE
! 317: S = QURTR*DMIN2
! 318: TTYPE = -11
! 319: END IF
! 320: ELSE IF( N0IN.GT.( N0+2 ) ) THEN
! 321: *
! 322: * Case 12, more than two eigenvalues deflated. No information.
! 323: *
! 324: S = ZERO
! 325: TTYPE = -12
! 326: END IF
! 327: *
! 328: TAU = S
! 329: RETURN
! 330: *
! 331: * End of DLASQ4
! 332: *
! 333: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasq4.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack