LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slaed6.f
Go to the documentation of this file.
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
 All Files Functions