1 SUBROUTINE SLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO ) 2* 3* -- LAPACK routine (instrumented to count operations, version 3.0) -- 4* Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, 5* Courant Institute, NAG Ltd., and Rice University 6* June 30, 1999 7* 8* .. Scalar Arguments .. 9 LOGICAL ORGATI 10 INTEGER INFO, KNITER 11 REAL FINIT, RHO, TAU 12* .. 13* .. Array Arguments .. 14 REAL D( 3 ), Z( 3 ) 15* .. 16* Common block to return operation count and iteration count 17* ITCNT is unchanged, OPS is only incremented 18* .. Common blocks .. 19 COMMON / LATIME / OPS, ITCNT 20* .. 21* .. Scalars in Common .. 22 REAL ITCNT, OPS 23* .. 24* 25* Purpose 26* ======= 27* 28* SLAED6 computes the positive or negative root (closest to the origin) 29* of 30* z(1) z(2) z(3) 31* f(x) = rho + --------- + ---------- + --------- 32* d(1)-x d(2)-x d(3)-x 33* 34* It is assumed that 35* 36* if ORGATI = .true. the root is between d(2) and d(3); 37* otherwise it is between d(1) and d(2) 38* 39* This routine will be called by SLAED4 when necessary. In most cases, 40* the root sought is the smallest in magnitude, though it might not be 41* in some extremely rare situations. 42* 43* Arguments 44* ========= 45* 46* KNITER (input) INTEGER 47* Refer to SLAED4 for its significance. 48* 49* ORGATI (input) LOGICAL 50* If ORGATI is true, the needed root is between d(2) and 51* d(3); otherwise it is between d(1) and d(2). See 52* SLAED4 for further details. 53* 54* RHO (input) REAL 55* Refer to the equation f(x) above. 56* 57* D (input) REAL array, dimension (3) 58* D satisfies d(1) < d(2) < d(3). 59* 60* Z (input) REAL array, dimension (3) 61* Each of the elements in z must be positive. 62* 63* FINIT (input) REAL 64* The value of f at 0. It is more accurate than the one 65* evaluated inside this routine (if someone wants to do 66* so). 67* 68* TAU (output) REAL 69* The root of the equation f(x). 70* 71* INFO (output) INTEGER 72* = 0: successful exit 73* > 0: if INFO = 1, failure to converge 74* 75* Further Details 76* =============== 77* 78* Based on contributions by 79* Ren-Cang Li, Computer Science Division, University of California 80* at Berkeley, USA 81* 82* ===================================================================== 83* 84* .. Parameters .. 85 INTEGER MAXIT 86 PARAMETER ( MAXIT = 20 ) 87 REAL ZERO, ONE, TWO, THREE, FOUR, EIGHT 88 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, 89 $ THREE = 3.0E0, FOUR = 4.0E0, EIGHT = 8.0E0 ) 90* .. 91* .. External Functions .. 92 REAL SLAMCH 93 EXTERNAL SLAMCH 94* .. 95* .. Local Arrays .. 96 REAL DSCALE( 3 ), ZSCALE( 3 ) 97* .. 98* .. Local Scalars .. 99 LOGICAL FIRST, SCALE 100 INTEGER I, ITER, NITER 101 REAL A, B, BASE, C, DDF, DF, EPS, ERRETM, ETA, F, 102 $ FC, SCLFAC, SCLINV, SMALL1, SMALL2, SMINV1, 103 $ SMINV2, TEMP, TEMP1, TEMP2, TEMP3, TEMP4 104* .. 105* .. Save statement .. 106 SAVE FIRST, SMALL1, SMINV1, SMALL2, SMINV2, EPS 107* .. 108* .. Intrinsic Functions .. 109 INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT 110* .. 111* .. Data statements .. 112 DATA FIRST / .TRUE. / 113* .. 114* .. Executable Statements .. 115* 116 INFO = 0 117* 118 NITER = 1 119 TAU = ZERO 120 IF( KNITER.EQ.2 ) THEN 121 IF( ORGATI ) THEN 122 TEMP = ( D( 3 )-D( 2 ) ) / TWO 123 C = RHO + Z( 1 ) / ( ( D( 1 )-D( 2 ) )-TEMP ) 124 A = C*( D( 2 )+D( 3 ) ) + Z( 2 ) + Z( 3 ) 125 B = C*D( 2 )*D( 3 ) + Z( 2 )*D( 3 ) + Z( 3 )*D( 2 ) 126 ELSE 127 TEMP = ( D( 1 )-D( 2 ) ) / TWO 128 C = RHO + Z( 3 ) / ( ( D( 3 )-D( 2 ) )-TEMP ) 129 A = C*( D( 1 )+D( 2 ) ) + Z( 1 ) + Z( 2 ) 130 B = C*D( 1 )*D( 2 ) + Z( 1 )*D( 2 ) + Z( 2 )*D( 1 ) 131 END IF 132 TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) ) 133 A = A / TEMP 134 B = B / TEMP 135 C = C / TEMP 136 OPS = OPS + 19 137 IF( C.EQ.ZERO ) THEN 138 TAU = B / A 139 OPS = OPS + 1 140 ELSE IF( A.LE.ZERO ) THEN 141 OPS = OPS + 8 142 TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 143 ELSE 144 OPS = OPS + 8 145 TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) 146 END IF 147 OPS = OPS + 9 148 TEMP = RHO + Z( 1 ) / ( D( 1 )-TAU ) + 149 $ Z( 2 ) / ( D( 2 )-TAU ) + Z( 3 ) / ( D( 3 )-TAU ) 150 IF( ABS( FINIT ).LE.ABS( TEMP ) ) 151 $ TAU = ZERO 152 END IF 153* 154* On first call to routine, get machine parameters for 155* possible scaling to avoid overflow 156* 157 IF( FIRST ) THEN 158 EPS = SLAMCH( 'Epsilon' ) 159 BASE = SLAMCH( 'Base' ) 160 SMALL1 = BASE**( INT( LOG( SLAMCH( 'SafMin' ) ) / LOG( BASE ) / 161 $ THREE ) ) 162 SMINV1 = ONE / SMALL1 163 SMALL2 = SMALL1*SMALL1 164 SMINV2 = SMINV1*SMINV1 165 FIRST = .FALSE. 166 END IF 167* 168* Determine if scaling of inputs necessary to avoid overflow 169* when computing 1/TEMP**3 170* 171 OPS = OPS + 2 172 IF( ORGATI ) THEN 173 TEMP = MIN( ABS( D( 2 )-TAU ), ABS( D( 3 )-TAU ) ) 174 ELSE 175 TEMP = MIN( ABS( D( 1 )-TAU ), ABS( D( 2 )-TAU ) ) 176 END IF 177 SCALE = .FALSE. 178 IF( TEMP.LE.SMALL1 ) THEN 179 SCALE = .TRUE. 180 IF( TEMP.LE.SMALL2 ) THEN 181* 182* Scale up by power of radix nearest 1/SAFMIN**(2/3) 183* 184 SCLFAC = SMINV2 185 SCLINV = SMALL2 186 ELSE 187* 188* Scale up by power of radix nearest 1/SAFMIN**(1/3) 189* 190 SCLFAC = SMINV1 191 SCLINV = SMALL1 192 END IF 193* 194* Scaling up safe because D, Z, TAU scaled elsewhere to be O(1) 195* 196 OPS = OPS + 7 197 DO 10 I = 1, 3 198 DSCALE( I ) = D( I )*SCLFAC 199 ZSCALE( I ) = Z( I )*SCLFAC 200 10 CONTINUE 201 TAU = TAU*SCLFAC 202 ELSE 203* 204* Copy D and Z to DSCALE and ZSCALE 205* 206 DO 20 I = 1, 3 207 DSCALE( I ) = D( I ) 208 ZSCALE( I ) = Z( I ) 209 20 CONTINUE 210 END IF 211* 212 FC = ZERO 213 DF = ZERO 214 DDF = ZERO 215 OPS = OPS + 11 216 DO 30 I = 1, 3 217 TEMP = ONE / ( DSCALE( I )-TAU ) 218 TEMP1 = ZSCALE( I )*TEMP 219 TEMP2 = TEMP1*TEMP 220 TEMP3 = TEMP2*TEMP 221 FC = FC + TEMP1 / DSCALE( I ) 222 DF = DF + TEMP2 223 DDF = DDF + TEMP3 224 30 CONTINUE 225 F = FINIT + TAU*FC 226* 227 IF( ABS( F ).LE.ZERO ) 228 $ GO TO 60 229* 230* Iteration begins 231* 232* It is not hard to see that 233* 234* 1) Iterations will go up monotonically 235* if FINIT < 0; 236* 237* 2) Iterations will go down monotonically 238* if FINIT > 0. 239* 240 ITER = NITER + 1 241* 242 DO 50 NITER = ITER, MAXIT 243* 244 OPS = OPS + 18 245 IF( ORGATI ) THEN 246 TEMP1 = DSCALE( 2 ) - TAU 247 TEMP2 = DSCALE( 3 ) - TAU 248 ELSE 249 TEMP1 = DSCALE( 1 ) - TAU 250 TEMP2 = DSCALE( 2 ) - TAU 251 END IF 252 A = ( TEMP1+TEMP2 )*F - TEMP1*TEMP2*DF 253 B = TEMP1*TEMP2*F 254 C = F - ( TEMP1+TEMP2 )*DF + TEMP1*TEMP2*DDF 255 TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) ) 256 A = A / TEMP 257 B = B / TEMP 258 C = C / TEMP 259 IF( C.EQ.ZERO ) THEN 260 OPS = OPS + 1 261 ETA = B / A 262 ELSE IF( A.LE.ZERO ) THEN 263 OPS = OPS + 8 264 ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 265 ELSE 266 OPS = OPS + 8 267 ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) 268 END IF 269 IF( F*ETA.GE.ZERO ) THEN 270 OPS = OPS + 1 271 ETA = -F / DF 272 END IF 273* 274 OPS = OPS + 1 275 TEMP = ETA + TAU 276 IF( ORGATI ) THEN 277 IF( ETA.GT.ZERO .AND. TEMP.GE.DSCALE( 3 ) ) THEN 278 OPS = OPS + 2 279 ETA = ( DSCALE( 3 )-TAU ) / TWO 280 END IF 281 IF( ETA.LT.ZERO .AND. TEMP.LE.DSCALE( 2 ) ) THEN 282 OPS = OPS + 2 283 ETA = ( DSCALE( 2 )-TAU ) / TWO 284 END IF 285 ELSE 286 IF( ETA.GT.ZERO .AND. TEMP.GE.DSCALE( 2 ) ) THEN 287 OPS = OPS + 2 288 ETA = ( DSCALE( 2 )-TAU ) / TWO 289 END IF 290 IF( ETA.LT.ZERO .AND. TEMP.LE.DSCALE( 1 ) ) THEN 291 OPS = OPS + 2 292 ETA = ( DSCALE( 1 )-TAU ) / TWO 293 END IF 294 END IF 295 OPS = OPS + 1 296 TAU = TAU + ETA 297* 298 FC = ZERO 299 ERRETM = ZERO 300 DF = ZERO 301 DDF = ZERO 302 OPS = OPS + 37 303 DO 40 I = 1, 3 304 TEMP = ONE / ( DSCALE( I )-TAU ) 305 TEMP1 = ZSCALE( I )*TEMP 306 TEMP2 = TEMP1*TEMP 307 TEMP3 = TEMP2*TEMP 308 TEMP4 = TEMP1 / DSCALE( I ) 309 FC = FC + TEMP4 310 ERRETM = ERRETM + ABS( TEMP4 ) 311 DF = DF + TEMP2 312 DDF = DDF + TEMP3 313 40 CONTINUE 314 F = FINIT + TAU*FC 315 ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) + 316 $ ABS( TAU )*DF 317 IF( ABS( F ).LE.EPS*ERRETM ) 318 $ GO TO 60 319 50 CONTINUE 320 INFO = 1 321 60 CONTINUE 322* 323* Undo scaling 324* 325 IF( SCALE ) THEN 326 OPS = OPS + 1 327 TAU = TAU*SCLINV 328 END IF 329 RETURN 330* 331* End of SLAED6 332* 333 END 334