![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLASQ4 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLASQ4 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasq4.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasq4.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasq4.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, 00022 * DN1, DN2, TAU, TTYPE, G ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER I0, N0, N0IN, PP, TTYPE 00026 * REAL DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU 00027 * .. 00028 * .. Array Arguments .. 00029 * REAL Z( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> SLASQ4 computes an approximation TAU to the smallest eigenvalue 00039 *> using values of d from the previous transform. 00040 *> \endverbatim 00041 * 00042 * Arguments: 00043 * ========== 00044 * 00045 *> \param[in] I0 00046 *> \verbatim 00047 *> I0 is INTEGER 00048 *> First index. 00049 *> \endverbatim 00050 *> 00051 *> \param[in] N0 00052 *> \verbatim 00053 *> N0 is INTEGER 00054 *> Last index. 00055 *> \endverbatim 00056 *> 00057 *> \param[in] Z 00058 *> \verbatim 00059 *> Z is REAL array, dimension ( 4*N ) 00060 *> Z holds the qd array. 00061 *> \endverbatim 00062 *> 00063 *> \param[in] PP 00064 *> \verbatim 00065 *> PP is INTEGER 00066 *> PP=0 for ping, PP=1 for pong. 00067 *> \endverbatim 00068 *> 00069 *> \param[in] N0IN 00070 *> \verbatim 00071 *> N0IN is INTEGER 00072 *> The value of N0 at start of EIGTEST. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] DMIN 00076 *> \verbatim 00077 *> DMIN is REAL 00078 *> Minimum value of d. 00079 *> \endverbatim 00080 *> 00081 *> \param[in] DMIN1 00082 *> \verbatim 00083 *> DMIN1 is REAL 00084 *> Minimum value of d, excluding D( N0 ). 00085 *> \endverbatim 00086 *> 00087 *> \param[in] DMIN2 00088 *> \verbatim 00089 *> DMIN2 is REAL 00090 *> Minimum value of d, excluding D( N0 ) and D( N0-1 ). 00091 *> \endverbatim 00092 *> 00093 *> \param[in] DN 00094 *> \verbatim 00095 *> DN is REAL 00096 *> d(N) 00097 *> \endverbatim 00098 *> 00099 *> \param[in] DN1 00100 *> \verbatim 00101 *> DN1 is REAL 00102 *> d(N-1) 00103 *> \endverbatim 00104 *> 00105 *> \param[in] DN2 00106 *> \verbatim 00107 *> DN2 is REAL 00108 *> d(N-2) 00109 *> \endverbatim 00110 *> 00111 *> \param[out] TAU 00112 *> \verbatim 00113 *> TAU is REAL 00114 *> This is the shift. 00115 *> \endverbatim 00116 *> 00117 *> \param[out] TTYPE 00118 *> \verbatim 00119 *> TTYPE is INTEGER 00120 *> Shift type. 00121 *> \endverbatim 00122 *> 00123 *> \param[in,out] G 00124 *> \verbatim 00125 *> G is REAL 00126 *> G is passed as an argument in order to save its value between 00127 *> calls to SLASQ4. 00128 *> \endverbatim 00129 * 00130 * Authors: 00131 * ======== 00132 * 00133 *> \author Univ. of Tennessee 00134 *> \author Univ. of California Berkeley 00135 *> \author Univ. of Colorado Denver 00136 *> \author NAG Ltd. 00137 * 00138 *> \date November 2011 00139 * 00140 *> \ingroup auxOTHERcomputational 00141 * 00142 *> \par Further Details: 00143 * ===================== 00144 *> 00145 *> \verbatim 00146 *> 00147 *> CNST1 = 9/16 00148 *> \endverbatim 00149 *> 00150 * ===================================================================== 00151 SUBROUTINE SLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, 00152 $ DN1, DN2, TAU, TTYPE, G ) 00153 * 00154 * -- LAPACK computational routine (version 3.4.0) -- 00155 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00156 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00157 * November 2011 00158 * 00159 * .. Scalar Arguments .. 00160 INTEGER I0, N0, N0IN, PP, TTYPE 00161 REAL DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU 00162 * .. 00163 * .. Array Arguments .. 00164 REAL Z( * ) 00165 * .. 00166 * 00167 * ===================================================================== 00168 * 00169 * .. Parameters .. 00170 REAL CNST1, CNST2, CNST3 00171 PARAMETER ( CNST1 = 0.5630E0, CNST2 = 1.010E0, 00172 $ CNST3 = 1.050E0 ) 00173 REAL QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD 00174 PARAMETER ( QURTR = 0.250E0, THIRD = 0.3330E0, 00175 $ HALF = 0.50E0, ZERO = 0.0E0, ONE = 1.0E0, 00176 $ TWO = 2.0E0, HUNDRD = 100.0E0 ) 00177 * .. 00178 * .. Local Scalars .. 00179 INTEGER I4, NN, NP 00180 REAL A2, B1, B2, GAM, GAP1, GAP2, S 00181 * .. 00182 * .. Intrinsic Functions .. 00183 INTRINSIC MAX, MIN, SQRT 00184 * .. 00185 * .. Executable Statements .. 00186 * 00187 * A negative DMIN forces the shift to take that absolute value 00188 * TTYPE records the type of shift. 00189 * 00190 IF( DMIN.LE.ZERO ) THEN 00191 TAU = -DMIN 00192 TTYPE = -1 00193 RETURN 00194 END IF 00195 * 00196 NN = 4*N0 + PP 00197 IF( N0IN.EQ.N0 ) THEN 00198 * 00199 * No eigenvalues deflated. 00200 * 00201 IF( DMIN.EQ.DN .OR. DMIN.EQ.DN1 ) THEN 00202 * 00203 B1 = SQRT( Z( NN-3 ) )*SQRT( Z( NN-5 ) ) 00204 B2 = SQRT( Z( NN-7 ) )*SQRT( Z( NN-9 ) ) 00205 A2 = Z( NN-7 ) + Z( NN-5 ) 00206 * 00207 * Cases 2 and 3. 00208 * 00209 IF( DMIN.EQ.DN .AND. DMIN1.EQ.DN1 ) THEN 00210 GAP2 = DMIN2 - A2 - DMIN2*QURTR 00211 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2 ) THEN 00212 GAP1 = A2 - DN - ( B2 / GAP2 )*B2 00213 ELSE 00214 GAP1 = A2 - DN - ( B1+B2 ) 00215 END IF 00216 IF( GAP1.GT.ZERO .AND. GAP1.GT.B1 ) THEN 00217 S = MAX( DN-( B1 / GAP1 )*B1, HALF*DMIN ) 00218 TTYPE = -2 00219 ELSE 00220 S = ZERO 00221 IF( DN.GT.B1 ) 00222 $ S = DN - B1 00223 IF( A2.GT.( B1+B2 ) ) 00224 $ S = MIN( S, A2-( B1+B2 ) ) 00225 S = MAX( S, THIRD*DMIN ) 00226 TTYPE = -3 00227 END IF 00228 ELSE 00229 * 00230 * Case 4. 00231 * 00232 TTYPE = -4 00233 S = QURTR*DMIN 00234 IF( DMIN.EQ.DN ) THEN 00235 GAM = DN 00236 A2 = ZERO 00237 IF( Z( NN-5 ) .GT. Z( NN-7 ) ) 00238 $ RETURN 00239 B2 = Z( NN-5 ) / Z( NN-7 ) 00240 NP = NN - 9 00241 ELSE 00242 NP = NN - 2*PP 00243 B2 = Z( NP-2 ) 00244 GAM = DN1 00245 IF( Z( NP-4 ) .GT. Z( NP-2 ) ) 00246 $ RETURN 00247 A2 = Z( NP-4 ) / Z( NP-2 ) 00248 IF( Z( NN-9 ) .GT. Z( NN-11 ) ) 00249 $ RETURN 00250 B2 = Z( NN-9 ) / Z( NN-11 ) 00251 NP = NN - 13 00252 END IF 00253 * 00254 * Approximate contribution to norm squared from I < NN-1. 00255 * 00256 A2 = A2 + B2 00257 DO 10 I4 = NP, 4*I0 - 1 + PP, -4 00258 IF( B2.EQ.ZERO ) 00259 $ GO TO 20 00260 B1 = B2 00261 IF( Z( I4 ) .GT. Z( I4-2 ) ) 00262 $ RETURN 00263 B2 = B2*( Z( I4 ) / Z( I4-2 ) ) 00264 A2 = A2 + B2 00265 IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) 00266 $ GO TO 20 00267 10 CONTINUE 00268 20 CONTINUE 00269 A2 = CNST3*A2 00270 * 00271 * Rayleigh quotient residual bound. 00272 * 00273 IF( A2.LT.CNST1 ) 00274 $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 ) 00275 END IF 00276 ELSE IF( DMIN.EQ.DN2 ) THEN 00277 * 00278 * Case 5. 00279 * 00280 TTYPE = -5 00281 S = QURTR*DMIN 00282 * 00283 * Compute contribution to norm squared from I > NN-2. 00284 * 00285 NP = NN - 2*PP 00286 B1 = Z( NP-2 ) 00287 B2 = Z( NP-6 ) 00288 GAM = DN2 00289 IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 ) 00290 $ RETURN 00291 A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 ) 00292 * 00293 * Approximate contribution to norm squared from I < NN-2. 00294 * 00295 IF( N0-I0.GT.2 ) THEN 00296 B2 = Z( NN-13 ) / Z( NN-15 ) 00297 A2 = A2 + B2 00298 DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4 00299 IF( B2.EQ.ZERO ) 00300 $ GO TO 40 00301 B1 = B2 00302 IF( Z( I4 ) .GT. Z( I4-2 ) ) 00303 $ RETURN 00304 B2 = B2*( Z( I4 ) / Z( I4-2 ) ) 00305 A2 = A2 + B2 00306 IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) 00307 $ GO TO 40 00308 30 CONTINUE 00309 40 CONTINUE 00310 A2 = CNST3*A2 00311 END IF 00312 * 00313 IF( A2.LT.CNST1 ) 00314 $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 ) 00315 ELSE 00316 * 00317 * Case 6, no information to guide us. 00318 * 00319 IF( TTYPE.EQ.-6 ) THEN 00320 G = G + THIRD*( ONE-G ) 00321 ELSE IF( TTYPE.EQ.-18 ) THEN 00322 G = QURTR*THIRD 00323 ELSE 00324 G = QURTR 00325 END IF 00326 S = G*DMIN 00327 TTYPE = -6 00328 END IF 00329 * 00330 ELSE IF( N0IN.EQ.( N0+1 ) ) THEN 00331 * 00332 * One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN. 00333 * 00334 IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN 00335 * 00336 * Cases 7 and 8. 00337 * 00338 TTYPE = -7 00339 S = THIRD*DMIN1 00340 IF( Z( NN-5 ).GT.Z( NN-7 ) ) 00341 $ RETURN 00342 B1 = Z( NN-5 ) / Z( NN-7 ) 00343 B2 = B1 00344 IF( B2.EQ.ZERO ) 00345 $ GO TO 60 00346 DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4 00347 A2 = B1 00348 IF( Z( I4 ).GT.Z( I4-2 ) ) 00349 $ RETURN 00350 B1 = B1*( Z( I4 ) / Z( I4-2 ) ) 00351 B2 = B2 + B1 00352 IF( HUNDRD*MAX( B1, A2 ).LT.B2 ) 00353 $ GO TO 60 00354 50 CONTINUE 00355 60 CONTINUE 00356 B2 = SQRT( CNST3*B2 ) 00357 A2 = DMIN1 / ( ONE+B2**2 ) 00358 GAP2 = HALF*DMIN2 - A2 00359 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN 00360 S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) ) 00361 ELSE 00362 S = MAX( S, A2*( ONE-CNST2*B2 ) ) 00363 TTYPE = -8 00364 END IF 00365 ELSE 00366 * 00367 * Case 9. 00368 * 00369 S = QURTR*DMIN1 00370 IF( DMIN1.EQ.DN1 ) 00371 $ S = HALF*DMIN1 00372 TTYPE = -9 00373 END IF 00374 * 00375 ELSE IF( N0IN.EQ.( N0+2 ) ) THEN 00376 * 00377 * Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN. 00378 * 00379 * Cases 10 and 11. 00380 * 00381 IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN 00382 TTYPE = -10 00383 S = THIRD*DMIN2 00384 IF( Z( NN-5 ).GT.Z( NN-7 ) ) 00385 $ RETURN 00386 B1 = Z( NN-5 ) / Z( NN-7 ) 00387 B2 = B1 00388 IF( B2.EQ.ZERO ) 00389 $ GO TO 80 00390 DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4 00391 IF( Z( I4 ).GT.Z( I4-2 ) ) 00392 $ RETURN 00393 B1 = B1*( Z( I4 ) / Z( I4-2 ) ) 00394 B2 = B2 + B1 00395 IF( HUNDRD*B1.LT.B2 ) 00396 $ GO TO 80 00397 70 CONTINUE 00398 80 CONTINUE 00399 B2 = SQRT( CNST3*B2 ) 00400 A2 = DMIN2 / ( ONE+B2**2 ) 00401 GAP2 = Z( NN-7 ) + Z( NN-9 ) - 00402 $ SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2 00403 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN 00404 S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) ) 00405 ELSE 00406 S = MAX( S, A2*( ONE-CNST2*B2 ) ) 00407 END IF 00408 ELSE 00409 S = QURTR*DMIN2 00410 TTYPE = -11 00411 END IF 00412 ELSE IF( N0IN.GT.( N0+2 ) ) THEN 00413 * 00414 * Case 12, more than two eigenvalues deflated. No information. 00415 * 00416 S = ZERO 00417 TTYPE = -12 00418 END IF 00419 * 00420 TAU = S 00421 RETURN 00422 * 00423 * End of SLASQ4 00424 * 00425 END