1 SUBROUTINE SLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, 2 $ ITER, NDIV, IEEE ) 3* 4* -- LAPACK auxiliary routine (instrumented to count ops, version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* May 17, 2000 8* 9* .. Scalar Arguments .. 10 LOGICAL IEEE 11 INTEGER I0, ITER, N0, NDIV, NFAIL, PP 12 REAL DESIG, DMIN, QMAX, SIGMA 13* .. 14* .. Array Arguments .. 15 REAL Z( * ) 16* .. 17* .. Common block to return operation count .. 18 COMMON / LATIME / OPS, ITCNT 19* .. 20* .. Scalars in Common .. 21 REAL ITCNT, OPS 22* .. 23* 24* Purpose 25* ======= 26* 27* SLASQ3 checks for deflation, computes a shift (TAU) and calls dqds. 28* In case of failure it changes shifts, and tries again until output 29* is positive. 30* 31* Arguments 32* ========= 33* 34* I0 (input) INTEGER 35* First index. 36* 37* N0 (input) INTEGER 38* Last index. 39* 40* Z (input) REAL array, dimension ( 4*N ) 41* Z holds the qd array. 42* 43* PP (input) INTEGER 44* PP=0 for ping, PP=1 for pong. 45* 46* DMIN (output) REAL 47* Minimum value of d. 48* 49* SIGMA (output) REAL 50* Sum of shifts used in current segment. 51* 52* DESIG (input/output) REAL 53* Lower order part of SIGMA 54* 55* QMAX (input) REAL 56* Maximum value of q. 57* 58* NFAIL (output) INTEGER 59* Number of times shift was too big. 60* 61* ITER (output) INTEGER 62* Number of iterations. 63* 64* NDIV (output) INTEGER 65* Number of divisions. 66* 67* TTYPE (output) INTEGER 68* Shift type. 69* 70* IEEE (input) LOGICAL 71* Flag for IEEE or non IEEE arithmetic (passed to SLASQ5). 72* 73* ===================================================================== 74* 75* .. Parameters .. 76 REAL CBIAS 77 PARAMETER ( CBIAS = 1.50E0 ) 78 REAL ZERO, QURTR, HALF, ONE, TWO, HUNDRD 79 PARAMETER ( ZERO = 0.0E0, QURTR = 0.250E0, HALF = 0.5E0, 80 $ ONE = 1.0E0, TWO = 2.0E0, HUNDRD = 100.0E0 ) 81* .. 82* .. Local Scalars .. 83 INTEGER IPN4, J4, N0IN, NN, TTYPE 84 REAL DMIN1, DMIN2, DN, DN1, DN2, EPS, S, SAFMIN, T, 85 $ TAU, TEMP, TOL, TOL2 86* .. 87* .. External Subroutines .. 88 EXTERNAL SLASQ4, SLASQ5, SLASQ6 89* .. 90* .. External Function .. 91 REAL SLAMCH 92 EXTERNAL SLAMCH 93* .. 94* .. Intrinsic Functions .. 95 INTRINSIC ABS, MIN, REAL, SQRT 96* .. 97* .. Save statement .. 98 SAVE TTYPE 99 SAVE DMIN1, DMIN2, DN, DN1, DN2, TAU 100* .. 101* .. Data statement .. 102 DATA TTYPE / 0 / 103 DATA DMIN1 / ZERO /, DMIN2 / ZERO /, DN / ZERO /, 104 $ DN1 / ZERO /, DN2 / ZERO /, TAU / ZERO / 105* .. 106* .. Executable Statements .. 107* 108 OPS = OPS + REAL( 2 ) 109 N0IN = N0 110 EPS = SLAMCH( 'Precision' ) 111 SAFMIN = SLAMCH( 'Safe minimum' ) 112 TOL = EPS*HUNDRD 113 TOL2 = TOL**2 114* 115* Check for deflation. 116* 117 10 CONTINUE 118* 119 IF( N0.LT.I0 ) 120 $ RETURN 121 IF( N0.EQ.I0 ) 122 $ GO TO 20 123 NN = 4*N0 + PP 124 IF( N0.EQ.( I0+1 ) ) 125 $ GO TO 40 126* 127* Check whether E(N0-1) is negligible, 1 eigenvalue. 128* 129 OPS = OPS + REAL( 3 ) 130 IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND. 131 $ Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) ) 132 $ GO TO 30 133* 134 20 CONTINUE 135* 136 OPS = OPS + REAL( 1 ) 137 Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA 138 N0 = N0 - 1 139 GO TO 10 140* 141* Check whether E(N0-2) is negligible, 2 eigenvalues. 142* 143 30 CONTINUE 144* 145 OPS = OPS + REAL( 2 ) 146 IF( Z( NN-9 ).GT.TOL2*SIGMA .AND. 147 $ Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) ) 148 $ GO TO 50 149* 150 40 CONTINUE 151* 152 IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN 153 S = Z( NN-3 ) 154 Z( NN-3 ) = Z( NN-7 ) 155 Z( NN-7 ) = S 156 END IF 157 OPS = OPS + REAL( 3 ) 158 IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2 ) THEN 159 OPS = OPS + REAL( 5 ) 160 T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) ) 161 S = Z( NN-3 )*( Z( NN-5 ) / T ) 162 IF( S.LE.T ) THEN 163 OPS = OPS + REAL( 7 ) 164 S = Z( NN-3 )*( Z( NN-5 ) / 165 $ ( T*( ONE+SQRT( ONE+S / T ) ) ) ) 166 ELSE 167 OPS = OPS + REAL( 6 ) 168 S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) ) 169 END IF 170 OPS = OPS + REAL( 4 ) 171 T = Z( NN-7 ) + ( S+Z( NN-5 ) ) 172 Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T ) 173 Z( NN-7 ) = T 174 END IF 175 Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA 176 Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA 177 N0 = N0 - 2 178 GO TO 10 179* 180 50 CONTINUE 181* 182* Reverse the qd-array, if warranted. 183* 184 IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN 185 OPS = OPS + REAL( 1 ) 186 IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN 187 IPN4 = 4*( I0+N0 ) 188 DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4 189 TEMP = Z( J4-3 ) 190 Z( J4-3 ) = Z( IPN4-J4-3 ) 191 Z( IPN4-J4-3 ) = TEMP 192 TEMP = Z( J4-2 ) 193 Z( J4-2 ) = Z( IPN4-J4-2 ) 194 Z( IPN4-J4-2 ) = TEMP 195 TEMP = Z( J4-1 ) 196 Z( J4-1 ) = Z( IPN4-J4-5 ) 197 Z( IPN4-J4-5 ) = TEMP 198 TEMP = Z( J4 ) 199 Z( J4 ) = Z( IPN4-J4-4 ) 200 Z( IPN4-J4-4 ) = TEMP 201 60 CONTINUE 202 IF( N0-I0.LE.4 ) THEN 203 Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 ) 204 Z( 4*N0-PP ) = Z( 4*I0-PP ) 205 END IF 206 DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) ) 207 Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ), 208 $ Z( 4*I0+PP+3 ) ) 209 Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ), 210 $ Z( 4*I0-PP+4 ) ) 211 QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) ) 212 DMIN = -ZERO 213 END IF 214 END IF 215* 216 70 CONTINUE 217* 218 IF( DMIN.LT.ZERO .OR. SAFMIN*QMAX.LT.MIN( Z( 4*N0+PP-1 ), 219 $ Z( 4*N0+PP-9 ), DMIN2+Z( 4*N0-PP ) ) ) THEN 220* 221* Choose a shift. 222* 223 CALL SLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, 224 $ DN2, TAU, TTYPE ) 225* 226* Call dqds until DMIN > 0. 227* 228 80 CONTINUE 229* 230 CALL SLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN, 231 $ DN1, DN2, IEEE ) 232* 233 NDIV = NDIV + ( N0-I0+2 ) 234 ITER = ITER + 1 235* 236* Check status. 237* 238 IF( DMIN.GE.ZERO .AND. DMIN1.GT.ZERO ) THEN 239* 240* Success. 241* 242 GO TO 100 243* 244 ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND. 245 $ Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND. 246 $ ABS( DN ).LT.TOL*SIGMA ) THEN 247* 248* Convergence hidden by negative DN. 249* 250 OPS = OPS + REAL( 2 ) 251 Z( 4*( N0-1 )-PP+2 ) = ZERO 252 DMIN = ZERO 253 GO TO 100 254 ELSE IF( DMIN.LT.ZERO ) THEN 255* 256* TAU too big. Select new TAU and try again. 257* 258 NFAIL = NFAIL + 1 259 IF( TTYPE.LT.-22 ) THEN 260* 261* Failed twice. Play it safe. 262* 263 TAU = ZERO 264 ELSE IF( DMIN1.GT.ZERO ) THEN 265* 266* Late failure. Gives excellent shift. 267* 268 OPS = OPS + REAL( 4 ) 269 TAU = ( TAU+DMIN )*( ONE-TWO*EPS ) 270 TTYPE = TTYPE - 11 271 ELSE 272* 273* Early failure. Divide by 4. 274* 275 OPS = OPS + REAL( 1 ) 276 TAU = QURTR*TAU 277 TTYPE = TTYPE - 12 278 END IF 279 GO TO 80 280 ELSE IF( DMIN.NE.DMIN ) THEN 281* 282* NaN. 283* 284 TAU = ZERO 285 GO TO 80 286 ELSE 287* 288* Possible underflow. Play it safe. 289* 290 GO TO 90 291 END IF 292 END IF 293* 294* Risk of underflow. 295* 296 90 CONTINUE 297 CALL SLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 ) 298 NDIV = NDIV + ( N0-I0+2 ) 299 ITER = ITER + 1 300 TAU = ZERO 301* 302 100 CONTINUE 303 OPS = OPS + REAL( 4 ) 304 IF( TAU.LT.SIGMA ) THEN 305 DESIG = DESIG + TAU 306 T = SIGMA + DESIG 307 DESIG = DESIG - ( T-SIGMA ) 308 ELSE 309 T = SIGMA + TAU 310 DESIG = SIGMA - ( T-TAU ) + DESIG 311 END IF 312 SIGMA = T 313* 314 RETURN 315* 316* End of SLASQ3 317* 318 END 319