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