 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * February 2007 00007 * 00008 * .. Scalar Arguments .. 00009 LOGICAL ORGATI 00010 INTEGER INFO, KNITER 00011 DOUBLE PRECISION FINIT, RHO, TAU 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION D( 3 ), Z( 3 ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DLAED6 computes the positive or negative root (closest to the origin) 00021 * of 00022 * z(1) z(2) z(3) 00023 * f(x) = rho + --------- + ---------- + --------- 00024 * d(1)-x d(2)-x d(3)-x 00025 * 00026 * It is assumed that 00027 * 00028 * if ORGATI = .true. the root is between d(2) and d(3); 00029 * otherwise it is between d(1) and d(2) 00030 * 00031 * This routine will be called by DLAED4 when necessary. In most cases, 00032 * the root sought is the smallest in magnitude, though it might not be 00033 * in some extremely rare situations. 00034 * 00035 * Arguments 00036 * ========= 00037 * 00038 * KNITER (input) INTEGER 00039 * Refer to DLAED4 for its significance. 00040 * 00041 * ORGATI (input) LOGICAL 00042 * If ORGATI is true, the needed root is between d(2) and 00043 * d(3); otherwise it is between d(1) and d(2). See 00044 * DLAED4 for further details. 00045 * 00046 * RHO (input) DOUBLE PRECISION 00047 * Refer to the equation f(x) above. 00048 * 00049 * D (input) DOUBLE PRECISION array, dimension (3) 00050 * D satisfies d(1) < d(2) < d(3). 00051 * 00052 * Z (input) DOUBLE PRECISION array, dimension (3) 00053 * Each of the elements in z must be positive. 00054 * 00055 * FINIT (input) DOUBLE PRECISION 00056 * The value of f at 0. It is more accurate than the one 00057 * evaluated inside this routine (if someone wants to do 00058 * so). 00059 * 00060 * TAU (output) DOUBLE PRECISION 00061 * The root of the equation f(x). 00062 * 00063 * INFO (output) INTEGER 00064 * = 0: successful exit 00065 * > 0: if INFO = 1, failure to converge 00066 * 00067 * Further Details 00068 * =============== 00069 * 00070 * 30/06/99: Based on contributions by 00071 * Ren-Cang Li, Computer Science Division, University of California 00072 * at Berkeley, USA 00073 * 00074 * 10/02/03: This version has a few statements commented out for thread 00075 * safety (machine parameters are computed on each entry). SJH. 00076 * 00077 * 05/10/06: Modified from a new version of Ren-Cang Li, use 00078 * Gragg-Thornton-Warner cubic convergent scheme for better stability. 00079 * 00080 * ===================================================================== 00081 * 00082 * .. Parameters .. 00083 INTEGER MAXIT 00084 PARAMETER ( MAXIT = 40 ) 00085 DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT 00086 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 00087 $ THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0 ) 00088 * .. 00089 * .. External Functions .. 00090 DOUBLE PRECISION DLAMCH 00091 EXTERNAL DLAMCH 00092 * .. 00093 * .. Local Arrays .. 00094 DOUBLE PRECISION DSCALE( 3 ), ZSCALE( 3 ) 00095 * .. 00096 * .. Local Scalars .. 00097 LOGICAL SCALE 00098 INTEGER I, ITER, NITER 00099 DOUBLE PRECISION A, B, BASE, C, DDF, DF, EPS, ERRETM, ETA, F, 00100 $ FC, SCLFAC, SCLINV, SMALL1, SMALL2, SMINV1, 00101 $ SMINV2, TEMP, TEMP1, TEMP2, TEMP3, TEMP4, 00102 $ LBD, UBD 00103 * .. 00104 * .. Intrinsic Functions .. 00105 INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT 00106 * .. 00107 * .. Executable Statements .. 00108 * 00109 INFO = 0 00110 * 00111 IF( ORGATI ) THEN 00112 LBD = D(2) 00113 UBD = D(3) 00114 ELSE 00115 LBD = D(1) 00116 UBD = D(2) 00117 END IF 00118 IF( FINIT .LT. ZERO )THEN 00119 LBD = ZERO 00120 ELSE 00121 UBD = ZERO 00122 END IF 00123 * 00124 NITER = 1 00125 TAU = ZERO 00126 IF( KNITER.EQ.2 ) THEN 00127 IF( ORGATI ) THEN 00128 TEMP = ( D( 3 )-D( 2 ) ) / TWO 00129 C = RHO + Z( 1 ) / ( ( D( 1 )-D( 2 ) )-TEMP ) 00130 A = C*( D( 2 )+D( 3 ) ) + Z( 2 ) + Z( 3 ) 00131 B = C*D( 2 )*D( 3 ) + Z( 2 )*D( 3 ) + Z( 3 )*D( 2 ) 00132 ELSE 00133 TEMP = ( D( 1 )-D( 2 ) ) / TWO 00134 C = RHO + Z( 3 ) / ( ( D( 3 )-D( 2 ) )-TEMP ) 00135 A = C*( D( 1 )+D( 2 ) ) + Z( 1 ) + Z( 2 ) 00136 B = C*D( 1 )*D( 2 ) + Z( 1 )*D( 2 ) + Z( 2 )*D( 1 ) 00137 END IF 00138 TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) ) 00139 A = A / TEMP 00140 B = B / TEMP 00141 C = C / TEMP 00142 IF( C.EQ.ZERO ) THEN 00143 TAU = B / A 00144 ELSE IF( A.LE.ZERO ) THEN 00145 TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 00146 ELSE 00147 TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) 00148 END IF 00149 IF( TAU .LT. LBD .OR. TAU .GT. UBD ) 00150 $ TAU = ( LBD+UBD )/TWO 00151 IF( D(1).EQ.TAU .OR. D(2).EQ.TAU .OR. D(3).EQ.TAU ) THEN 00152 TAU = ZERO 00153 ELSE 00154 TEMP = FINIT + TAU*Z(1)/( D(1)*( D( 1 )-TAU ) ) + 00155 $ TAU*Z(2)/( D(2)*( D( 2 )-TAU ) ) + 00156 $ TAU*Z(3)/( D(3)*( D( 3 )-TAU ) ) 00157 IF( TEMP .LE. ZERO )THEN 00158 LBD = TAU 00159 ELSE 00160 UBD = TAU 00161 END IF 00162 IF( ABS( FINIT ).LE.ABS( TEMP ) ) 00163 $ TAU = ZERO 00164 END IF 00165 END IF 00166 * 00167 * get machine parameters for possible scaling to avoid overflow 00168 * 00169 * modified by Sven: parameters SMALL1, SMINV1, SMALL2, 00170 * SMINV2, EPS are not SAVEd anymore between one call to the 00171 * others but recomputed at each call 00172 * 00173 EPS = DLAMCH( 'Epsilon' ) 00174 BASE = DLAMCH( 'Base' ) 00175 SMALL1 = BASE**( INT( LOG( DLAMCH( 'SafMin' ) ) / LOG( BASE ) / 00176 $ THREE ) ) 00177 SMINV1 = ONE / SMALL1 00178 SMALL2 = SMALL1*SMALL1 00179 SMINV2 = SMINV1*SMINV1 00180 * 00181 * Determine if scaling of inputs necessary to avoid overflow 00182 * when computing 1/TEMP**3 00183 * 00184 IF( ORGATI ) THEN 00185 TEMP = MIN( ABS( D( 2 )-TAU ), ABS( D( 3 )-TAU ) ) 00186 ELSE 00187 TEMP = MIN( ABS( D( 1 )-TAU ), ABS( D( 2 )-TAU ) ) 00188 END IF 00189 SCALE = .FALSE. 00190 IF( TEMP.LE.SMALL1 ) THEN 00191 SCALE = .TRUE. 00192 IF( TEMP.LE.SMALL2 ) THEN 00193 * 00194 * Scale up by power of radix nearest 1/SAFMIN**(2/3) 00195 * 00196 SCLFAC = SMINV2 00197 SCLINV = SMALL2 00198 ELSE 00199 * 00200 * Scale up by power of radix nearest 1/SAFMIN**(1/3) 00201 * 00202 SCLFAC = SMINV1 00203 SCLINV = SMALL1 00204 END IF 00205 * 00206 * Scaling up safe because D, Z, TAU scaled elsewhere to be O(1) 00207 * 00208 DO 10 I = 1, 3 00209 DSCALE( I ) = D( I )*SCLFAC 00210 ZSCALE( I ) = Z( I )*SCLFAC 00211 10 CONTINUE 00212 TAU = TAU*SCLFAC 00213 LBD = LBD*SCLFAC 00214 UBD = UBD*SCLFAC 00215 ELSE 00216 * 00217 * Copy D and Z to DSCALE and ZSCALE 00218 * 00219 DO 20 I = 1, 3 00220 DSCALE( I ) = D( I ) 00221 ZSCALE( I ) = Z( I ) 00222 20 CONTINUE 00223 END IF 00224 * 00225 FC = ZERO 00226 DF = ZERO 00227 DDF = ZERO 00228 DO 30 I = 1, 3 00229 TEMP = ONE / ( DSCALE( I )-TAU ) 00230 TEMP1 = ZSCALE( I )*TEMP 00231 TEMP2 = TEMP1*TEMP 00232 TEMP3 = TEMP2*TEMP 00233 FC = FC + TEMP1 / DSCALE( I ) 00234 DF = DF + TEMP2 00235 DDF = DDF + TEMP3 00236 30 CONTINUE 00237 F = FINIT + TAU*FC 00238 * 00239 IF( ABS( F ).LE.ZERO ) 00240 $ GO TO 60 00241 IF( F .LE. ZERO )THEN 00242 LBD = TAU 00243 ELSE 00244 UBD = TAU 00245 END IF 00246 * 00247 * Iteration begins -- Use Gragg-Thornton-Warner cubic convergent 00248 * scheme 00249 * 00250 * It is not hard to see that 00251 * 00252 * 1) Iterations will go up monotonically 00253 * if FINIT < 0; 00254 * 00255 * 2) Iterations will go down monotonically 00256 * if FINIT > 0. 00257 * 00258 ITER = NITER + 1 00259 * 00260 DO 50 NITER = ITER, MAXIT 00261 * 00262 IF( ORGATI ) THEN 00263 TEMP1 = DSCALE( 2 ) - TAU 00264 TEMP2 = DSCALE( 3 ) - TAU 00265 ELSE 00266 TEMP1 = DSCALE( 1 ) - TAU 00267 TEMP2 = DSCALE( 2 ) - TAU 00268 END IF 00269 A = ( TEMP1+TEMP2 )*F - TEMP1*TEMP2*DF 00270 B = TEMP1*TEMP2*F 00271 C = F - ( TEMP1+TEMP2 )*DF + TEMP1*TEMP2*DDF 00272 TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) ) 00273 A = A / TEMP 00274 B = B / TEMP 00275 C = C / TEMP 00276 IF( C.EQ.ZERO ) THEN 00277 ETA = B / A 00278 ELSE IF( A.LE.ZERO ) THEN 00279 ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 00280 ELSE 00281 ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) 00282 END IF 00283 IF( F*ETA.GE.ZERO ) THEN 00284 ETA = -F / DF 00285 END IF 00286 * 00287 TAU = TAU + ETA 00288 IF( TAU .LT. LBD .OR. TAU .GT. UBD ) 00289 $ TAU = ( LBD + UBD )/TWO 00290 * 00291 FC = ZERO 00292 ERRETM = ZERO 00293 DF = ZERO 00294 DDF = ZERO 00295 DO 40 I = 1, 3 00296 TEMP = ONE / ( DSCALE( I )-TAU ) 00297 TEMP1 = ZSCALE( I )*TEMP 00298 TEMP2 = TEMP1*TEMP 00299 TEMP3 = TEMP2*TEMP 00300 TEMP4 = TEMP1 / DSCALE( I ) 00301 FC = FC + TEMP4 00302 ERRETM = ERRETM + ABS( TEMP4 ) 00303 DF = DF + TEMP2 00304 DDF = DDF + TEMP3 00305 40 CONTINUE 00306 F = FINIT + TAU*FC 00307 ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) + 00308 $ ABS( TAU )*DF 00309 IF( ABS( F ).LE.EPS*ERRETM ) 00310 $ GO TO 60 00311 IF( F .LE. ZERO )THEN 00312 LBD = TAU 00313 ELSE 00314 UBD = TAU 00315 END IF 00316 50 CONTINUE 00317 INFO = 1 00318 60 CONTINUE 00319 * 00320 * Undo scaling 00321 * 00322 IF( SCALE ) 00323 $ TAU = TAU*SCLINV 00324 RETURN 00325 * 00326 * End of DLAED6 00327 * 00328 END
1.7.3 
