![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLASQ3 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLASQ3 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq3.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq3.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq3.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, 00022 * ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, 00023 * DN2, G, TAU ) 00024 * 00025 * .. Scalar Arguments .. 00026 * LOGICAL IEEE 00027 * INTEGER I0, ITER, N0, NDIV, NFAIL, PP 00028 * DOUBLE PRECISION DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, 00029 * $ QMAX, SIGMA, TAU 00030 * .. 00031 * .. Array Arguments .. 00032 * DOUBLE PRECISION Z( * ) 00033 * .. 00034 * 00035 * 00036 *> \par Purpose: 00037 * ============= 00038 *> 00039 *> \verbatim 00040 *> 00041 *> DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds. 00042 *> In case of failure it changes shifts, and tries again until output 00043 *> is positive. 00044 *> \endverbatim 00045 * 00046 * Arguments: 00047 * ========== 00048 * 00049 *> \param[in] I0 00050 *> \verbatim 00051 *> I0 is INTEGER 00052 *> First index. 00053 *> \endverbatim 00054 *> 00055 *> \param[in,out] N0 00056 *> \verbatim 00057 *> N0 is INTEGER 00058 *> Last index. 00059 *> \endverbatim 00060 *> 00061 *> \param[in] Z 00062 *> \verbatim 00063 *> Z is DOUBLE PRECISION array, dimension ( 4*N ) 00064 *> Z holds the qd array. 00065 *> \endverbatim 00066 *> 00067 *> \param[in,out] PP 00068 *> \verbatim 00069 *> PP is INTEGER 00070 *> PP=0 for ping, PP=1 for pong. 00071 *> PP=2 indicates that flipping was applied to the Z array 00072 *> and that the initial tests for deflation should not be 00073 *> performed. 00074 *> \endverbatim 00075 *> 00076 *> \param[out] DMIN 00077 *> \verbatim 00078 *> DMIN is DOUBLE PRECISION 00079 *> Minimum value of d. 00080 *> \endverbatim 00081 *> 00082 *> \param[out] SIGMA 00083 *> \verbatim 00084 *> SIGMA is DOUBLE PRECISION 00085 *> Sum of shifts used in current segment. 00086 *> \endverbatim 00087 *> 00088 *> \param[in,out] DESIG 00089 *> \verbatim 00090 *> DESIG is DOUBLE PRECISION 00091 *> Lower order part of SIGMA 00092 *> \endverbatim 00093 *> 00094 *> \param[in] QMAX 00095 *> \verbatim 00096 *> QMAX is DOUBLE PRECISION 00097 *> Maximum value of q. 00098 *> \endverbatim 00099 *> 00100 *> \param[out] NFAIL 00101 *> \verbatim 00102 *> NFAIL is INTEGER 00103 *> Number of times shift was too big. 00104 *> \endverbatim 00105 *> 00106 *> \param[out] ITER 00107 *> \verbatim 00108 *> ITER is INTEGER 00109 *> Number of iterations. 00110 *> \endverbatim 00111 *> 00112 *> \param[out] NDIV 00113 *> \verbatim 00114 *> NDIV is INTEGER 00115 *> Number of divisions. 00116 *> \endverbatim 00117 *> 00118 *> \param[in] IEEE 00119 *> \verbatim 00120 *> IEEE is LOGICAL 00121 *> Flag for IEEE or non IEEE arithmetic (passed to DLASQ5). 00122 *> \endverbatim 00123 *> 00124 *> \param[in,out] TTYPE 00125 *> \verbatim 00126 *> TTYPE is INTEGER 00127 *> Shift type. 00128 *> \endverbatim 00129 *> 00130 *> \param[in,out] DMIN1 00131 *> \verbatim 00132 *> DMIN1 is DOUBLE PRECISION 00133 *> \endverbatim 00134 *> 00135 *> \param[in,out] DMIN2 00136 *> \verbatim 00137 *> DMIN2 is DOUBLE PRECISION 00138 *> \endverbatim 00139 *> 00140 *> \param[in,out] DN 00141 *> \verbatim 00142 *> DN is DOUBLE PRECISION 00143 *> \endverbatim 00144 *> 00145 *> \param[in,out] DN1 00146 *> \verbatim 00147 *> DN1 is DOUBLE PRECISION 00148 *> \endverbatim 00149 *> 00150 *> \param[in,out] DN2 00151 *> \verbatim 00152 *> DN2 is DOUBLE PRECISION 00153 *> \endverbatim 00154 *> 00155 *> \param[in,out] G 00156 *> \verbatim 00157 *> G is DOUBLE PRECISION 00158 *> \endverbatim 00159 *> 00160 *> \param[in,out] TAU 00161 *> \verbatim 00162 *> TAU is DOUBLE PRECISION 00163 *> 00164 *> These are passed as arguments in order to save their values 00165 *> between calls to DLASQ3. 00166 *> \endverbatim 00167 * 00168 * Authors: 00169 * ======== 00170 * 00171 *> \author Univ. of Tennessee 00172 *> \author Univ. of California Berkeley 00173 *> \author Univ. of Colorado Denver 00174 *> \author NAG Ltd. 00175 * 00176 *> \date April 2012 00177 * 00178 *> \ingroup auxOTHERcomputational 00179 * 00180 * ===================================================================== 00181 SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, 00182 $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, 00183 $ DN2, G, TAU ) 00184 * 00185 * -- LAPACK computational routine (version 3.4.1) -- 00186 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00187 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00188 * April 2012 00189 * 00190 * .. Scalar Arguments .. 00191 LOGICAL IEEE 00192 INTEGER I0, ITER, N0, NDIV, NFAIL, PP 00193 DOUBLE PRECISION DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, 00194 $ QMAX, SIGMA, TAU 00195 * .. 00196 * .. Array Arguments .. 00197 DOUBLE PRECISION Z( * ) 00198 * .. 00199 * 00200 * ===================================================================== 00201 * 00202 * .. Parameters .. 00203 DOUBLE PRECISION CBIAS 00204 PARAMETER ( CBIAS = 1.50D0 ) 00205 DOUBLE PRECISION ZERO, QURTR, HALF, ONE, TWO, HUNDRD 00206 PARAMETER ( ZERO = 0.0D0, QURTR = 0.250D0, HALF = 0.5D0, 00207 $ ONE = 1.0D0, TWO = 2.0D0, HUNDRD = 100.0D0 ) 00208 * .. 00209 * .. Local Scalars .. 00210 INTEGER IPN4, J4, N0IN, NN, TTYPE 00211 DOUBLE PRECISION EPS, S, T, TEMP, TOL, TOL2 00212 * .. 00213 * .. External Subroutines .. 00214 EXTERNAL DLASQ4, DLASQ5, DLASQ6 00215 * .. 00216 * .. External Function .. 00217 DOUBLE PRECISION DLAMCH 00218 LOGICAL DISNAN 00219 EXTERNAL DISNAN, DLAMCH 00220 * .. 00221 * .. Intrinsic Functions .. 00222 INTRINSIC ABS, MAX, MIN, SQRT 00223 * .. 00224 * .. Executable Statements .. 00225 * 00226 N0IN = N0 00227 EPS = DLAMCH( 'Precision' ) 00228 TOL = EPS*HUNDRD 00229 TOL2 = TOL**2 00230 * 00231 * Check for deflation. 00232 * 00233 10 CONTINUE 00234 * 00235 IF( N0.LT.I0 ) 00236 $ RETURN 00237 IF( N0.EQ.I0 ) 00238 $ GO TO 20 00239 NN = 4*N0 + PP 00240 IF( N0.EQ.( I0+1 ) ) 00241 $ GO TO 40 00242 * 00243 * Check whether E(N0-1) is negligible, 1 eigenvalue. 00244 * 00245 IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND. 00246 $ Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) ) 00247 $ GO TO 30 00248 * 00249 20 CONTINUE 00250 * 00251 Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA 00252 N0 = N0 - 1 00253 GO TO 10 00254 * 00255 * Check whether E(N0-2) is negligible, 2 eigenvalues. 00256 * 00257 30 CONTINUE 00258 * 00259 IF( Z( NN-9 ).GT.TOL2*SIGMA .AND. 00260 $ Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) ) 00261 $ GO TO 50 00262 * 00263 40 CONTINUE 00264 * 00265 IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN 00266 S = Z( NN-3 ) 00267 Z( NN-3 ) = Z( NN-7 ) 00268 Z( NN-7 ) = S 00269 END IF 00270 IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2 ) THEN 00271 T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) ) 00272 S = Z( NN-3 )*( Z( NN-5 ) / T ) 00273 IF( S.LE.T ) THEN 00274 S = Z( NN-3 )*( Z( NN-5 ) / 00275 $ ( T*( ONE+SQRT( ONE+S / T ) ) ) ) 00276 ELSE 00277 S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) ) 00278 END IF 00279 T = Z( NN-7 ) + ( S+Z( NN-5 ) ) 00280 Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T ) 00281 Z( NN-7 ) = T 00282 END IF 00283 Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA 00284 Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA 00285 N0 = N0 - 2 00286 GO TO 10 00287 * 00288 50 CONTINUE 00289 IF( PP.EQ.2 ) 00290 $ PP = 0 00291 * 00292 * Reverse the qd-array, if warranted. 00293 * 00294 IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN 00295 IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN 00296 IPN4 = 4*( I0+N0 ) 00297 DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4 00298 TEMP = Z( J4-3 ) 00299 Z( J4-3 ) = Z( IPN4-J4-3 ) 00300 Z( IPN4-J4-3 ) = TEMP 00301 TEMP = Z( J4-2 ) 00302 Z( J4-2 ) = Z( IPN4-J4-2 ) 00303 Z( IPN4-J4-2 ) = TEMP 00304 TEMP = Z( J4-1 ) 00305 Z( J4-1 ) = Z( IPN4-J4-5 ) 00306 Z( IPN4-J4-5 ) = TEMP 00307 TEMP = Z( J4 ) 00308 Z( J4 ) = Z( IPN4-J4-4 ) 00309 Z( IPN4-J4-4 ) = TEMP 00310 60 CONTINUE 00311 IF( N0-I0.LE.4 ) THEN 00312 Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 ) 00313 Z( 4*N0-PP ) = Z( 4*I0-PP ) 00314 END IF 00315 DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) ) 00316 Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ), 00317 $ Z( 4*I0+PP+3 ) ) 00318 Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ), 00319 $ Z( 4*I0-PP+4 ) ) 00320 QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) ) 00321 DMIN = -ZERO 00322 END IF 00323 END IF 00324 * 00325 * Choose a shift. 00326 * 00327 CALL DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, 00328 $ DN2, TAU, TTYPE, G ) 00329 * 00330 * Call dqds until DMIN > 0. 00331 * 00332 70 CONTINUE 00333 * 00334 CALL DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, 00335 $ DN1, DN2, IEEE, EPS ) 00336 * 00337 NDIV = NDIV + ( N0-I0+2 ) 00338 ITER = ITER + 1 00339 * 00340 * Check status. 00341 * 00342 IF( DMIN.GE.ZERO .AND. DMIN1.GE.ZERO ) THEN 00343 * 00344 * Success. 00345 * 00346 GO TO 90 00347 * 00348 ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND. 00349 $ Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND. 00350 $ ABS( DN ).LT.TOL*SIGMA ) THEN 00351 * 00352 * Convergence hidden by negative DN. 00353 * 00354 Z( 4*( N0-1 )-PP+2 ) = ZERO 00355 DMIN = ZERO 00356 GO TO 90 00357 ELSE IF( DMIN.LT.ZERO ) THEN 00358 * 00359 * TAU too big. Select new TAU and try again. 00360 * 00361 NFAIL = NFAIL + 1 00362 IF( TTYPE.LT.-22 ) THEN 00363 * 00364 * Failed twice. Play it safe. 00365 * 00366 TAU = ZERO 00367 ELSE IF( DMIN1.GT.ZERO ) THEN 00368 * 00369 * Late failure. Gives excellent shift. 00370 * 00371 TAU = ( TAU+DMIN )*( ONE-TWO*EPS ) 00372 TTYPE = TTYPE - 11 00373 ELSE 00374 * 00375 * Early failure. Divide by 4. 00376 * 00377 TAU = QURTR*TAU 00378 TTYPE = TTYPE - 12 00379 END IF 00380 GO TO 70 00381 ELSE IF( DISNAN( DMIN ) ) THEN 00382 * 00383 * NaN. 00384 * 00385 IF( TAU.EQ.ZERO ) THEN 00386 GO TO 80 00387 ELSE 00388 TAU = ZERO 00389 GO TO 70 00390 END IF 00391 ELSE 00392 * 00393 * Possible underflow. Play it safe. 00394 * 00395 GO TO 80 00396 END IF 00397 * 00398 * Risk of underflow. 00399 * 00400 80 CONTINUE 00401 CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 ) 00402 NDIV = NDIV + ( N0-I0+2 ) 00403 ITER = ITER + 1 00404 TAU = ZERO 00405 * 00406 90 CONTINUE 00407 IF( TAU.LT.SIGMA ) THEN 00408 DESIG = DESIG + TAU 00409 T = SIGMA + DESIG 00410 DESIG = DESIG - ( T-SIGMA ) 00411 ELSE 00412 T = SIGMA + TAU 00413 DESIG = SIGMA - ( T-TAU ) + DESIG 00414 END IF 00415 SIGMA = T 00416 * 00417 RETURN 00418 * 00419 * End of DLASQ3 00420 * 00421 END