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