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