Annotation of rpl/lapack/lapack/dlasq3.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
! 2: $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
! 3: $ DN2, G, TAU )
! 4: *
! 5: * -- LAPACK routine (version 3.2) --
! 6: *
! 7: * -- Contributed by Osni Marques of the Lawrence Berkeley National --
! 8: * -- Laboratory and Beresford Parlett of the Univ. of California at --
! 9: * -- Berkeley --
! 10: * -- November 2008 --
! 11: *
! 12: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 13: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 14: *
! 15: * .. Scalar Arguments ..
! 16: LOGICAL IEEE
! 17: INTEGER I0, ITER, N0, NDIV, NFAIL, PP
! 18: DOUBLE PRECISION DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
! 19: $ QMAX, SIGMA, TAU
! 20: * ..
! 21: * .. Array Arguments ..
! 22: DOUBLE PRECISION Z( * )
! 23: * ..
! 24: *
! 25: * Purpose
! 26: * =======
! 27: *
! 28: * DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
! 29: * In case of failure it changes shifts, and tries again until output
! 30: * is positive.
! 31: *
! 32: * Arguments
! 33: * =========
! 34: *
! 35: * I0 (input) INTEGER
! 36: * First index.
! 37: *
! 38: * N0 (input) INTEGER
! 39: * Last index.
! 40: *
! 41: * Z (input) DOUBLE PRECISION array, dimension ( 4*N )
! 42: * Z holds the qd array.
! 43: *
! 44: * PP (input/output) INTEGER
! 45: * PP=0 for ping, PP=1 for pong.
! 46: * PP=2 indicates that flipping was applied to the Z array
! 47: * and that the initial tests for deflation should not be
! 48: * performed.
! 49: *
! 50: * DMIN (output) DOUBLE PRECISION
! 51: * Minimum value of d.
! 52: *
! 53: * SIGMA (output) DOUBLE PRECISION
! 54: * Sum of shifts used in current segment.
! 55: *
! 56: * DESIG (input/output) DOUBLE PRECISION
! 57: * Lower order part of SIGMA
! 58: *
! 59: * QMAX (input) DOUBLE PRECISION
! 60: * Maximum value of q.
! 61: *
! 62: * NFAIL (output) INTEGER
! 63: * Number of times shift was too big.
! 64: *
! 65: * ITER (output) INTEGER
! 66: * Number of iterations.
! 67: *
! 68: * NDIV (output) INTEGER
! 69: * Number of divisions.
! 70: *
! 71: * IEEE (input) LOGICAL
! 72: * Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
! 73: *
! 74: * TTYPE (input/output) INTEGER
! 75: * Shift type.
! 76: *
! 77: * DMIN1, DMIN2, DN, DN1, DN2, G, TAU (input/output) DOUBLE PRECISION
! 78: * These are passed as arguments in order to save their values
! 79: * between calls to DLASQ3.
! 80: *
! 81: * =====================================================================
! 82: *
! 83: * .. Parameters ..
! 84: DOUBLE PRECISION CBIAS
! 85: PARAMETER ( CBIAS = 1.50D0 )
! 86: DOUBLE PRECISION ZERO, QURTR, HALF, ONE, TWO, HUNDRD
! 87: PARAMETER ( ZERO = 0.0D0, QURTR = 0.250D0, HALF = 0.5D0,
! 88: $ ONE = 1.0D0, TWO = 2.0D0, HUNDRD = 100.0D0 )
! 89: * ..
! 90: * .. Local Scalars ..
! 91: INTEGER IPN4, J4, N0IN, NN, TTYPE
! 92: DOUBLE PRECISION EPS, S, T, TEMP, TOL, TOL2
! 93: * ..
! 94: * .. External Subroutines ..
! 95: EXTERNAL DLASQ4, DLASQ5, DLASQ6
! 96: * ..
! 97: * .. External Function ..
! 98: DOUBLE PRECISION DLAMCH
! 99: LOGICAL DISNAN
! 100: EXTERNAL DISNAN, DLAMCH
! 101: * ..
! 102: * .. Intrinsic Functions ..
! 103: INTRINSIC ABS, MAX, MIN, SQRT
! 104: * ..
! 105: * .. Executable Statements ..
! 106: *
! 107: N0IN = N0
! 108: EPS = DLAMCH( 'Precision' )
! 109: TOL = EPS*HUNDRD
! 110: TOL2 = TOL**2
! 111: *
! 112: * Check for deflation.
! 113: *
! 114: 10 CONTINUE
! 115: *
! 116: IF( N0.LT.I0 )
! 117: $ RETURN
! 118: IF( N0.EQ.I0 )
! 119: $ GO TO 20
! 120: NN = 4*N0 + PP
! 121: IF( N0.EQ.( I0+1 ) )
! 122: $ GO TO 40
! 123: *
! 124: * Check whether E(N0-1) is negligible, 1 eigenvalue.
! 125: *
! 126: IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
! 127: $ Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
! 128: $ GO TO 30
! 129: *
! 130: 20 CONTINUE
! 131: *
! 132: Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
! 133: N0 = N0 - 1
! 134: GO TO 10
! 135: *
! 136: * Check whether E(N0-2) is negligible, 2 eigenvalues.
! 137: *
! 138: 30 CONTINUE
! 139: *
! 140: IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
! 141: $ Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
! 142: $ GO TO 50
! 143: *
! 144: 40 CONTINUE
! 145: *
! 146: IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
! 147: S = Z( NN-3 )
! 148: Z( NN-3 ) = Z( NN-7 )
! 149: Z( NN-7 ) = S
! 150: END IF
! 151: IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2 ) THEN
! 152: T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
! 153: S = Z( NN-3 )*( Z( NN-5 ) / T )
! 154: IF( S.LE.T ) THEN
! 155: S = Z( NN-3 )*( Z( NN-5 ) /
! 156: $ ( T*( ONE+SQRT( ONE+S / T ) ) ) )
! 157: ELSE
! 158: S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
! 159: END IF
! 160: T = Z( NN-7 ) + ( S+Z( NN-5 ) )
! 161: Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
! 162: Z( NN-7 ) = T
! 163: END IF
! 164: Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
! 165: Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
! 166: N0 = N0 - 2
! 167: GO TO 10
! 168: *
! 169: 50 CONTINUE
! 170: IF( PP.EQ.2 )
! 171: $ PP = 0
! 172: *
! 173: * Reverse the qd-array, if warranted.
! 174: *
! 175: IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
! 176: IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
! 177: IPN4 = 4*( I0+N0 )
! 178: DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
! 179: TEMP = Z( J4-3 )
! 180: Z( J4-3 ) = Z( IPN4-J4-3 )
! 181: Z( IPN4-J4-3 ) = TEMP
! 182: TEMP = Z( J4-2 )
! 183: Z( J4-2 ) = Z( IPN4-J4-2 )
! 184: Z( IPN4-J4-2 ) = TEMP
! 185: TEMP = Z( J4-1 )
! 186: Z( J4-1 ) = Z( IPN4-J4-5 )
! 187: Z( IPN4-J4-5 ) = TEMP
! 188: TEMP = Z( J4 )
! 189: Z( J4 ) = Z( IPN4-J4-4 )
! 190: Z( IPN4-J4-4 ) = TEMP
! 191: 60 CONTINUE
! 192: IF( N0-I0.LE.4 ) THEN
! 193: Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
! 194: Z( 4*N0-PP ) = Z( 4*I0-PP )
! 195: END IF
! 196: DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
! 197: Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
! 198: $ Z( 4*I0+PP+3 ) )
! 199: Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
! 200: $ Z( 4*I0-PP+4 ) )
! 201: QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
! 202: DMIN = -ZERO
! 203: END IF
! 204: END IF
! 205: *
! 206: * Choose a shift.
! 207: *
! 208: CALL DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
! 209: $ DN2, TAU, TTYPE, G )
! 210: *
! 211: * Call dqds until DMIN > 0.
! 212: *
! 213: 70 CONTINUE
! 214: *
! 215: CALL DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
! 216: $ DN1, DN2, IEEE )
! 217: *
! 218: NDIV = NDIV + ( N0-I0+2 )
! 219: ITER = ITER + 1
! 220: *
! 221: * Check status.
! 222: *
! 223: IF( DMIN.GE.ZERO .AND. DMIN1.GT.ZERO ) THEN
! 224: *
! 225: * Success.
! 226: *
! 227: GO TO 90
! 228: *
! 229: ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND.
! 230: $ Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
! 231: $ ABS( DN ).LT.TOL*SIGMA ) THEN
! 232: *
! 233: * Convergence hidden by negative DN.
! 234: *
! 235: Z( 4*( N0-1 )-PP+2 ) = ZERO
! 236: DMIN = ZERO
! 237: GO TO 90
! 238: ELSE IF( DMIN.LT.ZERO ) THEN
! 239: *
! 240: * TAU too big. Select new TAU and try again.
! 241: *
! 242: NFAIL = NFAIL + 1
! 243: IF( TTYPE.LT.-22 ) THEN
! 244: *
! 245: * Failed twice. Play it safe.
! 246: *
! 247: TAU = ZERO
! 248: ELSE IF( DMIN1.GT.ZERO ) THEN
! 249: *
! 250: * Late failure. Gives excellent shift.
! 251: *
! 252: TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
! 253: TTYPE = TTYPE - 11
! 254: ELSE
! 255: *
! 256: * Early failure. Divide by 4.
! 257: *
! 258: TAU = QURTR*TAU
! 259: TTYPE = TTYPE - 12
! 260: END IF
! 261: GO TO 70
! 262: ELSE IF( DISNAN( DMIN ) ) THEN
! 263: *
! 264: * NaN.
! 265: *
! 266: IF( TAU.EQ.ZERO ) THEN
! 267: GO TO 80
! 268: ELSE
! 269: TAU = ZERO
! 270: GO TO 70
! 271: END IF
! 272: ELSE
! 273: *
! 274: * Possible underflow. Play it safe.
! 275: *
! 276: GO TO 80
! 277: END IF
! 278: *
! 279: * Risk of underflow.
! 280: *
! 281: 80 CONTINUE
! 282: CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
! 283: NDIV = NDIV + ( N0-I0+2 )
! 284: ITER = ITER + 1
! 285: TAU = ZERO
! 286: *
! 287: 90 CONTINUE
! 288: IF( TAU.LT.SIGMA ) THEN
! 289: DESIG = DESIG + TAU
! 290: T = SIGMA + DESIG
! 291: DESIG = DESIG - ( T-SIGMA )
! 292: ELSE
! 293: T = SIGMA + TAU
! 294: DESIG = SIGMA - ( T-TAU ) + DESIG
! 295: END IF
! 296: SIGMA = T
! 297: *
! 298: RETURN
! 299: *
! 300: * End of DLASQ3
! 301: *
! 302: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasq3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack