Annotation of rpl/lapack/lapack/dlaed6.f, revision 1.1
1.1 ! bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaed6.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack