![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, 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: * February 2007 7: * 8: * .. Scalar Arguments .. 9: LOGICAL ORGATI 10: INTEGER INFO, KNITER 11: DOUBLE PRECISION FINIT, RHO, TAU 12: * .. 13: * .. Array Arguments .. 14: DOUBLE PRECISION D( 3 ), Z( 3 ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * DLAED6 computes the positive or negative root (closest to the origin) 21: * of 22: * z(1) z(2) z(3) 23: * f(x) = rho + --------- + ---------- + --------- 24: * d(1)-x d(2)-x d(3)-x 25: * 26: * It is assumed that 27: * 28: * if ORGATI = .true. the root is between d(2) and d(3); 29: * otherwise it is between d(1) and d(2) 30: * 31: * This routine will be called by DLAED4 when necessary. In most cases, 32: * the root sought is the smallest in magnitude, though it might not be 33: * in some extremely rare situations. 34: * 35: * Arguments 36: * ========= 37: * 38: * KNITER (input) INTEGER 39: * Refer to DLAED4 for its significance. 40: * 41: * ORGATI (input) LOGICAL 42: * If ORGATI is true, the needed root is between d(2) and 43: * d(3); otherwise it is between d(1) and d(2). See 44: * DLAED4 for further details. 45: * 46: * RHO (input) DOUBLE PRECISION 47: * Refer to the equation f(x) above. 48: * 49: * D (input) DOUBLE PRECISION array, dimension (3) 50: * D satisfies d(1) < d(2) < d(3). 51: * 52: * Z (input) DOUBLE PRECISION array, dimension (3) 53: * Each of the elements in z must be positive. 54: * 55: * FINIT (input) DOUBLE PRECISION 56: * The value of f at 0. It is more accurate than the one 57: * evaluated inside this routine (if someone wants to do 58: * so). 59: * 60: * TAU (output) DOUBLE PRECISION 61: * The root of the equation f(x). 62: * 63: * INFO (output) INTEGER 64: * = 0: successful exit 65: * > 0: if INFO = 1, failure to converge 66: * 67: * Further Details 68: * =============== 69: * 70: * 30/06/99: Based on contributions by 71: * Ren-Cang Li, Computer Science Division, University of California 72: * at Berkeley, USA 73: * 74: * 10/02/03: This version has a few statements commented out for thread 75: * safety (machine parameters are computed on each entry). SJH. 76: * 77: * 05/10/06: Modified from a new version of Ren-Cang Li, use 78: * Gragg-Thornton-Warner cubic convergent scheme for better stability. 79: * 80: * ===================================================================== 81: * 82: * .. Parameters .. 83: INTEGER MAXIT 84: PARAMETER ( MAXIT = 40 ) 85: DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT 86: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 87: $ THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0 ) 88: * .. 89: * .. External Functions .. 90: DOUBLE PRECISION DLAMCH 91: EXTERNAL DLAMCH 92: * .. 93: * .. Local Arrays .. 94: DOUBLE PRECISION DSCALE( 3 ), ZSCALE( 3 ) 95: * .. 96: * .. Local Scalars .. 97: LOGICAL SCALE 98: INTEGER I, ITER, NITER 99: DOUBLE PRECISION A, B, BASE, C, DDF, DF, EPS, ERRETM, ETA, F, 100: $ FC, SCLFAC, SCLINV, SMALL1, SMALL2, SMINV1, 101: $ SMINV2, TEMP, TEMP1, TEMP2, TEMP3, TEMP4, 102: $ LBD, UBD 103: * .. 104: * .. Intrinsic Functions .. 105: INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT 106: * .. 107: * .. Executable Statements .. 108: * 109: INFO = 0 110: * 111: IF( ORGATI ) THEN 112: LBD = D(2) 113: UBD = D(3) 114: ELSE 115: LBD = D(1) 116: UBD = D(2) 117: END IF 118: IF( FINIT .LT. ZERO )THEN 119: LBD = ZERO 120: ELSE 121: UBD = ZERO 122: END IF 123: * 124: NITER = 1 125: TAU = ZERO 126: IF( KNITER.EQ.2 ) THEN 127: IF( ORGATI ) THEN 128: TEMP = ( D( 3 )-D( 2 ) ) / TWO 129: C = RHO + Z( 1 ) / ( ( D( 1 )-D( 2 ) )-TEMP ) 130: A = C*( D( 2 )+D( 3 ) ) + Z( 2 ) + Z( 3 ) 131: B = C*D( 2 )*D( 3 ) + Z( 2 )*D( 3 ) + Z( 3 )*D( 2 ) 132: ELSE 133: TEMP = ( D( 1 )-D( 2 ) ) / TWO 134: C = RHO + Z( 3 ) / ( ( D( 3 )-D( 2 ) )-TEMP ) 135: A = C*( D( 1 )+D( 2 ) ) + Z( 1 ) + Z( 2 ) 136: B = C*D( 1 )*D( 2 ) + Z( 1 )*D( 2 ) + Z( 2 )*D( 1 ) 137: END IF 138: TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) ) 139: A = A / TEMP 140: B = B / TEMP 141: C = C / TEMP 142: IF( C.EQ.ZERO ) THEN 143: TAU = B / A 144: ELSE IF( A.LE.ZERO ) THEN 145: TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 146: ELSE 147: TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) 148: END IF 149: IF( TAU .LT. LBD .OR. TAU .GT. UBD ) 150: $ TAU = ( LBD+UBD )/TWO 151: IF( D(1).EQ.TAU .OR. D(2).EQ.TAU .OR. D(3).EQ.TAU ) THEN 152: TAU = ZERO 153: ELSE 154: TEMP = FINIT + TAU*Z(1)/( D(1)*( D( 1 )-TAU ) ) + 155: $ TAU*Z(2)/( D(2)*( D( 2 )-TAU ) ) + 156: $ TAU*Z(3)/( D(3)*( D( 3 )-TAU ) ) 157: IF( TEMP .LE. ZERO )THEN 158: LBD = TAU 159: ELSE 160: UBD = TAU 161: END IF 162: IF( ABS( FINIT ).LE.ABS( TEMP ) ) 163: $ TAU = ZERO 164: END IF 165: END IF 166: * 167: * get machine parameters for possible scaling to avoid overflow 168: * 169: * modified by Sven: parameters SMALL1, SMINV1, SMALL2, 170: * SMINV2, EPS are not SAVEd anymore between one call to the 171: * others but recomputed at each call 172: * 173: EPS = DLAMCH( 'Epsilon' ) 174: BASE = DLAMCH( 'Base' ) 175: SMALL1 = BASE**( INT( LOG( DLAMCH( 'SafMin' ) ) / LOG( BASE ) / 176: $ THREE ) ) 177: SMINV1 = ONE / SMALL1 178: SMALL2 = SMALL1*SMALL1 179: SMINV2 = SMINV1*SMINV1 180: * 181: * Determine if scaling of inputs necessary to avoid overflow 182: * when computing 1/TEMP**3 183: * 184: IF( ORGATI ) THEN 185: TEMP = MIN( ABS( D( 2 )-TAU ), ABS( D( 3 )-TAU ) ) 186: ELSE 187: TEMP = MIN( ABS( D( 1 )-TAU ), ABS( D( 2 )-TAU ) ) 188: END IF 189: SCALE = .FALSE. 190: IF( TEMP.LE.SMALL1 ) THEN 191: SCALE = .TRUE. 192: IF( TEMP.LE.SMALL2 ) THEN 193: * 194: * Scale up by power of radix nearest 1/SAFMIN**(2/3) 195: * 196: SCLFAC = SMINV2 197: SCLINV = SMALL2 198: ELSE 199: * 200: * Scale up by power of radix nearest 1/SAFMIN**(1/3) 201: * 202: SCLFAC = SMINV1 203: SCLINV = SMALL1 204: END IF 205: * 206: * Scaling up safe because D, Z, TAU scaled elsewhere to be O(1) 207: * 208: DO 10 I = 1, 3 209: DSCALE( I ) = D( I )*SCLFAC 210: ZSCALE( I ) = Z( I )*SCLFAC 211: 10 CONTINUE 212: TAU = TAU*SCLFAC 213: LBD = LBD*SCLFAC 214: UBD = UBD*SCLFAC 215: ELSE 216: * 217: * Copy D and Z to DSCALE and ZSCALE 218: * 219: DO 20 I = 1, 3 220: DSCALE( I ) = D( I ) 221: ZSCALE( I ) = Z( I ) 222: 20 CONTINUE 223: END IF 224: * 225: FC = ZERO 226: DF = ZERO 227: DDF = ZERO 228: DO 30 I = 1, 3 229: TEMP = ONE / ( DSCALE( I )-TAU ) 230: TEMP1 = ZSCALE( I )*TEMP 231: TEMP2 = TEMP1*TEMP 232: TEMP3 = TEMP2*TEMP 233: FC = FC + TEMP1 / DSCALE( I ) 234: DF = DF + TEMP2 235: DDF = DDF + TEMP3 236: 30 CONTINUE 237: F = FINIT + TAU*FC 238: * 239: IF( ABS( F ).LE.ZERO ) 240: $ GO TO 60 241: IF( F .LE. ZERO )THEN 242: LBD = TAU 243: ELSE 244: UBD = TAU 245: END IF 246: * 247: * Iteration begins -- Use Gragg-Thornton-Warner cubic convergent 248: * scheme 249: * 250: * It is not hard to see that 251: * 252: * 1) Iterations will go up monotonically 253: * if FINIT < 0; 254: * 255: * 2) Iterations will go down monotonically 256: * if FINIT > 0. 257: * 258: ITER = NITER + 1 259: * 260: DO 50 NITER = ITER, MAXIT 261: * 262: IF( ORGATI ) THEN 263: TEMP1 = DSCALE( 2 ) - TAU 264: TEMP2 = DSCALE( 3 ) - TAU 265: ELSE 266: TEMP1 = DSCALE( 1 ) - TAU 267: TEMP2 = DSCALE( 2 ) - TAU 268: END IF 269: A = ( TEMP1+TEMP2 )*F - TEMP1*TEMP2*DF 270: B = TEMP1*TEMP2*F 271: C = F - ( TEMP1+TEMP2 )*DF + TEMP1*TEMP2*DDF 272: TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) ) 273: A = A / TEMP 274: B = B / TEMP 275: C = C / TEMP 276: IF( C.EQ.ZERO ) THEN 277: ETA = B / A 278: ELSE IF( A.LE.ZERO ) THEN 279: ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 280: ELSE 281: ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) 282: END IF 283: IF( F*ETA.GE.ZERO ) THEN 284: ETA = -F / DF 285: END IF 286: * 287: TAU = TAU + ETA 288: IF( TAU .LT. LBD .OR. TAU .GT. UBD ) 289: $ TAU = ( LBD + UBD )/TWO 290: * 291: FC = ZERO 292: ERRETM = ZERO 293: DF = ZERO 294: DDF = ZERO 295: DO 40 I = 1, 3 296: TEMP = ONE / ( DSCALE( I )-TAU ) 297: TEMP1 = ZSCALE( I )*TEMP 298: TEMP2 = TEMP1*TEMP 299: TEMP3 = TEMP2*TEMP 300: TEMP4 = TEMP1 / DSCALE( I ) 301: FC = FC + TEMP4 302: ERRETM = ERRETM + ABS( TEMP4 ) 303: DF = DF + TEMP2 304: DDF = DDF + TEMP3 305: 40 CONTINUE 306: F = FINIT + TAU*FC 307: ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) + 308: $ ABS( TAU )*DF 309: IF( ABS( F ).LE.EPS*ERRETM ) 310: $ GO TO 60 311: IF( F .LE. ZERO )THEN 312: LBD = TAU 313: ELSE 314: UBD = TAU 315: END IF 316: 50 CONTINUE 317: INFO = 1 318: 60 CONTINUE 319: * 320: * Undo scaling 321: * 322: IF( SCALE ) 323: $ TAU = TAU*SCLINV 324: RETURN 325: * 326: * End of DLAED6 327: * 328: END