![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
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