![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
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