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