 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, 00002 $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, 00003 $ DN2, G, TAU ) 00004 * 00005 * -- LAPACK routine (version 3.2.2) -- 00006 * 00007 * -- Contributed by Osni Marques of the Lawrence Berkeley National -- 00008 * -- Laboratory and Beresford Parlett of the Univ. of California at -- 00009 * -- Berkeley -- 00010 * -- June 2010 -- 00011 * 00012 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00013 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00014 * 00015 * .. Scalar Arguments .. 00016 LOGICAL IEEE 00017 INTEGER I0, ITER, N0, NDIV, NFAIL, PP 00018 REAL DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, 00019 $ QMAX, SIGMA, TAU 00020 * .. 00021 * .. Array Arguments .. 00022 REAL Z( * ) 00023 * .. 00024 * 00025 * Purpose 00026 * ======= 00027 * 00028 * SLASQ3 checks for deflation, computes a shift (TAU) and calls dqds. 00029 * In case of failure it changes shifts, and tries again until output 00030 * is positive. 00031 * 00032 * Arguments 00033 * ========= 00034 * 00035 * I0 (input) INTEGER 00036 * First index. 00037 * 00038 * N0 (input/output) INTEGER 00039 * Last index. 00040 * 00041 * Z (input) REAL array, dimension ( 4*N ) 00042 * Z holds the qd array. 00043 * 00044 * PP (input/output) INTEGER 00045 * PP=0 for ping, PP=1 for pong. 00046 * PP=2 indicates that flipping was applied to the Z array 00047 * and that the initial tests for deflation should not be 00048 * performed. 00049 * 00050 * DMIN (output) REAL 00051 * Minimum value of d. 00052 * 00053 * SIGMA (output) REAL 00054 * Sum of shifts used in current segment. 00055 * 00056 * DESIG (input/output) REAL 00057 * Lower order part of SIGMA 00058 * 00059 * QMAX (input) REAL 00060 * Maximum value of q. 00061 * 00062 * NFAIL (output) INTEGER 00063 * Number of times shift was too big. 00064 * 00065 * ITER (output) INTEGER 00066 * Number of iterations. 00067 * 00068 * NDIV (output) INTEGER 00069 * Number of divisions. 00070 * 00071 * IEEE (input) LOGICAL 00072 * Flag for IEEE or non IEEE arithmetic (passed to SLASQ5). 00073 * 00074 * TTYPE (input/output) INTEGER 00075 * Shift type. 00076 * 00077 * DMIN1 (input/output) REAL 00078 * 00079 * DMIN2 (input/output) REAL 00080 * 00081 * DN (input/output) REAL 00082 * 00083 * DN1 (input/output) REAL 00084 * 00085 * DN2 (input/output) REAL 00086 * 00087 * G (input/output) REAL 00088 * 00089 * TAU (input/output) REAL 00090 * 00091 * These are passed as arguments in order to save their values 00092 * between calls to SLASQ3. 00093 * 00094 * ===================================================================== 00095 * 00096 * .. Parameters .. 00097 REAL CBIAS 00098 PARAMETER ( CBIAS = 1.50E0 ) 00099 REAL ZERO, QURTR, HALF, ONE, TWO, HUNDRD 00100 PARAMETER ( ZERO = 0.0E0, QURTR = 0.250E0, HALF = 0.5E0, 00101 $ ONE = 1.0E0, TWO = 2.0E0, HUNDRD = 100.0E0 ) 00102 * .. 00103 * .. Local Scalars .. 00104 INTEGER IPN4, J4, N0IN, NN, TTYPE 00105 REAL EPS, S, T, TEMP, TOL, TOL2 00106 * .. 00107 * .. External Subroutines .. 00108 EXTERNAL SLASQ4, SLASQ5, SLASQ6 00109 * .. 00110 * .. External Function .. 00111 REAL SLAMCH 00112 LOGICAL SISNAN 00113 EXTERNAL SISNAN, SLAMCH 00114 * .. 00115 * .. Intrinsic Functions .. 00116 INTRINSIC ABS, MAX, MIN, SQRT 00117 * .. 00118 * .. Executable Statements .. 00119 * 00120 N0IN = N0 00121 EPS = SLAMCH( 'Precision' ) 00122 TOL = EPS*HUNDRD 00123 TOL2 = TOL**2 00124 * 00125 * Check for deflation. 00126 * 00127 10 CONTINUE 00128 * 00129 IF( N0.LT.I0 ) 00130 $ RETURN 00131 IF( N0.EQ.I0 ) 00132 $ GO TO 20 00133 NN = 4*N0 + PP 00134 IF( N0.EQ.( I0+1 ) ) 00135 $ GO TO 40 00136 * 00137 * Check whether E(N0-1) is negligible, 1 eigenvalue. 00138 * 00139 IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND. 00140 $ Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) ) 00141 $ GO TO 30 00142 * 00143 20 CONTINUE 00144 * 00145 Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA 00146 N0 = N0 - 1 00147 GO TO 10 00148 * 00149 * Check whether E(N0-2) is negligible, 2 eigenvalues. 00150 * 00151 30 CONTINUE 00152 * 00153 IF( Z( NN-9 ).GT.TOL2*SIGMA .AND. 00154 $ Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) ) 00155 $ GO TO 50 00156 * 00157 40 CONTINUE 00158 * 00159 IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN 00160 S = Z( NN-3 ) 00161 Z( NN-3 ) = Z( NN-7 ) 00162 Z( NN-7 ) = S 00163 END IF 00164 IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2 ) THEN 00165 T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) ) 00166 S = Z( NN-3 )*( Z( NN-5 ) / T ) 00167 IF( S.LE.T ) THEN 00168 S = Z( NN-3 )*( Z( NN-5 ) / 00169 $ ( T*( ONE+SQRT( ONE+S / T ) ) ) ) 00170 ELSE 00171 S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) ) 00172 END IF 00173 T = Z( NN-7 ) + ( S+Z( NN-5 ) ) 00174 Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T ) 00175 Z( NN-7 ) = T 00176 END IF 00177 Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA 00178 Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA 00179 N0 = N0 - 2 00180 GO TO 10 00181 * 00182 50 CONTINUE 00183 IF( PP.EQ.2 ) 00184 $ PP = 0 00185 * 00186 * Reverse the qd-array, if warranted. 00187 * 00188 IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN 00189 IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN 00190 IPN4 = 4*( I0+N0 ) 00191 DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4 00192 TEMP = Z( J4-3 ) 00193 Z( J4-3 ) = Z( IPN4-J4-3 ) 00194 Z( IPN4-J4-3 ) = TEMP 00195 TEMP = Z( J4-2 ) 00196 Z( J4-2 ) = Z( IPN4-J4-2 ) 00197 Z( IPN4-J4-2 ) = TEMP 00198 TEMP = Z( J4-1 ) 00199 Z( J4-1 ) = Z( IPN4-J4-5 ) 00200 Z( IPN4-J4-5 ) = TEMP 00201 TEMP = Z( J4 ) 00202 Z( J4 ) = Z( IPN4-J4-4 ) 00203 Z( IPN4-J4-4 ) = TEMP 00204 60 CONTINUE 00205 IF( N0-I0.LE.4 ) THEN 00206 Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 ) 00207 Z( 4*N0-PP ) = Z( 4*I0-PP ) 00208 END IF 00209 DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) ) 00210 Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ), 00211 $ Z( 4*I0+PP+3 ) ) 00212 Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ), 00213 $ Z( 4*I0-PP+4 ) ) 00214 QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) ) 00215 DMIN = -ZERO 00216 END IF 00217 END IF 00218 * 00219 * Choose a shift. 00220 * 00221 CALL SLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, 00222 $ DN2, TAU, TTYPE, G ) 00223 * 00224 * Call dqds until DMIN > 0. 00225 * 00226 70 CONTINUE 00227 * 00228 CALL SLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN, 00229 $ DN1, DN2, IEEE ) 00230 * 00231 NDIV = NDIV + ( N0-I0+2 ) 00232 ITER = ITER + 1 00233 * 00234 * Check status. 00235 * 00236 IF( DMIN.GE.ZERO .AND. DMIN1.GT.ZERO ) THEN 00237 * 00238 * Success. 00239 * 00240 GO TO 90 00241 * 00242 ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND. 00243 $ Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND. 00244 $ ABS( DN ).LT.TOL*SIGMA ) THEN 00245 * 00246 * Convergence hidden by negative DN. 00247 * 00248 Z( 4*( N0-1 )-PP+2 ) = ZERO 00249 DMIN = ZERO 00250 GO TO 90 00251 ELSE IF( DMIN.LT.ZERO ) THEN 00252 * 00253 * TAU too big. Select new TAU and try again. 00254 * 00255 NFAIL = NFAIL + 1 00256 IF( TTYPE.LT.-22 ) THEN 00257 * 00258 * Failed twice. Play it safe. 00259 * 00260 TAU = ZERO 00261 ELSE IF( DMIN1.GT.ZERO ) THEN 00262 * 00263 * Late failure. Gives excellent shift. 00264 * 00265 TAU = ( TAU+DMIN )*( ONE-TWO*EPS ) 00266 TTYPE = TTYPE - 11 00267 ELSE 00268 * 00269 * Early failure. Divide by 4. 00270 * 00271 TAU = QURTR*TAU 00272 TTYPE = TTYPE - 12 00273 END IF 00274 GO TO 70 00275 ELSE IF( SISNAN( DMIN ) ) THEN 00276 * 00277 * NaN. 00278 * 00279 IF( TAU.EQ.ZERO ) THEN 00280 GO TO 80 00281 ELSE 00282 TAU = ZERO 00283 GO TO 70 00284 END IF 00285 ELSE 00286 * 00287 * Possible underflow. Play it safe. 00288 * 00289 GO TO 80 00290 END IF 00291 * 00292 * Risk of underflow. 00293 * 00294 80 CONTINUE 00295 CALL SLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 ) 00296 NDIV = NDIV + ( N0-I0+2 ) 00297 ITER = ITER + 1 00298 TAU = ZERO 00299 * 00300 90 CONTINUE 00301 IF( TAU.LT.SIGMA ) THEN 00302 DESIG = DESIG + TAU 00303 T = SIGMA + DESIG 00304 DESIG = DESIG - ( T-SIGMA ) 00305 ELSE 00306 T = SIGMA + TAU 00307 DESIG = SIGMA - ( T-TAU ) + DESIG 00308 END IF 00309 SIGMA = T 00310 * 00311 RETURN 00312 * 00313 * End of SLASQ3 00314 * 00315 END
1.7.3 
