001: SUBROUTINE SLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C ) 002: * 003: * -- LAPACK auxiliary routine (version 3.2) -- 004: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 005: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 006: * November 2006 007: * 008: * .. Scalar Arguments .. 009: INTEGER J, JOB 010: REAL C, GAMMA, S, SEST, SESTPR 011: * .. 012: * .. Array Arguments .. 013: REAL W( J ), X( J ) 014: * .. 015: * 016: * Purpose 017: * ======= 018: * 019: * SLAIC1 applies one step of incremental condition estimation in 020: * its simplest version: 021: * 022: * Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j 023: * lower triangular matrix L, such that 024: * twonorm(L*x) = sest 025: * Then SLAIC1 computes sestpr, s, c such that 026: * the vector 027: * [ s*x ] 028: * xhat = [ c ] 029: * is an approximate singular vector of 030: * [ L 0 ] 031: * Lhat = [ w' gamma ] 032: * in the sense that 033: * twonorm(Lhat*xhat) = sestpr. 034: * 035: * Depending on JOB, an estimate for the largest or smallest singular 036: * value is computed. 037: * 038: * Note that [s c]' and sestpr**2 is an eigenpair of the system 039: * 040: * diag(sest*sest, 0) + [alpha gamma] * [ alpha ] 041: * [ gamma ] 042: * 043: * where alpha = x'*w. 044: * 045: * Arguments 046: * ========= 047: * 048: * JOB (input) INTEGER 049: * = 1: an estimate for the largest singular value is computed. 050: * = 2: an estimate for the smallest singular value is computed. 051: * 052: * J (input) INTEGER 053: * Length of X and W 054: * 055: * X (input) REAL array, dimension (J) 056: * The j-vector x. 057: * 058: * SEST (input) REAL 059: * Estimated singular value of j by j matrix L 060: * 061: * W (input) REAL array, dimension (J) 062: * The j-vector w. 063: * 064: * GAMMA (input) REAL 065: * The diagonal element gamma. 066: * 067: * SESTPR (output) REAL 068: * Estimated singular value of (j+1) by (j+1) matrix Lhat. 069: * 070: * S (output) REAL 071: * Sine needed in forming xhat. 072: * 073: * C (output) REAL 074: * Cosine needed in forming xhat. 075: * 076: * ===================================================================== 077: * 078: * .. Parameters .. 079: REAL ZERO, ONE, TWO 080: PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 ) 081: REAL HALF, FOUR 082: PARAMETER ( HALF = 0.5E0, FOUR = 4.0E0 ) 083: * .. 084: * .. Local Scalars .. 085: REAL ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS, 086: $ NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2 087: * .. 088: * .. Intrinsic Functions .. 089: INTRINSIC ABS, MAX, SIGN, SQRT 090: * .. 091: * .. External Functions .. 092: REAL SDOT, SLAMCH 093: EXTERNAL SDOT, SLAMCH 094: * .. 095: * .. Executable Statements .. 096: * 097: EPS = SLAMCH( 'Epsilon' ) 098: ALPHA = SDOT( J, X, 1, W, 1 ) 099: * 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 SLAIC1 292: * 293: END 294: