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