|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE SLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C ) 00002 * 00003 * -- LAPACK auxiliary routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 INTEGER J, JOB 00010 REAL C, GAMMA, S, SEST, SESTPR 00011 * .. 00012 * .. Array Arguments .. 00013 REAL W( J ), X( J ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * SLAIC1 applies one step of incremental condition estimation in 00020 * its simplest version: 00021 * 00022 * Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j 00023 * lower triangular matrix L, such that 00024 * twonorm(L*x) = sest 00025 * Then SLAIC1 computes sestpr, s, c such that 00026 * the vector 00027 * [ s*x ] 00028 * xhat = [ c ] 00029 * is an approximate singular vector of 00030 * [ L 0 ] 00031 * Lhat = [ w' gamma ] 00032 * in the sense that 00033 * twonorm(Lhat*xhat) = sestpr. 00034 * 00035 * Depending on JOB, an estimate for the largest or smallest singular 00036 * value is computed. 00037 * 00038 * Note that [s c]' and sestpr**2 is an eigenpair of the system 00039 * 00040 * diag(sest*sest, 0) + [alpha gamma] * [ alpha ] 00041 * [ gamma ] 00042 * 00043 * where alpha = x'*w. 00044 * 00045 * Arguments 00046 * ========= 00047 * 00048 * JOB (input) INTEGER 00049 * = 1: an estimate for the largest singular value is computed. 00050 * = 2: an estimate for the smallest singular value is computed. 00051 * 00052 * J (input) INTEGER 00053 * Length of X and W 00054 * 00055 * X (input) REAL array, dimension (J) 00056 * The j-vector x. 00057 * 00058 * SEST (input) REAL 00059 * Estimated singular value of j by j matrix L 00060 * 00061 * W (input) REAL array, dimension (J) 00062 * The j-vector w. 00063 * 00064 * GAMMA (input) REAL 00065 * The diagonal element gamma. 00066 * 00067 * SESTPR (output) REAL 00068 * Estimated singular value of (j+1) by (j+1) matrix Lhat. 00069 * 00070 * S (output) REAL 00071 * Sine needed in forming xhat. 00072 * 00073 * C (output) REAL 00074 * Cosine needed in forming xhat. 00075 * 00076 * ===================================================================== 00077 * 00078 * .. Parameters .. 00079 REAL ZERO, ONE, TWO 00080 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 ) 00081 REAL HALF, FOUR 00082 PARAMETER ( HALF = 0.5E0, FOUR = 4.0E0 ) 00083 * .. 00084 * .. Local Scalars .. 00085 REAL ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS, 00086 $ NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2 00087 * .. 00088 * .. Intrinsic Functions .. 00089 INTRINSIC ABS, MAX, SIGN, SQRT 00090 * .. 00091 * .. External Functions .. 00092 REAL SDOT, SLAMCH 00093 EXTERNAL SDOT, SLAMCH 00094 * .. 00095 * .. Executable Statements .. 00096 * 00097 EPS = SLAMCH( 'Epsilon' ) 00098 ALPHA = SDOT( J, X, 1, W, 1 ) 00099 * 00100 ABSALP = ABS( ALPHA ) 00101 ABSGAM = ABS( GAMMA ) 00102 ABSEST = ABS( SEST ) 00103 * 00104 IF( JOB.EQ.1 ) THEN 00105 * 00106 * Estimating largest singular value 00107 * 00108 * special cases 00109 * 00110 IF( SEST.EQ.ZERO ) THEN 00111 S1 = MAX( ABSGAM, ABSALP ) 00112 IF( S1.EQ.ZERO ) THEN 00113 S = ZERO 00114 C = ONE 00115 SESTPR = ZERO 00116 ELSE 00117 S = ALPHA / S1 00118 C = GAMMA / S1 00119 TMP = SQRT( S*S+C*C ) 00120 S = S / TMP 00121 C = C / TMP 00122 SESTPR = S1*TMP 00123 END IF 00124 RETURN 00125 ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN 00126 S = ONE 00127 C = ZERO 00128 TMP = MAX( ABSEST, ABSALP ) 00129 S1 = ABSEST / TMP 00130 S2 = ABSALP / TMP 00131 SESTPR = TMP*SQRT( S1*S1+S2*S2 ) 00132 RETURN 00133 ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN 00134 S1 = ABSGAM 00135 S2 = ABSEST 00136 IF( S1.LE.S2 ) THEN 00137 S = ONE 00138 C = ZERO 00139 SESTPR = S2 00140 ELSE 00141 S = ZERO 00142 C = ONE 00143 SESTPR = S1 00144 END IF 00145 RETURN 00146 ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN 00147 S1 = ABSGAM 00148 S2 = ABSALP 00149 IF( S1.LE.S2 ) THEN 00150 TMP = S1 / S2 00151 S = SQRT( ONE+TMP*TMP ) 00152 SESTPR = S2*S 00153 C = ( GAMMA / S2 ) / S 00154 S = SIGN( ONE, ALPHA ) / S 00155 ELSE 00156 TMP = S2 / S1 00157 C = SQRT( ONE+TMP*TMP ) 00158 SESTPR = S1*C 00159 S = ( ALPHA / S1 ) / C 00160 C = SIGN( ONE, GAMMA ) / C 00161 END IF 00162 RETURN 00163 ELSE 00164 * 00165 * normal case 00166 * 00167 ZETA1 = ALPHA / ABSEST 00168 ZETA2 = GAMMA / ABSEST 00169 * 00170 B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF 00171 C = ZETA1*ZETA1 00172 IF( B.GT.ZERO ) THEN 00173 T = C / ( B+SQRT( B*B+C ) ) 00174 ELSE 00175 T = SQRT( B*B+C ) - B 00176 END IF 00177 * 00178 SINE = -ZETA1 / T 00179 COSINE = -ZETA2 / ( ONE+T ) 00180 TMP = SQRT( SINE*SINE+COSINE*COSINE ) 00181 S = SINE / TMP 00182 C = COSINE / TMP 00183 SESTPR = SQRT( T+ONE )*ABSEST 00184 RETURN 00185 END IF 00186 * 00187 ELSE IF( JOB.EQ.2 ) THEN 00188 * 00189 * Estimating smallest singular value 00190 * 00191 * special cases 00192 * 00193 IF( SEST.EQ.ZERO ) THEN 00194 SESTPR = ZERO 00195 IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN 00196 SINE = ONE 00197 COSINE = ZERO 00198 ELSE 00199 SINE = -GAMMA 00200 COSINE = ALPHA 00201 END IF 00202 S1 = MAX( ABS( SINE ), ABS( COSINE ) ) 00203 S = SINE / S1 00204 C = COSINE / S1 00205 TMP = SQRT( S*S+C*C ) 00206 S = S / TMP 00207 C = C / TMP 00208 RETURN 00209 ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN 00210 S = ZERO 00211 C = ONE 00212 SESTPR = ABSGAM 00213 RETURN 00214 ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN 00215 S1 = ABSGAM 00216 S2 = ABSEST 00217 IF( S1.LE.S2 ) THEN 00218 S = ZERO 00219 C = ONE 00220 SESTPR = S1 00221 ELSE 00222 S = ONE 00223 C = ZERO 00224 SESTPR = S2 00225 END IF 00226 RETURN 00227 ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN 00228 S1 = ABSGAM 00229 S2 = ABSALP 00230 IF( S1.LE.S2 ) THEN 00231 TMP = S1 / S2 00232 C = SQRT( ONE+TMP*TMP ) 00233 SESTPR = ABSEST*( TMP / C ) 00234 S = -( GAMMA / S2 ) / C 00235 C = SIGN( ONE, ALPHA ) / C 00236 ELSE 00237 TMP = S2 / S1 00238 S = SQRT( ONE+TMP*TMP ) 00239 SESTPR = ABSEST / S 00240 C = ( ALPHA / S1 ) / S 00241 S = -SIGN( ONE, GAMMA ) / S 00242 END IF 00243 RETURN 00244 ELSE 00245 * 00246 * normal case 00247 * 00248 ZETA1 = ALPHA / ABSEST 00249 ZETA2 = GAMMA / ABSEST 00250 * 00251 NORMA = MAX( ONE+ZETA1*ZETA1+ABS( ZETA1*ZETA2 ), 00252 $ ABS( ZETA1*ZETA2 )+ZETA2*ZETA2 ) 00253 * 00254 * See if root is closer to zero or to ONE 00255 * 00256 TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 ) 00257 IF( TEST.GE.ZERO ) THEN 00258 * 00259 * root is close to zero, compute directly 00260 * 00261 B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF 00262 C = ZETA2*ZETA2 00263 T = C / ( B+SQRT( ABS( B*B-C ) ) ) 00264 SINE = ZETA1 / ( ONE-T ) 00265 COSINE = -ZETA2 / T 00266 SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST 00267 ELSE 00268 * 00269 * root is closer to ONE, shift by that amount 00270 * 00271 B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF 00272 C = ZETA1*ZETA1 00273 IF( B.GE.ZERO ) THEN 00274 T = -C / ( B+SQRT( B*B+C ) ) 00275 ELSE 00276 T = B - SQRT( B*B+C ) 00277 END IF 00278 SINE = -ZETA1 / T 00279 COSINE = -ZETA2 / ( ONE+T ) 00280 SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST 00281 END IF 00282 TMP = SQRT( SINE*SINE+COSINE*COSINE ) 00283 S = SINE / TMP 00284 C = COSINE / TMP 00285 RETURN 00286 * 00287 END IF 00288 END IF 00289 RETURN 00290 * 00291 * End of SLAIC1 00292 * 00293 END
1.7.3 
