![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLAIC1 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLAIC1 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaic1.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaic1.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaic1.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER J, JOB 00025 * REAL C, GAMMA, S, SEST, SESTPR 00026 * .. 00027 * .. Array Arguments .. 00028 * REAL W( J ), X( J ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> SLAIC1 applies one step of incremental condition estimation in 00038 *> its simplest version: 00039 *> 00040 *> Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j 00041 *> lower triangular matrix L, such that 00042 *> twonorm(L*x) = sest 00043 *> Then SLAIC1 computes sestpr, s, c such that 00044 *> the vector 00045 *> [ s*x ] 00046 *> xhat = [ c ] 00047 *> is an approximate singular vector of 00048 *> [ L 0 ] 00049 *> Lhat = [ w**T gamma ] 00050 *> in the sense that 00051 *> twonorm(Lhat*xhat) = sestpr. 00052 *> 00053 *> Depending on JOB, an estimate for the largest or smallest singular 00054 *> value is computed. 00055 *> 00056 *> Note that [s c]**T and sestpr**2 is an eigenpair of the system 00057 *> 00058 *> diag(sest*sest, 0) + [alpha gamma] * [ alpha ] 00059 *> [ gamma ] 00060 *> 00061 *> where alpha = x**T*w. 00062 *> \endverbatim 00063 * 00064 * Arguments: 00065 * ========== 00066 * 00067 *> \param[in] JOB 00068 *> \verbatim 00069 *> JOB is INTEGER 00070 *> = 1: an estimate for the largest singular value is computed. 00071 *> = 2: an estimate for the smallest singular value is computed. 00072 *> \endverbatim 00073 *> 00074 *> \param[in] J 00075 *> \verbatim 00076 *> J is INTEGER 00077 *> Length of X and W 00078 *> \endverbatim 00079 *> 00080 *> \param[in] X 00081 *> \verbatim 00082 *> X is REAL array, dimension (J) 00083 *> The j-vector x. 00084 *> \endverbatim 00085 *> 00086 *> \param[in] SEST 00087 *> \verbatim 00088 *> SEST is REAL 00089 *> Estimated singular value of j by j matrix L 00090 *> \endverbatim 00091 *> 00092 *> \param[in] W 00093 *> \verbatim 00094 *> W is REAL array, dimension (J) 00095 *> The j-vector w. 00096 *> \endverbatim 00097 *> 00098 *> \param[in] GAMMA 00099 *> \verbatim 00100 *> GAMMA is REAL 00101 *> The diagonal element gamma. 00102 *> \endverbatim 00103 *> 00104 *> \param[out] SESTPR 00105 *> \verbatim 00106 *> SESTPR is REAL 00107 *> Estimated singular value of (j+1) by (j+1) matrix Lhat. 00108 *> \endverbatim 00109 *> 00110 *> \param[out] S 00111 *> \verbatim 00112 *> S is REAL 00113 *> Sine needed in forming xhat. 00114 *> \endverbatim 00115 *> 00116 *> \param[out] C 00117 *> \verbatim 00118 *> C is REAL 00119 *> Cosine needed in forming xhat. 00120 *> \endverbatim 00121 * 00122 * Authors: 00123 * ======== 00124 * 00125 *> \author Univ. of Tennessee 00126 *> \author Univ. of California Berkeley 00127 *> \author Univ. of Colorado Denver 00128 *> \author NAG Ltd. 00129 * 00130 *> \date November 2011 00131 * 00132 *> \ingroup realOTHERauxiliary 00133 * 00134 * ===================================================================== 00135 SUBROUTINE SLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C ) 00136 * 00137 * -- LAPACK auxiliary routine (version 3.4.0) -- 00138 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00139 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00140 * November 2011 00141 * 00142 * .. Scalar Arguments .. 00143 INTEGER J, JOB 00144 REAL C, GAMMA, S, SEST, SESTPR 00145 * .. 00146 * .. Array Arguments .. 00147 REAL W( J ), X( J ) 00148 * .. 00149 * 00150 * ===================================================================== 00151 * 00152 * .. Parameters .. 00153 REAL ZERO, ONE, TWO 00154 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 ) 00155 REAL HALF, FOUR 00156 PARAMETER ( HALF = 0.5E0, FOUR = 4.0E0 ) 00157 * .. 00158 * .. Local Scalars .. 00159 REAL ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS, 00160 $ NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2 00161 * .. 00162 * .. Intrinsic Functions .. 00163 INTRINSIC ABS, MAX, SIGN, SQRT 00164 * .. 00165 * .. External Functions .. 00166 REAL SDOT, SLAMCH 00167 EXTERNAL SDOT, SLAMCH 00168 * .. 00169 * .. Executable Statements .. 00170 * 00171 EPS = SLAMCH( 'Epsilon' ) 00172 ALPHA = SDOT( J, X, 1, W, 1 ) 00173 * 00174 ABSALP = ABS( ALPHA ) 00175 ABSGAM = ABS( GAMMA ) 00176 ABSEST = ABS( SEST ) 00177 * 00178 IF( JOB.EQ.1 ) THEN 00179 * 00180 * Estimating largest singular value 00181 * 00182 * special cases 00183 * 00184 IF( SEST.EQ.ZERO ) THEN 00185 S1 = MAX( ABSGAM, ABSALP ) 00186 IF( S1.EQ.ZERO ) THEN 00187 S = ZERO 00188 C = ONE 00189 SESTPR = ZERO 00190 ELSE 00191 S = ALPHA / S1 00192 C = GAMMA / S1 00193 TMP = SQRT( S*S+C*C ) 00194 S = S / TMP 00195 C = C / TMP 00196 SESTPR = S1*TMP 00197 END IF 00198 RETURN 00199 ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN 00200 S = ONE 00201 C = ZERO 00202 TMP = MAX( ABSEST, ABSALP ) 00203 S1 = ABSEST / TMP 00204 S2 = ABSALP / TMP 00205 SESTPR = TMP*SQRT( S1*S1+S2*S2 ) 00206 RETURN 00207 ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN 00208 S1 = ABSGAM 00209 S2 = ABSEST 00210 IF( S1.LE.S2 ) THEN 00211 S = ONE 00212 C = ZERO 00213 SESTPR = S2 00214 ELSE 00215 S = ZERO 00216 C = ONE 00217 SESTPR = S1 00218 END IF 00219 RETURN 00220 ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN 00221 S1 = ABSGAM 00222 S2 = ABSALP 00223 IF( S1.LE.S2 ) THEN 00224 TMP = S1 / S2 00225 S = SQRT( ONE+TMP*TMP ) 00226 SESTPR = S2*S 00227 C = ( GAMMA / S2 ) / S 00228 S = SIGN( ONE, ALPHA ) / S 00229 ELSE 00230 TMP = S2 / S1 00231 C = SQRT( ONE+TMP*TMP ) 00232 SESTPR = S1*C 00233 S = ( ALPHA / S1 ) / C 00234 C = SIGN( ONE, GAMMA ) / C 00235 END IF 00236 RETURN 00237 ELSE 00238 * 00239 * normal case 00240 * 00241 ZETA1 = ALPHA / ABSEST 00242 ZETA2 = GAMMA / ABSEST 00243 * 00244 B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF 00245 C = ZETA1*ZETA1 00246 IF( B.GT.ZERO ) THEN 00247 T = C / ( B+SQRT( B*B+C ) ) 00248 ELSE 00249 T = SQRT( B*B+C ) - B 00250 END IF 00251 * 00252 SINE = -ZETA1 / T 00253 COSINE = -ZETA2 / ( ONE+T ) 00254 TMP = SQRT( SINE*SINE+COSINE*COSINE ) 00255 S = SINE / TMP 00256 C = COSINE / TMP 00257 SESTPR = SQRT( T+ONE )*ABSEST 00258 RETURN 00259 END IF 00260 * 00261 ELSE IF( JOB.EQ.2 ) THEN 00262 * 00263 * Estimating smallest singular value 00264 * 00265 * special cases 00266 * 00267 IF( SEST.EQ.ZERO ) THEN 00268 SESTPR = ZERO 00269 IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN 00270 SINE = ONE 00271 COSINE = ZERO 00272 ELSE 00273 SINE = -GAMMA 00274 COSINE = ALPHA 00275 END IF 00276 S1 = MAX( ABS( SINE ), ABS( COSINE ) ) 00277 S = SINE / S1 00278 C = COSINE / S1 00279 TMP = SQRT( S*S+C*C ) 00280 S = S / TMP 00281 C = C / TMP 00282 RETURN 00283 ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN 00284 S = ZERO 00285 C = ONE 00286 SESTPR = ABSGAM 00287 RETURN 00288 ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN 00289 S1 = ABSGAM 00290 S2 = ABSEST 00291 IF( S1.LE.S2 ) THEN 00292 S = ZERO 00293 C = ONE 00294 SESTPR = S1 00295 ELSE 00296 S = ONE 00297 C = ZERO 00298 SESTPR = S2 00299 END IF 00300 RETURN 00301 ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN 00302 S1 = ABSGAM 00303 S2 = ABSALP 00304 IF( S1.LE.S2 ) THEN 00305 TMP = S1 / S2 00306 C = SQRT( ONE+TMP*TMP ) 00307 SESTPR = ABSEST*( TMP / C ) 00308 S = -( GAMMA / S2 ) / C 00309 C = SIGN( ONE, ALPHA ) / C 00310 ELSE 00311 TMP = S2 / S1 00312 S = SQRT( ONE+TMP*TMP ) 00313 SESTPR = ABSEST / S 00314 C = ( ALPHA / S1 ) / S 00315 S = -SIGN( ONE, GAMMA ) / S 00316 END IF 00317 RETURN 00318 ELSE 00319 * 00320 * normal case 00321 * 00322 ZETA1 = ALPHA / ABSEST 00323 ZETA2 = GAMMA / ABSEST 00324 * 00325 NORMA = MAX( ONE+ZETA1*ZETA1+ABS( ZETA1*ZETA2 ), 00326 $ ABS( ZETA1*ZETA2 )+ZETA2*ZETA2 ) 00327 * 00328 * See if root is closer to zero or to ONE 00329 * 00330 TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 ) 00331 IF( TEST.GE.ZERO ) THEN 00332 * 00333 * root is close to zero, compute directly 00334 * 00335 B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF 00336 C = ZETA2*ZETA2 00337 T = C / ( B+SQRT( ABS( B*B-C ) ) ) 00338 SINE = ZETA1 / ( ONE-T ) 00339 COSINE = -ZETA2 / T 00340 SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST 00341 ELSE 00342 * 00343 * root is closer to ONE, shift by that amount 00344 * 00345 B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF 00346 C = ZETA1*ZETA1 00347 IF( B.GE.ZERO ) THEN 00348 T = -C / ( B+SQRT( B*B+C ) ) 00349 ELSE 00350 T = B - SQRT( B*B+C ) 00351 END IF 00352 SINE = -ZETA1 / T 00353 COSINE = -ZETA2 / ( ONE+T ) 00354 SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST 00355 END IF 00356 TMP = SQRT( SINE*SINE+COSINE*COSINE ) 00357 S = SINE / TMP 00358 C = COSINE / TMP 00359 RETURN 00360 * 00361 END IF 00362 END IF 00363 RETURN 00364 * 00365 * End of SLAIC1 00366 * 00367 END