 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, 00002 $ DN1, DN2, TAU, TTYPE, G ) 00003 * 00004 * -- LAPACK routine (version 3.3.1) -- 00005 * 00006 * -- Contributed by Osni Marques of the Lawrence Berkeley National -- 00007 * -- Laboratory and Beresford Parlett of the Univ. of California at -- 00008 * -- Berkeley -- 00009 * -- November 2008 -- 00010 * 00011 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00012 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00013 * 00014 * .. Scalar Arguments .. 00015 INTEGER I0, N0, N0IN, PP, TTYPE 00016 DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU 00017 * .. 00018 * .. Array Arguments .. 00019 DOUBLE PRECISION Z( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * DLASQ4 computes an approximation TAU to the smallest eigenvalue 00026 * using values of d from the previous transform. 00027 * 00028 * Arguments 00029 * ========= 00030 * 00031 * I0 (input) INTEGER 00032 * First index. 00033 * 00034 * N0 (input) INTEGER 00035 * Last index. 00036 * 00037 * Z (input) DOUBLE PRECISION array, dimension ( 4*N ) 00038 * Z holds the qd array. 00039 * 00040 * PP (input) INTEGER 00041 * PP=0 for ping, PP=1 for pong. 00042 * 00043 * NOIN (input) INTEGER 00044 * The value of N0 at start of EIGTEST. 00045 * 00046 * DMIN (input) DOUBLE PRECISION 00047 * Minimum value of d. 00048 * 00049 * DMIN1 (input) DOUBLE PRECISION 00050 * Minimum value of d, excluding D( N0 ). 00051 * 00052 * DMIN2 (input) DOUBLE PRECISION 00053 * Minimum value of d, excluding D( N0 ) and D( N0-1 ). 00054 * 00055 * DN (input) DOUBLE PRECISION 00056 * d(N) 00057 * 00058 * DN1 (input) DOUBLE PRECISION 00059 * d(N-1) 00060 * 00061 * DN2 (input) DOUBLE PRECISION 00062 * d(N-2) 00063 * 00064 * TAU (output) DOUBLE PRECISION 00065 * This is the shift. 00066 * 00067 * TTYPE (output) INTEGER 00068 * Shift type. 00069 * 00070 * G (input/output) REAL 00071 * G is passed as an argument in order to save its value between 00072 * calls to DLASQ4. 00073 * 00074 * Further Details 00075 * =============== 00076 * CNST1 = 9/16 00077 * 00078 * ===================================================================== 00079 * 00080 * .. Parameters .. 00081 DOUBLE PRECISION CNST1, CNST2, CNST3 00082 PARAMETER ( CNST1 = 0.5630D0, CNST2 = 1.010D0, 00083 $ CNST3 = 1.050D0 ) 00084 DOUBLE PRECISION QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD 00085 PARAMETER ( QURTR = 0.250D0, THIRD = 0.3330D0, 00086 $ HALF = 0.50D0, ZERO = 0.0D0, ONE = 1.0D0, 00087 $ TWO = 2.0D0, HUNDRD = 100.0D0 ) 00088 * .. 00089 * .. Local Scalars .. 00090 INTEGER I4, NN, NP 00091 DOUBLE PRECISION A2, B1, B2, GAM, GAP1, GAP2, S 00092 * .. 00093 * .. Intrinsic Functions .. 00094 INTRINSIC MAX, MIN, SQRT 00095 * .. 00096 * .. Executable Statements .. 00097 * 00098 * A negative DMIN forces the shift to take that absolute value 00099 * TTYPE records the type of shift. 00100 * 00101 IF( DMIN.LE.ZERO ) THEN 00102 TAU = -DMIN 00103 TTYPE = -1 00104 RETURN 00105 END IF 00106 * 00107 NN = 4*N0 + PP 00108 IF( N0IN.EQ.N0 ) THEN 00109 * 00110 * No eigenvalues deflated. 00111 * 00112 IF( DMIN.EQ.DN .OR. DMIN.EQ.DN1 ) THEN 00113 * 00114 B1 = SQRT( Z( NN-3 ) )*SQRT( Z( NN-5 ) ) 00115 B2 = SQRT( Z( NN-7 ) )*SQRT( Z( NN-9 ) ) 00116 A2 = Z( NN-7 ) + Z( NN-5 ) 00117 * 00118 * Cases 2 and 3. 00119 * 00120 IF( DMIN.EQ.DN .AND. DMIN1.EQ.DN1 ) THEN 00121 GAP2 = DMIN2 - A2 - DMIN2*QURTR 00122 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2 ) THEN 00123 GAP1 = A2 - DN - ( B2 / GAP2 )*B2 00124 ELSE 00125 GAP1 = A2 - DN - ( B1+B2 ) 00126 END IF 00127 IF( GAP1.GT.ZERO .AND. GAP1.GT.B1 ) THEN 00128 S = MAX( DN-( B1 / GAP1 )*B1, HALF*DMIN ) 00129 TTYPE = -2 00130 ELSE 00131 S = ZERO 00132 IF( DN.GT.B1 ) 00133 $ S = DN - B1 00134 IF( A2.GT.( B1+B2 ) ) 00135 $ S = MIN( S, A2-( B1+B2 ) ) 00136 S = MAX( S, THIRD*DMIN ) 00137 TTYPE = -3 00138 END IF 00139 ELSE 00140 * 00141 * Case 4. 00142 * 00143 TTYPE = -4 00144 S = QURTR*DMIN 00145 IF( DMIN.EQ.DN ) THEN 00146 GAM = DN 00147 A2 = ZERO 00148 IF( Z( NN-5 ) .GT. Z( NN-7 ) ) 00149 $ RETURN 00150 B2 = Z( NN-5 ) / Z( NN-7 ) 00151 NP = NN - 9 00152 ELSE 00153 NP = NN - 2*PP 00154 B2 = Z( NP-2 ) 00155 GAM = DN1 00156 IF( Z( NP-4 ) .GT. Z( NP-2 ) ) 00157 $ RETURN 00158 A2 = Z( NP-4 ) / Z( NP-2 ) 00159 IF( Z( NN-9 ) .GT. Z( NN-11 ) ) 00160 $ RETURN 00161 B2 = Z( NN-9 ) / Z( NN-11 ) 00162 NP = NN - 13 00163 END IF 00164 * 00165 * Approximate contribution to norm squared from I < NN-1. 00166 * 00167 A2 = A2 + B2 00168 DO 10 I4 = NP, 4*I0 - 1 + PP, -4 00169 IF( B2.EQ.ZERO ) 00170 $ GO TO 20 00171 B1 = B2 00172 IF( Z( I4 ) .GT. Z( I4-2 ) ) 00173 $ RETURN 00174 B2 = B2*( Z( I4 ) / Z( I4-2 ) ) 00175 A2 = A2 + B2 00176 IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) 00177 $ GO TO 20 00178 10 CONTINUE 00179 20 CONTINUE 00180 A2 = CNST3*A2 00181 * 00182 * Rayleigh quotient residual bound. 00183 * 00184 IF( A2.LT.CNST1 ) 00185 $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 ) 00186 END IF 00187 ELSE IF( DMIN.EQ.DN2 ) THEN 00188 * 00189 * Case 5. 00190 * 00191 TTYPE = -5 00192 S = QURTR*DMIN 00193 * 00194 * Compute contribution to norm squared from I > NN-2. 00195 * 00196 NP = NN - 2*PP 00197 B1 = Z( NP-2 ) 00198 B2 = Z( NP-6 ) 00199 GAM = DN2 00200 IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 ) 00201 $ RETURN 00202 A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 ) 00203 * 00204 * Approximate contribution to norm squared from I < NN-2. 00205 * 00206 IF( N0-I0.GT.2 ) THEN 00207 B2 = Z( NN-13 ) / Z( NN-15 ) 00208 A2 = A2 + B2 00209 DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4 00210 IF( B2.EQ.ZERO ) 00211 $ GO TO 40 00212 B1 = B2 00213 IF( Z( I4 ) .GT. Z( I4-2 ) ) 00214 $ RETURN 00215 B2 = B2*( Z( I4 ) / Z( I4-2 ) ) 00216 A2 = A2 + B2 00217 IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) 00218 $ GO TO 40 00219 30 CONTINUE 00220 40 CONTINUE 00221 A2 = CNST3*A2 00222 END IF 00223 * 00224 IF( A2.LT.CNST1 ) 00225 $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 ) 00226 ELSE 00227 * 00228 * Case 6, no information to guide us. 00229 * 00230 IF( TTYPE.EQ.-6 ) THEN 00231 G = G + THIRD*( ONE-G ) 00232 ELSE IF( TTYPE.EQ.-18 ) THEN 00233 G = QURTR*THIRD 00234 ELSE 00235 G = QURTR 00236 END IF 00237 S = G*DMIN 00238 TTYPE = -6 00239 END IF 00240 * 00241 ELSE IF( N0IN.EQ.( N0+1 ) ) THEN 00242 * 00243 * One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN. 00244 * 00245 IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN 00246 * 00247 * Cases 7 and 8. 00248 * 00249 TTYPE = -7 00250 S = THIRD*DMIN1 00251 IF( Z( NN-5 ).GT.Z( NN-7 ) ) 00252 $ RETURN 00253 B1 = Z( NN-5 ) / Z( NN-7 ) 00254 B2 = B1 00255 IF( B2.EQ.ZERO ) 00256 $ GO TO 60 00257 DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4 00258 A2 = B1 00259 IF( Z( I4 ).GT.Z( I4-2 ) ) 00260 $ RETURN 00261 B1 = B1*( Z( I4 ) / Z( I4-2 ) ) 00262 B2 = B2 + B1 00263 IF( HUNDRD*MAX( B1, A2 ).LT.B2 ) 00264 $ GO TO 60 00265 50 CONTINUE 00266 60 CONTINUE 00267 B2 = SQRT( CNST3*B2 ) 00268 A2 = DMIN1 / ( ONE+B2**2 ) 00269 GAP2 = HALF*DMIN2 - A2 00270 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN 00271 S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) ) 00272 ELSE 00273 S = MAX( S, A2*( ONE-CNST2*B2 ) ) 00274 TTYPE = -8 00275 END IF 00276 ELSE 00277 * 00278 * Case 9. 00279 * 00280 S = QURTR*DMIN1 00281 IF( DMIN1.EQ.DN1 ) 00282 $ S = HALF*DMIN1 00283 TTYPE = -9 00284 END IF 00285 * 00286 ELSE IF( N0IN.EQ.( N0+2 ) ) THEN 00287 * 00288 * Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN. 00289 * 00290 * Cases 10 and 11. 00291 * 00292 IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN 00293 TTYPE = -10 00294 S = THIRD*DMIN2 00295 IF( Z( NN-5 ).GT.Z( NN-7 ) ) 00296 $ RETURN 00297 B1 = Z( NN-5 ) / Z( NN-7 ) 00298 B2 = B1 00299 IF( B2.EQ.ZERO ) 00300 $ GO TO 80 00301 DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4 00302 IF( Z( I4 ).GT.Z( I4-2 ) ) 00303 $ RETURN 00304 B1 = B1*( Z( I4 ) / Z( I4-2 ) ) 00305 B2 = B2 + B1 00306 IF( HUNDRD*B1.LT.B2 ) 00307 $ GO TO 80 00308 70 CONTINUE 00309 80 CONTINUE 00310 B2 = SQRT( CNST3*B2 ) 00311 A2 = DMIN2 / ( ONE+B2**2 ) 00312 GAP2 = Z( NN-7 ) + Z( NN-9 ) - 00313 $ SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2 00314 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN 00315 S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) ) 00316 ELSE 00317 S = MAX( S, A2*( ONE-CNST2*B2 ) ) 00318 END IF 00319 ELSE 00320 S = QURTR*DMIN2 00321 TTYPE = -11 00322 END IF 00323 ELSE IF( N0IN.GT.( N0+2 ) ) THEN 00324 * 00325 * Case 12, more than two eigenvalues deflated. No information. 00326 * 00327 S = ZERO 00328 TTYPE = -12 00329 END IF 00330 * 00331 TAU = S 00332 RETURN 00333 * 00334 * End of DLASQ4 00335 * 00336 END
1.7.3 
