SRC/dlaed6.f

   1 *> \brief \b DLAED6 used by sstedc. Computes one Newton step in solution of the secular equation.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLAED6 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed6.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed6.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed6.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       LOGICAL            ORGATI
  25 *       INTEGER            INFO, KNITER
  26 *       DOUBLE PRECISION   FINIT, RHO, TAU
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   D( 3 ), Z( 3 )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> DLAED6 computes the positive or negative root (closest to the origin)
  39 *> of
  40 *>                  z(1)        z(2)        z(3)
  41 *> f(x) =   rho + --------- + ---------- + ---------
  42 *>                 d(1)-x      d(2)-x      d(3)-x
  43 *>
  44 *> It is assumed that
  45 *>
  46 *>       if ORGATI = .true. the root is between d(2) and d(3);
  47 *>       otherwise it is between d(1) and d(2)
  48 *>
  49 *> This routine will be called by DLAED4 when necessary. In most cases,
  50 *> the root sought is the smallest in magnitude, though it might not be
  51 *> in some extremely rare situations.
  52 *> \endverbatim
  53 *
  54 *  Arguments:
  55 *  ==========
  56 *
  57 *> \param[in] KNITER
  58 *> \verbatim
  59 *>          KNITER is INTEGER
  60 *>               Refer to DLAED4 for its significance.
  61 *> \endverbatim
  62 *>
  63 *> \param[in] ORGATI
  64 *> \verbatim
  65 *>          ORGATI is LOGICAL
  66 *>               If ORGATI is true, the needed root is between d(2) and
  67 *>               d(3); otherwise it is between d(1) and d(2).  See
  68 *>               DLAED4 for further details.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] RHO
  72 *> \verbatim
  73 *>          RHO is DOUBLE PRECISION
  74 *>               Refer to the equation f(x) above.
  75 *> \endverbatim
  76 *>
  77 *> \param[in] D
  78 *> \verbatim
  79 *>          D is DOUBLE PRECISION array, dimension (3)
  80 *>               D satisfies d(1) < d(2) < d(3).
  81 *> \endverbatim
  82 *>
  83 *> \param[in] Z
  84 *> \verbatim
  85 *>          Z is DOUBLE PRECISION array, dimension (3)
  86 *>               Each of the elements in z must be positive.
  87 *> \endverbatim
  88 *>
  89 *> \param[in] FINIT
  90 *> \verbatim
  91 *>          FINIT is DOUBLE PRECISION
  92 *>               The value of f at 0. It is more accurate than the one
  93 *>               evaluated inside this routine (if someone wants to do
  94 *>               so).
  95 *> \endverbatim
  96 *>
  97 *> \param[out] TAU
  98 *> \verbatim
  99 *>          TAU is DOUBLE PRECISION
 100 *>               The root of the equation f(x).
 101 *> \endverbatim
 102 *>
 103 *> \param[out] INFO
 104 *> \verbatim
 105 *>          INFO is INTEGER
 106 *>               = 0: successful exit
 107 *>               > 0: if INFO = 1, failure to converge
 108 *> \endverbatim
 109 *
 110 *  Authors:
 111 *  ========
 112 *
 113 *> \author Univ. of Tennessee
 114 *> \author Univ. of California Berkeley
 115 *> \author Univ. of Colorado Denver
 116 *> \author NAG Ltd.
 117 *
 118 *> \date November 2015
 119 *
 120 *> \ingroup auxOTHERcomputational
 121 *
 122 *> \par Further Details:
 123 *  =====================
 124 *>
 125 *> \verbatim
 126 *>
 127 *>  10/02/03: This version has a few statements commented out for thread
 128 *>  safety (machine parameters are computed on each entry). SJH.
 129 *>
 130 *>  05/10/06: Modified from a new version of Ren-Cang Li, use
 131 *>     Gragg-Thornton-Warner cubic convergent scheme for better stability.
 132 *> \endverbatim
 133 *
 134 *> \par Contributors:
 135 *  ==================
 136 *>
 137 *>     Ren-Cang Li, Computer Science Division, University of California
 138 *>     at Berkeley, USA
 139 *>
 140 *  =====================================================================
 141       SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
 142 *
 143 *  -- LAPACK computational routine (version 3.6.0) --
 144 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 145 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 146 *     November 2015
 147 *
 148 *     .. Scalar Arguments ..
 149       LOGICAL            ORGATI
 150       INTEGER            INFO, KNITER
 151       DOUBLE PRECISION   FINIT, RHO, TAU
 152 *     ..
 153 *     .. Array Arguments ..
 154       DOUBLE PRECISION   D( 3 ), Z( 3 )
 155 *     ..
 156 *
 157 *  =====================================================================
 158 *
 159 *     .. Parameters ..
 160       INTEGER            MAXIT
 161       PARAMETER          ( MAXIT = 40 )
 162       DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, EIGHT
 163       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
 164      $                   THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0 )
 165 *     ..
 166 *     .. External Functions ..
 167       DOUBLE PRECISION   DLAMCH
 168       EXTERNAL           DLAMCH
 169 *     ..
 170 *     .. Local Arrays ..
 171       DOUBLE PRECISION   DSCALE( 3 ), ZSCALE( 3 )
 172 *     ..
 173 *     .. Local Scalars ..
 174       LOGICAL            SCALE
 175       INTEGER            I, ITER, NITER
 176       DOUBLE PRECISION   A, B, BASE, C, DDF, DF, EPS, ERRETM, ETA, F,
 177      $                   FC, SCLFAC, SCLINV, SMALL1, SMALL2, SMINV1,
 178      $                   SMINV2, TEMP, TEMP1, TEMP2, TEMP3, TEMP4,
 179      $                   LBD, UBD
 180 *     ..
 181 *     .. Intrinsic Functions ..
 182       INTRINSIC          ABS, INT, LOG, MAX, MIN, SQRT
 183 *     ..
 184 *     .. Executable Statements ..
 185 *
 186       INFO = 0
 187 *
 188       IF( ORGATI ) THEN
 189          LBD = D(2)
 190          UBD = D(3)
 191       ELSE
 192          LBD = D(1)
 193          UBD = D(2)
 194       END IF
 195       IF( FINIT .LT. ZERO )THEN
 196          LBD = ZERO
 197       ELSE
 198          UBD = ZERO
 199       END IF
 200 *
 201       NITER = 1
 202       TAU = ZERO
 203       IF( KNITER.EQ.2 ) THEN
 204          IF( ORGATI ) THEN
 205             TEMP = ( D( 3 )-D( 2 ) ) / TWO
 206             C = RHO + Z( 1 ) / ( ( D( 1 )-D( 2 ) )-TEMP )
 207             A = C*( D( 2 )+D( 3 ) ) + Z( 2 ) + Z( 3 )
 208             B = C*D( 2 )*D( 3 ) + Z( 2 )*D( 3 ) + Z( 3 )*D( 2 )
 209          ELSE
 210             TEMP = ( D( 1 )-D( 2 ) ) / TWO
 211             C = RHO + Z( 3 ) / ( ( D( 3 )-D( 2 ) )-TEMP )
 212             A = C*( D( 1 )+D( 2 ) ) + Z( 1 ) + Z( 2 )
 213             B = C*D( 1 )*D( 2 ) + Z( 1 )*D( 2 ) + Z( 2 )*D( 1 )
 214          END IF
 215          TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
 216          A = A / TEMP
 217          B = B / TEMP
 218          C = C / TEMP
 219          IF( C.EQ.ZERO ) THEN
 220             TAU = B / A
 221          ELSE IF( A.LE.ZERO ) THEN
 222             TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
 223          ELSE
 224             TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
 225          END IF
 226          IF( TAU .LT. LBD .OR. TAU .GT. UBD )
 227      $      TAU = ( LBD+UBD )/TWO
 228          IF( D(1).EQ.TAU .OR. D(2).EQ.TAU .OR. D(3).EQ.TAU ) THEN
 229             TAU = ZERO
 230          ELSE
 231             TEMP = FINIT + TAU*Z(1)/( D(1)*( D( 1 )-TAU ) ) +
 232      $                     TAU*Z(2)/( D(2)*( D( 2 )-TAU ) ) +
 233      $                     TAU*Z(3)/( D(3)*( D( 3 )-TAU ) )
 234             IF( TEMP .LE. ZERO )THEN
 235                LBD = TAU
 236             ELSE
 237                UBD = TAU
 238             END IF
 239             IF( ABS( FINIT ).LE.ABS( TEMP ) )
 240      $         TAU = ZERO
 241          END IF
 242       END IF
 243 *
 244 *     get machine parameters for possible scaling to avoid overflow
 245 *
 246 *     modified by Sven: parameters SMALL1, SMINV1, SMALL2,
 247 *     SMINV2, EPS are not SAVEd anymore between one call to the
 248 *     others but recomputed at each call
 249 *
 250       EPS = DLAMCH( 'Epsilon' )
 251       BASE = DLAMCH( 'Base' )
 252       SMALL1 = BASE**( INT( LOG( DLAMCH( 'SafMin' ) ) / LOG( BASE ) /
 253      $         THREE ) )
 254       SMINV1 = ONE / SMALL1
 255       SMALL2 = SMALL1*SMALL1
 256       SMINV2 = SMINV1*SMINV1
 257 *
 258 *     Determine if scaling of inputs necessary to avoid overflow
 259 *     when computing 1/TEMP**3
 260 *
 261       IF( ORGATI ) THEN
 262          TEMP = MIN( ABS( D( 2 )-TAU ), ABS( D( 3 )-TAU ) )
 263       ELSE
 264          TEMP = MIN( ABS( D( 1 )-TAU ), ABS( D( 2 )-TAU ) )
 265       END IF
 266       SCALE = .FALSE.
 267       IF( TEMP.LE.SMALL1 ) THEN
 268          SCALE = .TRUE.
 269          IF( TEMP.LE.SMALL2 ) THEN
 270 *
 271 *        Scale up by power of radix nearest 1/SAFMIN**(2/3)
 272 *
 273             SCLFAC = SMINV2
 274             SCLINV = SMALL2
 275          ELSE
 276 *
 277 *        Scale up by power of radix nearest 1/SAFMIN**(1/3)
 278 *
 279             SCLFAC = SMINV1
 280             SCLINV = SMALL1
 281          END IF
 282 *
 283 *        Scaling up safe because D, Z, TAU scaled elsewhere to be O(1)
 284 *
 285          DO 10 I = 1, 3
 286             DSCALE( I ) = D( I )*SCLFAC
 287             ZSCALE( I ) = Z( I )*SCLFAC
 288    10    CONTINUE
 289          TAU = TAU*SCLFAC
 290          LBD = LBD*SCLFAC
 291          UBD = UBD*SCLFAC
 292       ELSE
 293 *
 294 *        Copy D and Z to DSCALE and ZSCALE
 295 *
 296          DO 20 I = 1, 3
 297             DSCALE( I ) = D( I )
 298             ZSCALE( I ) = Z( I )
 299    20    CONTINUE
 300       END IF
 301 *
 302       FC = ZERO
 303       DF = ZERO
 304       DDF = ZERO
 305       DO 30 I = 1, 3
 306          TEMP = ONE / ( DSCALE( I )-TAU )
 307          TEMP1 = ZSCALE( I )*TEMP
 308          TEMP2 = TEMP1*TEMP
 309          TEMP3 = TEMP2*TEMP
 310          FC = FC + TEMP1 / DSCALE( I )
 311          DF = DF + TEMP2
 312          DDF = DDF + TEMP3
 313    30 CONTINUE
 314       F = FINIT + TAU*FC
 315 *
 316       IF( ABS( F ).LE.ZERO )
 317      $   GO TO 60
 318       IF( F .LE. ZERO )THEN
 319          LBD = TAU
 320       ELSE
 321          UBD = TAU
 322       END IF
 323 *
 324 *        Iteration begins -- Use Gragg-Thornton-Warner cubic convergent
 325 *                            scheme
 326 *
 327 *     It is not hard to see that
 328 *
 329 *           1) Iterations will go up monotonically
 330 *              if FINIT < 0;
 331 *
 332 *           2) Iterations will go down monotonically
 333 *              if FINIT > 0.
 334 *
 335       ITER = NITER + 1
 336 *
 337       DO 50 NITER = ITER, MAXIT
 338 *
 339          IF( ORGATI ) THEN
 340             TEMP1 = DSCALE( 2 ) - TAU
 341             TEMP2 = DSCALE( 3 ) - TAU
 342          ELSE
 343             TEMP1 = DSCALE( 1 ) - TAU
 344             TEMP2 = DSCALE( 2 ) - TAU
 345          END IF
 346          A = ( TEMP1+TEMP2 )*F - TEMP1*TEMP2*DF
 347          B = TEMP1*TEMP2*F
 348          C = F - ( TEMP1+TEMP2 )*DF + TEMP1*TEMP2*DDF
 349          TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
 350          A = A / TEMP
 351          B = B / TEMP
 352          C = C / TEMP
 353          IF( C.EQ.ZERO ) THEN
 354             ETA = B / A
 355          ELSE IF( A.LE.ZERO ) THEN
 356             ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
 357          ELSE
 358             ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
 359          END IF
 360          IF( F*ETA.GE.ZERO ) THEN
 361             ETA = -F / DF
 362          END IF
 363 *
 364          TAU = TAU + ETA
 365          IF( TAU .LT. LBD .OR. TAU .GT. UBD )
 366      $      TAU = ( LBD + UBD )/TWO
 367 *
 368          FC = ZERO
 369          ERRETM = ZERO
 370          DF = ZERO
 371          DDF = ZERO
 372          DO 40 I = 1, 3
 373             IF ( ( DSCALE( I )-TAU ).NE.ZERO ) THEN
 374                TEMP = ONE / ( DSCALE( I )-TAU )
 375                TEMP1 = ZSCALE( I )*TEMP
 376                TEMP2 = TEMP1*TEMP
 377                TEMP3 = TEMP2*TEMP
 378                TEMP4 = TEMP1 / DSCALE( I )
 379                FC = FC + TEMP4
 380                ERRETM = ERRETM + ABS( TEMP4 )
 381                DF = DF + TEMP2
 382                DDF = DDF + TEMP3
 383             ELSE
 384                GO TO 60
 385             END IF
 386    40    CONTINUE
 387          F = FINIT + TAU*FC
 388          ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) +
 389      $            ABS( TAU )*DF
 390          IF( ( ABS( F ).LE.FOUR*EPS*ERRETM ) .OR.
 391      $      ( (UBD-LBD).LE.FOUR*EPS*ABS(TAU) )  )
 392      $      GO TO 60
 393          IF( F .LE. ZERO )THEN
 394             LBD = TAU
 395          ELSE
 396             UBD = TAU
 397          END IF
 398    50 CONTINUE
 399       INFO = 1
 400    60 CONTINUE
 401 *
 402 *     Undo scaling
 403 *
 404       IF( SCALE )
 405      $   TAU = TAU*SCLINV
 406       RETURN
 407 *
 408 *     End of DLAED6
 409 *
 410       END