[22] | 1 | *> \brief \b SLAIC1 |
---|
| 2 | * |
---|
| 3 | * =========== DOCUMENTATION =========== |
---|
| 4 | * |
---|
| 5 | * Online html documentation available at |
---|
| 6 | * http://www.netlib.org/lapack/explore-html/ |
---|
| 7 | * |
---|
| 8 | *> \htmlonly |
---|
| 9 | *> Download SLAIC1 + dependencies |
---|
| 10 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaic1.f"> |
---|
| 11 | *> [TGZ]</a> |
---|
| 12 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaic1.f"> |
---|
| 13 | *> [ZIP]</a> |
---|
| 14 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaic1.f"> |
---|
| 15 | *> [TXT]</a> |
---|
| 16 | *> \endhtmlonly |
---|
| 17 | * |
---|
| 18 | * Definition: |
---|
| 19 | * =========== |
---|
| 20 | * |
---|
| 21 | * SUBROUTINE SLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C ) |
---|
| 22 | * |
---|
| 23 | * .. Scalar Arguments .. |
---|
| 24 | * INTEGER J, JOB |
---|
| 25 | * REAL C, GAMMA, S, SEST, SESTPR |
---|
| 26 | * .. |
---|
| 27 | * .. Array Arguments .. |
---|
| 28 | * REAL W( J ), X( J ) |
---|
| 29 | * .. |
---|
| 30 | * |
---|
| 31 | * |
---|
| 32 | *> \par Purpose: |
---|
| 33 | * ============= |
---|
| 34 | *> |
---|
| 35 | *> \verbatim |
---|
| 36 | *> |
---|
| 37 | *> SLAIC1 applies one step of incremental condition estimation in |
---|
| 38 | *> its simplest version: |
---|
| 39 | *> |
---|
| 40 | *> Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j |
---|
| 41 | *> lower triangular matrix L, such that |
---|
| 42 | *> twonorm(L*x) = sest |
---|
| 43 | *> Then SLAIC1 computes sestpr, s, c such that |
---|
| 44 | *> the vector |
---|
| 45 | *> [ s*x ] |
---|
| 46 | *> xhat = [ c ] |
---|
| 47 | *> is an approximate singular vector of |
---|
| 48 | *> [ L 0 ] |
---|
| 49 | *> Lhat = [ w**T gamma ] |
---|
| 50 | *> in the sense that |
---|
| 51 | *> twonorm(Lhat*xhat) = sestpr. |
---|
| 52 | *> |
---|
| 53 | *> Depending on JOB, an estimate for the largest or smallest singular |
---|
| 54 | *> value is computed. |
---|
| 55 | *> |
---|
| 56 | *> Note that [s c]**T and sestpr**2 is an eigenpair of the system |
---|
| 57 | *> |
---|
| 58 | *> diag(sest*sest, 0) + [alpha gamma] * [ alpha ] |
---|
| 59 | *> [ gamma ] |
---|
| 60 | *> |
---|
| 61 | *> where alpha = x**T*w. |
---|
| 62 | *> \endverbatim |
---|
| 63 | * |
---|
| 64 | * Arguments: |
---|
| 65 | * ========== |
---|
| 66 | * |
---|
| 67 | *> \param[in] JOB |
---|
| 68 | *> \verbatim |
---|
| 69 | *> JOB is INTEGER |
---|
| 70 | *> = 1: an estimate for the largest singular value is computed. |
---|
| 71 | *> = 2: an estimate for the smallest singular value is computed. |
---|
| 72 | *> \endverbatim |
---|
| 73 | *> |
---|
| 74 | *> \param[in] J |
---|
| 75 | *> \verbatim |
---|
| 76 | *> J is INTEGER |
---|
| 77 | *> Length of X and W |
---|
| 78 | *> \endverbatim |
---|
| 79 | *> |
---|
| 80 | *> \param[in] X |
---|
| 81 | *> \verbatim |
---|
| 82 | *> X is REAL array, dimension (J) |
---|
| 83 | *> The j-vector x. |
---|
| 84 | *> \endverbatim |
---|
| 85 | *> |
---|
| 86 | *> \param[in] SEST |
---|
| 87 | *> \verbatim |
---|
| 88 | *> SEST is REAL |
---|
| 89 | *> Estimated singular value of j by j matrix L |
---|
| 90 | *> \endverbatim |
---|
| 91 | *> |
---|
| 92 | *> \param[in] W |
---|
| 93 | *> \verbatim |
---|
| 94 | *> W is REAL array, dimension (J) |
---|
| 95 | *> The j-vector w. |
---|
| 96 | *> \endverbatim |
---|
| 97 | *> |
---|
| 98 | *> \param[in] GAMMA |
---|
| 99 | *> \verbatim |
---|
| 100 | *> GAMMA is REAL |
---|
| 101 | *> The diagonal element gamma. |
---|
| 102 | *> \endverbatim |
---|
| 103 | *> |
---|
| 104 | *> \param[out] SESTPR |
---|
| 105 | *> \verbatim |
---|
| 106 | *> SESTPR is REAL |
---|
| 107 | *> Estimated singular value of (j+1) by (j+1) matrix Lhat. |
---|
| 108 | *> \endverbatim |
---|
| 109 | *> |
---|
| 110 | *> \param[out] S |
---|
| 111 | *> \verbatim |
---|
| 112 | *> S is REAL |
---|
| 113 | *> Sine needed in forming xhat. |
---|
| 114 | *> \endverbatim |
---|
| 115 | *> |
---|
| 116 | *> \param[out] C |
---|
| 117 | *> \verbatim |
---|
| 118 | *> C is REAL |
---|
| 119 | *> Cosine needed in forming xhat. |
---|
| 120 | *> \endverbatim |
---|
| 121 | * |
---|
| 122 | * Authors: |
---|
| 123 | * ======== |
---|
| 124 | * |
---|
| 125 | *> \author Univ. of Tennessee |
---|
| 126 | *> \author Univ. of California Berkeley |
---|
| 127 | *> \author Univ. of Colorado Denver |
---|
| 128 | *> \author NAG Ltd. |
---|
| 129 | * |
---|
| 130 | *> \date November 2011 |
---|
| 131 | * |
---|
| 132 | *> \ingroup realOTHERauxiliary |
---|
| 133 | * |
---|
| 134 | * ===================================================================== |
---|
| 135 | SUBROUTINE SLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C ) |
---|
| 136 | * |
---|
| 137 | * -- LAPACK auxiliary routine (version 3.4.0) -- |
---|
| 138 | * -- LAPACK is a software package provided by Univ. of Tennessee, -- |
---|
| 139 | * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
---|
| 140 | * November 2011 |
---|
| 141 | * |
---|
| 142 | * .. Scalar Arguments .. |
---|
| 143 | INTEGER J, JOB |
---|
| 144 | REAL C, GAMMA, S, SEST, SESTPR |
---|
| 145 | * .. |
---|
| 146 | * .. Array Arguments .. |
---|
| 147 | REAL W( J ), X( J ) |
---|
| 148 | * .. |
---|
| 149 | * |
---|
| 150 | * ===================================================================== |
---|
| 151 | * |
---|
| 152 | * .. Parameters .. |
---|
| 153 | REAL ZERO, ONE, TWO |
---|
| 154 | PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 ) |
---|
| 155 | REAL HALF, FOUR |
---|
| 156 | PARAMETER ( HALF = 0.5E0, FOUR = 4.0E0 ) |
---|
| 157 | * .. |
---|
| 158 | * .. Local Scalars .. |
---|
| 159 | REAL ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS, |
---|
| 160 | $ NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2 |
---|
| 161 | * .. |
---|
| 162 | * .. Intrinsic Functions .. |
---|
| 163 | INTRINSIC ABS, MAX, SIGN, SQRT |
---|
| 164 | * .. |
---|
| 165 | * .. External Functions .. |
---|
| 166 | REAL SDOT, SLAMCH |
---|
| 167 | EXTERNAL SDOT, SLAMCH |
---|
| 168 | * .. |
---|
| 169 | * .. Executable Statements .. |
---|
| 170 | * |
---|
| 171 | EPS = SLAMCH( 'Epsilon' ) |
---|
| 172 | ALPHA = SDOT( J, X, 1, W, 1 ) |
---|
| 173 | * |
---|
| 174 | ABSALP = ABS( ALPHA ) |
---|
| 175 | ABSGAM = ABS( GAMMA ) |
---|
| 176 | ABSEST = ABS( SEST ) |
---|
| 177 | * |
---|
| 178 | IF( JOB.EQ.1 ) THEN |
---|
| 179 | * |
---|
| 180 | * Estimating largest singular value |
---|
| 181 | * |
---|
| 182 | * special cases |
---|
| 183 | * |
---|
| 184 | IF( SEST.EQ.ZERO ) THEN |
---|
| 185 | S1 = MAX( ABSGAM, ABSALP ) |
---|
| 186 | IF( S1.EQ.ZERO ) THEN |
---|
| 187 | S = ZERO |
---|
| 188 | C = ONE |
---|
| 189 | SESTPR = ZERO |
---|
| 190 | ELSE |
---|
| 191 | S = ALPHA / S1 |
---|
| 192 | C = GAMMA / S1 |
---|
| 193 | TMP = SQRT( S*S+C*C ) |
---|
| 194 | S = S / TMP |
---|
| 195 | C = C / TMP |
---|
| 196 | SESTPR = S1*TMP |
---|
| 197 | END IF |
---|
| 198 | RETURN |
---|
| 199 | ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN |
---|
| 200 | S = ONE |
---|
| 201 | C = ZERO |
---|
| 202 | TMP = MAX( ABSEST, ABSALP ) |
---|
| 203 | S1 = ABSEST / TMP |
---|
| 204 | S2 = ABSALP / TMP |
---|
| 205 | SESTPR = TMP*SQRT( S1*S1+S2*S2 ) |
---|
| 206 | RETURN |
---|
| 207 | ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN |
---|
| 208 | S1 = ABSGAM |
---|
| 209 | S2 = ABSEST |
---|
| 210 | IF( S1.LE.S2 ) THEN |
---|
| 211 | S = ONE |
---|
| 212 | C = ZERO |
---|
| 213 | SESTPR = S2 |
---|
| 214 | ELSE |
---|
| 215 | S = ZERO |
---|
| 216 | C = ONE |
---|
| 217 | SESTPR = S1 |
---|
| 218 | END IF |
---|
| 219 | RETURN |
---|
| 220 | ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN |
---|
| 221 | S1 = ABSGAM |
---|
| 222 | S2 = ABSALP |
---|
| 223 | IF( S1.LE.S2 ) THEN |
---|
| 224 | TMP = S1 / S2 |
---|
| 225 | S = SQRT( ONE+TMP*TMP ) |
---|
| 226 | SESTPR = S2*S |
---|
| 227 | C = ( GAMMA / S2 ) / S |
---|
| 228 | S = SIGN( ONE, ALPHA ) / S |
---|
| 229 | ELSE |
---|
| 230 | TMP = S2 / S1 |
---|
| 231 | C = SQRT( ONE+TMP*TMP ) |
---|
| 232 | SESTPR = S1*C |
---|
| 233 | S = ( ALPHA / S1 ) / C |
---|
| 234 | C = SIGN( ONE, GAMMA ) / C |
---|
| 235 | END IF |
---|
| 236 | RETURN |
---|
| 237 | ELSE |
---|
| 238 | * |
---|
| 239 | * normal case |
---|
| 240 | * |
---|
| 241 | ZETA1 = ALPHA / ABSEST |
---|
| 242 | ZETA2 = GAMMA / ABSEST |
---|
| 243 | * |
---|
| 244 | B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF |
---|
| 245 | C = ZETA1*ZETA1 |
---|
| 246 | IF( B.GT.ZERO ) THEN |
---|
| 247 | T = C / ( B+SQRT( B*B+C ) ) |
---|
| 248 | ELSE |
---|
| 249 | T = SQRT( B*B+C ) - B |
---|
| 250 | END IF |
---|
| 251 | * |
---|
| 252 | SINE = -ZETA1 / T |
---|
| 253 | COSINE = -ZETA2 / ( ONE+T ) |
---|
| 254 | TMP = SQRT( SINE*SINE+COSINE*COSINE ) |
---|
| 255 | S = SINE / TMP |
---|
| 256 | C = COSINE / TMP |
---|
| 257 | SESTPR = SQRT( T+ONE )*ABSEST |
---|
| 258 | RETURN |
---|
| 259 | END IF |
---|
| 260 | * |
---|
| 261 | ELSE IF( JOB.EQ.2 ) THEN |
---|
| 262 | * |
---|
| 263 | * Estimating smallest singular value |
---|
| 264 | * |
---|
| 265 | * special cases |
---|
| 266 | * |
---|
| 267 | IF( SEST.EQ.ZERO ) THEN |
---|
| 268 | SESTPR = ZERO |
---|
| 269 | IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN |
---|
| 270 | SINE = ONE |
---|
| 271 | COSINE = ZERO |
---|
| 272 | ELSE |
---|
| 273 | SINE = -GAMMA |
---|
| 274 | COSINE = ALPHA |
---|
| 275 | END IF |
---|
| 276 | S1 = MAX( ABS( SINE ), ABS( COSINE ) ) |
---|
| 277 | S = SINE / S1 |
---|
| 278 | C = COSINE / S1 |
---|
| 279 | TMP = SQRT( S*S+C*C ) |
---|
| 280 | S = S / TMP |
---|
| 281 | C = C / TMP |
---|
| 282 | RETURN |
---|
| 283 | ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN |
---|
| 284 | S = ZERO |
---|
| 285 | C = ONE |
---|
| 286 | SESTPR = ABSGAM |
---|
| 287 | RETURN |
---|
| 288 | ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN |
---|
| 289 | S1 = ABSGAM |
---|
| 290 | S2 = ABSEST |
---|
| 291 | IF( S1.LE.S2 ) THEN |
---|
| 292 | S = ZERO |
---|
| 293 | C = ONE |
---|
| 294 | SESTPR = S1 |
---|
| 295 | ELSE |
---|
| 296 | S = ONE |
---|
| 297 | C = ZERO |
---|
| 298 | SESTPR = S2 |
---|
| 299 | END IF |
---|
| 300 | RETURN |
---|
| 301 | ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN |
---|
| 302 | S1 = ABSGAM |
---|
| 303 | S2 = ABSALP |
---|
| 304 | IF( S1.LE.S2 ) THEN |
---|
| 305 | TMP = S1 / S2 |
---|
| 306 | C = SQRT( ONE+TMP*TMP ) |
---|
| 307 | SESTPR = ABSEST*( TMP / C ) |
---|
| 308 | S = -( GAMMA / S2 ) / C |
---|
| 309 | C = SIGN( ONE, ALPHA ) / C |
---|
| 310 | ELSE |
---|
| 311 | TMP = S2 / S1 |
---|
| 312 | S = SQRT( ONE+TMP*TMP ) |
---|
| 313 | SESTPR = ABSEST / S |
---|
| 314 | C = ( ALPHA / S1 ) / S |
---|
| 315 | S = -SIGN( ONE, GAMMA ) / S |
---|
| 316 | END IF |
---|
| 317 | RETURN |
---|
| 318 | ELSE |
---|
| 319 | * |
---|
| 320 | * normal case |
---|
| 321 | * |
---|
| 322 | ZETA1 = ALPHA / ABSEST |
---|
| 323 | ZETA2 = GAMMA / ABSEST |
---|
| 324 | * |
---|
| 325 | NORMA = MAX( ONE+ZETA1*ZETA1+ABS( ZETA1*ZETA2 ), |
---|
| 326 | $ ABS( ZETA1*ZETA2 )+ZETA2*ZETA2 ) |
---|
| 327 | * |
---|
| 328 | * See if root is closer to zero or to ONE |
---|
| 329 | * |
---|
| 330 | TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 ) |
---|
| 331 | IF( TEST.GE.ZERO ) THEN |
---|
| 332 | * |
---|
| 333 | * root is close to zero, compute directly |
---|
| 334 | * |
---|
| 335 | B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF |
---|
| 336 | C = ZETA2*ZETA2 |
---|
| 337 | T = C / ( B+SQRT( ABS( B*B-C ) ) ) |
---|
| 338 | SINE = ZETA1 / ( ONE-T ) |
---|
| 339 | COSINE = -ZETA2 / T |
---|
| 340 | SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST |
---|
| 341 | ELSE |
---|
| 342 | * |
---|
| 343 | * root is closer to ONE, shift by that amount |
---|
| 344 | * |
---|
| 345 | B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF |
---|
| 346 | C = ZETA1*ZETA1 |
---|
| 347 | IF( B.GE.ZERO ) THEN |
---|
| 348 | T = -C / ( B+SQRT( B*B+C ) ) |
---|
| 349 | ELSE |
---|
| 350 | T = B - SQRT( B*B+C ) |
---|
| 351 | END IF |
---|
| 352 | SINE = -ZETA1 / T |
---|
| 353 | COSINE = -ZETA2 / ( ONE+T ) |
---|
| 354 | SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST |
---|
| 355 | END IF |
---|
| 356 | TMP = SQRT( SINE*SINE+COSINE*COSINE ) |
---|
| 357 | S = SINE / TMP |
---|
| 358 | C = COSINE / TMP |
---|
| 359 | RETURN |
---|
| 360 | * |
---|
| 361 | END IF |
---|
| 362 | END IF |
---|
| 363 | RETURN |
---|
| 364 | * |
---|
| 365 | * End of SLAIC1 |
---|
| 366 | * |
---|
| 367 | END |
---|