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