Mercurial > octave
comparison libcruft/lapack/slazq4.f @ 7789:82be108cc558
First attempt at single precision tyeps
* * *
corrections to qrupdate single precision routines
* * *
prefer demotion to single over promotion to double
* * *
Add single precision support to log2 function
* * *
Trivial PROJECT file update
* * *
Cache optimized hermitian/transpose methods
* * *
Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Sun, 27 Apr 2008 22:34:17 +0200 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| 7788:45f5faba05a2 | 7789:82be108cc558 |
|---|---|
| 1 SUBROUTINE SLAZQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, | |
| 2 $ DN1, DN2, TAU, TTYPE, G ) | |
| 3 * | |
| 4 * -- LAPACK auxiliary routine (version 3.1) -- | |
| 5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. | |
| 6 * November 2006 | |
| 7 * | |
| 8 * .. Scalar Arguments .. | |
| 9 INTEGER I0, N0, N0IN, PP, TTYPE | |
| 10 REAL DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU | |
| 11 * .. | |
| 12 * .. Array Arguments .. | |
| 13 REAL Z( * ) | |
| 14 * .. | |
| 15 * | |
| 16 * Purpose | |
| 17 * ======= | |
| 18 * | |
| 19 * SLAZQ4 computes an approximation TAU to the smallest eigenvalue | |
| 20 * using values of d from the previous transform. | |
| 21 * | |
| 22 * I0 (input) INTEGER | |
| 23 * First index. | |
| 24 * | |
| 25 * N0 (input) INTEGER | |
| 26 * Last index. | |
| 27 * | |
| 28 * Z (input) REAL array, dimension ( 4*N ) | |
| 29 * Z holds the qd array. | |
| 30 * | |
| 31 * PP (input) INTEGER | |
| 32 * PP=0 for ping, PP=1 for pong. | |
| 33 * | |
| 34 * N0IN (input) INTEGER | |
| 35 * The value of N0 at start of EIGTEST. | |
| 36 * | |
| 37 * DMIN (input) REAL | |
| 38 * Minimum value of d. | |
| 39 * | |
| 40 * DMIN1 (input) REAL | |
| 41 * Minimum value of d, excluding D( N0 ). | |
| 42 * | |
| 43 * DMIN2 (input) REAL | |
| 44 * Minimum value of d, excluding D( N0 ) and D( N0-1 ). | |
| 45 * | |
| 46 * DN (input) REAL | |
| 47 * d(N) | |
| 48 * | |
| 49 * DN1 (input) REAL | |
| 50 * d(N-1) | |
| 51 * | |
| 52 * DN2 (input) REAL | |
| 53 * d(N-2) | |
| 54 * | |
| 55 * TAU (output) REAL | |
| 56 * This is the shift. | |
| 57 * | |
| 58 * TTYPE (output) INTEGER | |
| 59 * Shift type. | |
| 60 * | |
| 61 * G (input/output) REAL | |
| 62 * G is passed as an argument in order to save its value between | |
| 63 * calls to SLAZQ4 | |
| 64 * | |
| 65 * Further Details | |
| 66 * =============== | |
| 67 * CNST1 = 9/16 | |
| 68 * | |
| 69 * This is a thread safe version of SLASQ4, which passes G through the | |
| 70 * argument list in place of declaring G in a SAVE statment. | |
| 71 * | |
| 72 * ===================================================================== | |
| 73 * | |
| 74 * .. Parameters .. | |
| 75 REAL CNST1, CNST2, CNST3 | |
| 76 PARAMETER ( CNST1 = 0.5630E0, CNST2 = 1.010E0, | |
| 77 $ CNST3 = 1.050E0 ) | |
| 78 REAL QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD | |
| 79 PARAMETER ( QURTR = 0.250E0, THIRD = 0.3330E0, | |
| 80 $ HALF = 0.50E0, ZERO = 0.0E0, ONE = 1.0E0, | |
| 81 $ TWO = 2.0E0, HUNDRD = 100.0E0 ) | |
| 82 * .. | |
| 83 * .. Local Scalars .. | |
| 84 INTEGER I4, NN, NP | |
| 85 REAL A2, B1, B2, GAM, GAP1, GAP2, S | |
| 86 * .. | |
| 87 * .. Intrinsic Functions .. | |
| 88 INTRINSIC MAX, MIN, SQRT | |
| 89 * .. | |
| 90 * .. Executable Statements .. | |
| 91 * | |
| 92 * A negative DMIN forces the shift to take that absolute value | |
| 93 * TTYPE records the type of shift. | |
| 94 * | |
| 95 IF( DMIN.LE.ZERO ) THEN | |
| 96 TAU = -DMIN | |
| 97 TTYPE = -1 | |
| 98 RETURN | |
| 99 END IF | |
| 100 * | |
| 101 NN = 4*N0 + PP | |
| 102 IF( N0IN.EQ.N0 ) THEN | |
| 103 * | |
| 104 * No eigenvalues deflated. | |
| 105 * | |
| 106 IF( DMIN.EQ.DN .OR. DMIN.EQ.DN1 ) THEN | |
| 107 * | |
| 108 B1 = SQRT( Z( NN-3 ) )*SQRT( Z( NN-5 ) ) | |
| 109 B2 = SQRT( Z( NN-7 ) )*SQRT( Z( NN-9 ) ) | |
| 110 A2 = Z( NN-7 ) + Z( NN-5 ) | |
| 111 * | |
| 112 * Cases 2 and 3. | |
| 113 * | |
| 114 IF( DMIN.EQ.DN .AND. DMIN1.EQ.DN1 ) THEN | |
| 115 GAP2 = DMIN2 - A2 - DMIN2*QURTR | |
| 116 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2 ) THEN | |
| 117 GAP1 = A2 - DN - ( B2 / GAP2 )*B2 | |
| 118 ELSE | |
| 119 GAP1 = A2 - DN - ( B1+B2 ) | |
| 120 END IF | |
| 121 IF( GAP1.GT.ZERO .AND. GAP1.GT.B1 ) THEN | |
| 122 S = MAX( DN-( B1 / GAP1 )*B1, HALF*DMIN ) | |
| 123 TTYPE = -2 | |
| 124 ELSE | |
| 125 S = ZERO | |
| 126 IF( DN.GT.B1 ) | |
| 127 $ S = DN - B1 | |
| 128 IF( A2.GT.( B1+B2 ) ) | |
| 129 $ S = MIN( S, A2-( B1+B2 ) ) | |
| 130 S = MAX( S, THIRD*DMIN ) | |
| 131 TTYPE = -3 | |
| 132 END IF | |
| 133 ELSE | |
| 134 * | |
| 135 * Case 4. | |
| 136 * | |
| 137 TTYPE = -4 | |
| 138 S = QURTR*DMIN | |
| 139 IF( DMIN.EQ.DN ) THEN | |
| 140 GAM = DN | |
| 141 A2 = ZERO | |
| 142 IF( Z( NN-5 ) .GT. Z( NN-7 ) ) | |
| 143 $ RETURN | |
| 144 B2 = Z( NN-5 ) / Z( NN-7 ) | |
| 145 NP = NN - 9 | |
| 146 ELSE | |
| 147 NP = NN - 2*PP | |
| 148 B2 = Z( NP-2 ) | |
| 149 GAM = DN1 | |
| 150 IF( Z( NP-4 ) .GT. Z( NP-2 ) ) | |
| 151 $ RETURN | |
| 152 A2 = Z( NP-4 ) / Z( NP-2 ) | |
| 153 IF( Z( NN-9 ) .GT. Z( NN-11 ) ) | |
| 154 $ RETURN | |
| 155 B2 = Z( NN-9 ) / Z( NN-11 ) | |
| 156 NP = NN - 13 | |
| 157 END IF | |
| 158 * | |
| 159 * Approximate contribution to norm squared from I < NN-1. | |
| 160 * | |
| 161 A2 = A2 + B2 | |
| 162 DO 10 I4 = NP, 4*I0 - 1 + PP, -4 | |
| 163 IF( B2.EQ.ZERO ) | |
| 164 $ GO TO 20 | |
| 165 B1 = B2 | |
| 166 IF( Z( I4 ) .GT. Z( I4-2 ) ) | |
| 167 $ RETURN | |
| 168 B2 = B2*( Z( I4 ) / Z( I4-2 ) ) | |
| 169 A2 = A2 + B2 | |
| 170 IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) | |
| 171 $ GO TO 20 | |
| 172 10 CONTINUE | |
| 173 20 CONTINUE | |
| 174 A2 = CNST3*A2 | |
| 175 * | |
| 176 * Rayleigh quotient residual bound. | |
| 177 * | |
| 178 IF( A2.LT.CNST1 ) | |
| 179 $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 ) | |
| 180 END IF | |
| 181 ELSE IF( DMIN.EQ.DN2 ) THEN | |
| 182 * | |
| 183 * Case 5. | |
| 184 * | |
| 185 TTYPE = -5 | |
| 186 S = QURTR*DMIN | |
| 187 * | |
| 188 * Compute contribution to norm squared from I > NN-2. | |
| 189 * | |
| 190 NP = NN - 2*PP | |
| 191 B1 = Z( NP-2 ) | |
| 192 B2 = Z( NP-6 ) | |
| 193 GAM = DN2 | |
| 194 IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 ) | |
| 195 $ RETURN | |
| 196 A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 ) | |
| 197 * | |
| 198 * Approximate contribution to norm squared from I < NN-2. | |
| 199 * | |
| 200 IF( N0-I0.GT.2 ) THEN | |
| 201 B2 = Z( NN-13 ) / Z( NN-15 ) | |
| 202 A2 = A2 + B2 | |
| 203 DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4 | |
| 204 IF( B2.EQ.ZERO ) | |
| 205 $ GO TO 40 | |
| 206 B1 = B2 | |
| 207 IF( Z( I4 ) .GT. Z( I4-2 ) ) | |
| 208 $ RETURN | |
| 209 B2 = B2*( Z( I4 ) / Z( I4-2 ) ) | |
| 210 A2 = A2 + B2 | |
| 211 IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) | |
| 212 $ GO TO 40 | |
| 213 30 CONTINUE | |
| 214 40 CONTINUE | |
| 215 A2 = CNST3*A2 | |
| 216 END IF | |
| 217 * | |
| 218 IF( A2.LT.CNST1 ) | |
| 219 $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 ) | |
| 220 ELSE | |
| 221 * | |
| 222 * Case 6, no information to guide us. | |
| 223 * | |
| 224 IF( TTYPE.EQ.-6 ) THEN | |
| 225 G = G + THIRD*( ONE-G ) | |
| 226 ELSE IF( TTYPE.EQ.-18 ) THEN | |
| 227 G = QURTR*THIRD | |
| 228 ELSE | |
| 229 G = QURTR | |
| 230 END IF | |
| 231 S = G*DMIN | |
| 232 TTYPE = -6 | |
| 233 END IF | |
| 234 * | |
| 235 ELSE IF( N0IN.EQ.( N0+1 ) ) THEN | |
| 236 * | |
| 237 * One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN. | |
| 238 * | |
| 239 IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN | |
| 240 * | |
| 241 * Cases 7 and 8. | |
| 242 * | |
| 243 TTYPE = -7 | |
| 244 S = THIRD*DMIN1 | |
| 245 IF( Z( NN-5 ).GT.Z( NN-7 ) ) | |
| 246 $ RETURN | |
| 247 B1 = Z( NN-5 ) / Z( NN-7 ) | |
| 248 B2 = B1 | |
| 249 IF( B2.EQ.ZERO ) | |
| 250 $ GO TO 60 | |
| 251 DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4 | |
| 252 A2 = B1 | |
| 253 IF( Z( I4 ).GT.Z( I4-2 ) ) | |
| 254 $ RETURN | |
| 255 B1 = B1*( Z( I4 ) / Z( I4-2 ) ) | |
| 256 B2 = B2 + B1 | |
| 257 IF( HUNDRD*MAX( B1, A2 ).LT.B2 ) | |
| 258 $ GO TO 60 | |
| 259 50 CONTINUE | |
| 260 60 CONTINUE | |
| 261 B2 = SQRT( CNST3*B2 ) | |
| 262 A2 = DMIN1 / ( ONE+B2**2 ) | |
| 263 GAP2 = HALF*DMIN2 - A2 | |
| 264 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN | |
| 265 S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) ) | |
| 266 ELSE | |
| 267 S = MAX( S, A2*( ONE-CNST2*B2 ) ) | |
| 268 TTYPE = -8 | |
| 269 END IF | |
| 270 ELSE | |
| 271 * | |
| 272 * Case 9. | |
| 273 * | |
| 274 S = QURTR*DMIN1 | |
| 275 IF( DMIN1.EQ.DN1 ) | |
| 276 $ S = HALF*DMIN1 | |
| 277 TTYPE = -9 | |
| 278 END IF | |
| 279 * | |
| 280 ELSE IF( N0IN.EQ.( N0+2 ) ) THEN | |
| 281 * | |
| 282 * Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN. | |
| 283 * | |
| 284 * Cases 10 and 11. | |
| 285 * | |
| 286 IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN | |
| 287 TTYPE = -10 | |
| 288 S = THIRD*DMIN2 | |
| 289 IF( Z( NN-5 ).GT.Z( NN-7 ) ) | |
| 290 $ RETURN | |
| 291 B1 = Z( NN-5 ) / Z( NN-7 ) | |
| 292 B2 = B1 | |
| 293 IF( B2.EQ.ZERO ) | |
| 294 $ GO TO 80 | |
| 295 DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4 | |
| 296 IF( Z( I4 ).GT.Z( I4-2 ) ) | |
| 297 $ RETURN | |
| 298 B1 = B1*( Z( I4 ) / Z( I4-2 ) ) | |
| 299 B2 = B2 + B1 | |
| 300 IF( HUNDRD*B1.LT.B2 ) | |
| 301 $ GO TO 80 | |
| 302 70 CONTINUE | |
| 303 80 CONTINUE | |
| 304 B2 = SQRT( CNST3*B2 ) | |
| 305 A2 = DMIN2 / ( ONE+B2**2 ) | |
| 306 GAP2 = Z( NN-7 ) + Z( NN-9 ) - | |
| 307 $ SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2 | |
| 308 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN | |
| 309 S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) ) | |
| 310 ELSE | |
| 311 S = MAX( S, A2*( ONE-CNST2*B2 ) ) | |
| 312 END IF | |
| 313 ELSE | |
| 314 S = QURTR*DMIN2 | |
| 315 TTYPE = -11 | |
| 316 END IF | |
| 317 ELSE IF( N0IN.GT.( N0+2 ) ) THEN | |
| 318 * | |
| 319 * Case 12, more than two eigenvalues deflated. No information. | |
| 320 * | |
| 321 S = ZERO | |
| 322 TTYPE = -12 | |
| 323 END IF | |
| 324 * | |
| 325 TAU = S | |
| 326 RETURN | |
| 327 * | |
| 328 * End of SLAZQ4 | |
| 329 * | |
| 330 END |