Annotation of rpl/lapack/lapack/dlaic1.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
! 2: *
! 3: * -- LAPACK auxiliary 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: * November 2006
! 7: *
! 8: * .. Scalar Arguments ..
! 9: INTEGER J, JOB
! 10: DOUBLE PRECISION C, GAMMA, S, SEST, SESTPR
! 11: * ..
! 12: * .. Array Arguments ..
! 13: DOUBLE PRECISION W( J ), X( J )
! 14: * ..
! 15: *
! 16: * Purpose
! 17: * =======
! 18: *
! 19: * DLAIC1 applies one step of incremental condition estimation in
! 20: * its simplest version:
! 21: *
! 22: * Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
! 23: * lower triangular matrix L, such that
! 24: * twonorm(L*x) = sest
! 25: * Then DLAIC1 computes sestpr, s, c such that
! 26: * the vector
! 27: * [ s*x ]
! 28: * xhat = [ c ]
! 29: * is an approximate singular vector of
! 30: * [ L 0 ]
! 31: * Lhat = [ w' gamma ]
! 32: * in the sense that
! 33: * twonorm(Lhat*xhat) = sestpr.
! 34: *
! 35: * Depending on JOB, an estimate for the largest or smallest singular
! 36: * value is computed.
! 37: *
! 38: * Note that [s c]' and sestpr**2 is an eigenpair of the system
! 39: *
! 40: * diag(sest*sest, 0) + [alpha gamma] * [ alpha ]
! 41: * [ gamma ]
! 42: *
! 43: * where alpha = x'*w.
! 44: *
! 45: * Arguments
! 46: * =========
! 47: *
! 48: * JOB (input) INTEGER
! 49: * = 1: an estimate for the largest singular value is computed.
! 50: * = 2: an estimate for the smallest singular value is computed.
! 51: *
! 52: * J (input) INTEGER
! 53: * Length of X and W
! 54: *
! 55: * X (input) DOUBLE PRECISION array, dimension (J)
! 56: * The j-vector x.
! 57: *
! 58: * SEST (input) DOUBLE PRECISION
! 59: * Estimated singular value of j by j matrix L
! 60: *
! 61: * W (input) DOUBLE PRECISION array, dimension (J)
! 62: * The j-vector w.
! 63: *
! 64: * GAMMA (input) DOUBLE PRECISION
! 65: * The diagonal element gamma.
! 66: *
! 67: * SESTPR (output) DOUBLE PRECISION
! 68: * Estimated singular value of (j+1) by (j+1) matrix Lhat.
! 69: *
! 70: * S (output) DOUBLE PRECISION
! 71: * Sine needed in forming xhat.
! 72: *
! 73: * C (output) DOUBLE PRECISION
! 74: * Cosine needed in forming xhat.
! 75: *
! 76: * =====================================================================
! 77: *
! 78: * .. Parameters ..
! 79: DOUBLE PRECISION ZERO, ONE, TWO
! 80: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
! 81: DOUBLE PRECISION HALF, FOUR
! 82: PARAMETER ( HALF = 0.5D0, FOUR = 4.0D0 )
! 83: * ..
! 84: * .. Local Scalars ..
! 85: DOUBLE PRECISION ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS,
! 86: $ NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2
! 87: * ..
! 88: * .. Intrinsic Functions ..
! 89: INTRINSIC ABS, MAX, SIGN, SQRT
! 90: * ..
! 91: * .. External Functions ..
! 92: DOUBLE PRECISION DDOT, DLAMCH
! 93: EXTERNAL DDOT, DLAMCH
! 94: * ..
! 95: * .. Executable Statements ..
! 96: *
! 97: EPS = DLAMCH( 'Epsilon' )
! 98: ALPHA = DDOT( J, X, 1, W, 1 )
! 99: *
! 100: ABSALP = ABS( ALPHA )
! 101: ABSGAM = ABS( GAMMA )
! 102: ABSEST = ABS( SEST )
! 103: *
! 104: IF( JOB.EQ.1 ) THEN
! 105: *
! 106: * Estimating largest singular value
! 107: *
! 108: * special cases
! 109: *
! 110: IF( SEST.EQ.ZERO ) THEN
! 111: S1 = MAX( ABSGAM, ABSALP )
! 112: IF( S1.EQ.ZERO ) THEN
! 113: S = ZERO
! 114: C = ONE
! 115: SESTPR = ZERO
! 116: ELSE
! 117: S = ALPHA / S1
! 118: C = GAMMA / S1
! 119: TMP = SQRT( S*S+C*C )
! 120: S = S / TMP
! 121: C = C / TMP
! 122: SESTPR = S1*TMP
! 123: END IF
! 124: RETURN
! 125: ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
! 126: S = ONE
! 127: C = ZERO
! 128: TMP = MAX( ABSEST, ABSALP )
! 129: S1 = ABSEST / TMP
! 130: S2 = ABSALP / TMP
! 131: SESTPR = TMP*SQRT( S1*S1+S2*S2 )
! 132: RETURN
! 133: ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
! 134: S1 = ABSGAM
! 135: S2 = ABSEST
! 136: IF( S1.LE.S2 ) THEN
! 137: S = ONE
! 138: C = ZERO
! 139: SESTPR = S2
! 140: ELSE
! 141: S = ZERO
! 142: C = ONE
! 143: SESTPR = S1
! 144: END IF
! 145: RETURN
! 146: ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
! 147: S1 = ABSGAM
! 148: S2 = ABSALP
! 149: IF( S1.LE.S2 ) THEN
! 150: TMP = S1 / S2
! 151: S = SQRT( ONE+TMP*TMP )
! 152: SESTPR = S2*S
! 153: C = ( GAMMA / S2 ) / S
! 154: S = SIGN( ONE, ALPHA ) / S
! 155: ELSE
! 156: TMP = S2 / S1
! 157: C = SQRT( ONE+TMP*TMP )
! 158: SESTPR = S1*C
! 159: S = ( ALPHA / S1 ) / C
! 160: C = SIGN( ONE, GAMMA ) / C
! 161: END IF
! 162: RETURN
! 163: ELSE
! 164: *
! 165: * normal case
! 166: *
! 167: ZETA1 = ALPHA / ABSEST
! 168: ZETA2 = GAMMA / ABSEST
! 169: *
! 170: B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF
! 171: C = ZETA1*ZETA1
! 172: IF( B.GT.ZERO ) THEN
! 173: T = C / ( B+SQRT( B*B+C ) )
! 174: ELSE
! 175: T = SQRT( B*B+C ) - B
! 176: END IF
! 177: *
! 178: SINE = -ZETA1 / T
! 179: COSINE = -ZETA2 / ( ONE+T )
! 180: TMP = SQRT( SINE*SINE+COSINE*COSINE )
! 181: S = SINE / TMP
! 182: C = COSINE / TMP
! 183: SESTPR = SQRT( T+ONE )*ABSEST
! 184: RETURN
! 185: END IF
! 186: *
! 187: ELSE IF( JOB.EQ.2 ) THEN
! 188: *
! 189: * Estimating smallest singular value
! 190: *
! 191: * special cases
! 192: *
! 193: IF( SEST.EQ.ZERO ) THEN
! 194: SESTPR = ZERO
! 195: IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN
! 196: SINE = ONE
! 197: COSINE = ZERO
! 198: ELSE
! 199: SINE = -GAMMA
! 200: COSINE = ALPHA
! 201: END IF
! 202: S1 = MAX( ABS( SINE ), ABS( COSINE ) )
! 203: S = SINE / S1
! 204: C = COSINE / S1
! 205: TMP = SQRT( S*S+C*C )
! 206: S = S / TMP
! 207: C = C / TMP
! 208: RETURN
! 209: ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
! 210: S = ZERO
! 211: C = ONE
! 212: SESTPR = ABSGAM
! 213: RETURN
! 214: ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
! 215: S1 = ABSGAM
! 216: S2 = ABSEST
! 217: IF( S1.LE.S2 ) THEN
! 218: S = ZERO
! 219: C = ONE
! 220: SESTPR = S1
! 221: ELSE
! 222: S = ONE
! 223: C = ZERO
! 224: SESTPR = S2
! 225: END IF
! 226: RETURN
! 227: ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
! 228: S1 = ABSGAM
! 229: S2 = ABSALP
! 230: IF( S1.LE.S2 ) THEN
! 231: TMP = S1 / S2
! 232: C = SQRT( ONE+TMP*TMP )
! 233: SESTPR = ABSEST*( TMP / C )
! 234: S = -( GAMMA / S2 ) / C
! 235: C = SIGN( ONE, ALPHA ) / C
! 236: ELSE
! 237: TMP = S2 / S1
! 238: S = SQRT( ONE+TMP*TMP )
! 239: SESTPR = ABSEST / S
! 240: C = ( ALPHA / S1 ) / S
! 241: S = -SIGN( ONE, GAMMA ) / S
! 242: END IF
! 243: RETURN
! 244: ELSE
! 245: *
! 246: * normal case
! 247: *
! 248: ZETA1 = ALPHA / ABSEST
! 249: ZETA2 = GAMMA / ABSEST
! 250: *
! 251: NORMA = MAX( ONE+ZETA1*ZETA1+ABS( ZETA1*ZETA2 ),
! 252: $ ABS( ZETA1*ZETA2 )+ZETA2*ZETA2 )
! 253: *
! 254: * See if root is closer to zero or to ONE
! 255: *
! 256: TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 )
! 257: IF( TEST.GE.ZERO ) THEN
! 258: *
! 259: * root is close to zero, compute directly
! 260: *
! 261: B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF
! 262: C = ZETA2*ZETA2
! 263: T = C / ( B+SQRT( ABS( B*B-C ) ) )
! 264: SINE = ZETA1 / ( ONE-T )
! 265: COSINE = -ZETA2 / T
! 266: SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST
! 267: ELSE
! 268: *
! 269: * root is closer to ONE, shift by that amount
! 270: *
! 271: B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF
! 272: C = ZETA1*ZETA1
! 273: IF( B.GE.ZERO ) THEN
! 274: T = -C / ( B+SQRT( B*B+C ) )
! 275: ELSE
! 276: T = B - SQRT( B*B+C )
! 277: END IF
! 278: SINE = -ZETA1 / T
! 279: COSINE = -ZETA2 / ( ONE+T )
! 280: SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST
! 281: END IF
! 282: TMP = SQRT( SINE*SINE+COSINE*COSINE )
! 283: S = SINE / TMP
! 284: C = COSINE / TMP
! 285: RETURN
! 286: *
! 287: END IF
! 288: END IF
! 289: RETURN
! 290: *
! 291: * End of DLAIC1
! 292: *
! 293: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaic1.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack