*> \brief \b DLAED6 used by DSTEDC. Computes one Newton step in solution of the secular equation. * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download DLAED6 + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO ) * * .. Scalar Arguments .. * LOGICAL ORGATI * INTEGER INFO, KNITER * DOUBLE PRECISION FINIT, RHO, TAU * .. * .. Array Arguments .. * DOUBLE PRECISION D( 3 ), Z( 3 ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> DLAED6 computes the positive or negative root (closest to the origin) *> of *> z(1) z(2) z(3) *> f(x) = rho + --------- + ---------- + --------- *> d(1)-x d(2)-x d(3)-x *> *> It is assumed that *> *> if ORGATI = .true. the root is between d(2) and d(3); *> otherwise it is between d(1) and d(2) *> *> This routine will be called by DLAED4 when necessary. In most cases, *> the root sought is the smallest in magnitude, though it might not be *> in some extremely rare situations. *> \endverbatim * * Arguments: * ========== * *> \param[in] KNITER *> \verbatim *> KNITER is INTEGER *> Refer to DLAED4 for its significance. *> \endverbatim *> *> \param[in] ORGATI *> \verbatim *> ORGATI is LOGICAL *> If ORGATI is true, the needed root is between d(2) and *> d(3); otherwise it is between d(1) and d(2). See *> DLAED4 for further details. *> \endverbatim *> *> \param[in] RHO *> \verbatim *> RHO is DOUBLE PRECISION *> Refer to the equation f(x) above. *> \endverbatim *> *> \param[in] D *> \verbatim *> D is DOUBLE PRECISION array, dimension (3) *> D satisfies d(1) < d(2) < d(3). *> \endverbatim *> *> \param[in] Z *> \verbatim *> Z is DOUBLE PRECISION array, dimension (3) *> Each of the elements in z must be positive. *> \endverbatim *> *> \param[in] FINIT *> \verbatim *> FINIT is DOUBLE PRECISION *> The value of f at 0. It is more accurate than the one *> evaluated inside this routine (if someone wants to do *> so). *> \endverbatim *> *> \param[out] TAU *> \verbatim *> TAU is DOUBLE PRECISION *> The root of the equation f(x). *> \endverbatim *> *> \param[out] INFO *> \verbatim *> INFO is INTEGER *> = 0: successful exit *> > 0: if INFO = 1, failure to converge *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \ingroup auxOTHERcomputational * *> \par Further Details: * ===================== *> *> \verbatim *> *> 10/02/03: This version has a few statements commented out for thread *> safety (machine parameters are computed on each entry). SJH. *> *> 05/10/06: Modified from a new version of Ren-Cang Li, use *> Gragg-Thornton-Warner cubic convergent scheme for better stability. *> \endverbatim * *> \par Contributors: * ================== *> *> Ren-Cang Li, Computer Science Division, University of California *> at Berkeley, USA *> * ===================================================================== SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO ) * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * * .. Scalar Arguments .. LOGICAL ORGATI INTEGER INFO, KNITER DOUBLE PRECISION FINIT, RHO, TAU * .. * .. Array Arguments .. DOUBLE PRECISION D( 3 ), Z( 3 ) * .. * * ===================================================================== * * .. Parameters .. INTEGER MAXIT PARAMETER ( MAXIT = 40 ) DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, \$ THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0 ) * .. * .. External Functions .. DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH * .. * .. Local Arrays .. DOUBLE PRECISION DSCALE( 3 ), ZSCALE( 3 ) * .. * .. Local Scalars .. LOGICAL SCALE INTEGER I, ITER, NITER DOUBLE PRECISION A, B, BASE, C, DDF, DF, EPS, ERRETM, ETA, F, \$ FC, SCLFAC, SCLINV, SMALL1, SMALL2, SMINV1, \$ SMINV2, TEMP, TEMP1, TEMP2, TEMP3, TEMP4, \$ LBD, UBD * .. * .. Intrinsic Functions .. INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT * .. * .. Executable Statements .. * INFO = 0 * IF( ORGATI ) THEN LBD = D(2) UBD = D(3) ELSE LBD = D(1) UBD = D(2) END IF IF( FINIT .LT. ZERO )THEN LBD = ZERO ELSE UBD = ZERO END IF * NITER = 1 TAU = ZERO IF( KNITER.EQ.2 ) THEN IF( ORGATI ) THEN TEMP = ( D( 3 )-D( 2 ) ) / TWO C = RHO + Z( 1 ) / ( ( D( 1 )-D( 2 ) )-TEMP ) A = C*( D( 2 )+D( 3 ) ) + Z( 2 ) + Z( 3 ) B = C*D( 2 )*D( 3 ) + Z( 2 )*D( 3 ) + Z( 3 )*D( 2 ) ELSE TEMP = ( D( 1 )-D( 2 ) ) / TWO C = RHO + Z( 3 ) / ( ( D( 3 )-D( 2 ) )-TEMP ) A = C*( D( 1 )+D( 2 ) ) + Z( 1 ) + Z( 2 ) B = C*D( 1 )*D( 2 ) + Z( 1 )*D( 2 ) + Z( 2 )*D( 1 ) END IF TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) ) A = A / TEMP B = B / TEMP C = C / TEMP IF( C.EQ.ZERO ) THEN TAU = B / A ELSE IF( A.LE.ZERO ) THEN TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) ELSE TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) END IF IF( TAU .LT. LBD .OR. TAU .GT. UBD ) \$ TAU = ( LBD+UBD )/TWO IF( D(1).EQ.TAU .OR. D(2).EQ.TAU .OR. D(3).EQ.TAU ) THEN TAU = ZERO ELSE TEMP = FINIT + TAU*Z(1)/( D(1)*( D( 1 )-TAU ) ) + \$ TAU*Z(2)/( D(2)*( D( 2 )-TAU ) ) + \$ TAU*Z(3)/( D(3)*( D( 3 )-TAU ) ) IF( TEMP .LE. ZERO )THEN LBD = TAU ELSE UBD = TAU END IF IF( ABS( FINIT ).LE.ABS( TEMP ) ) \$ TAU = ZERO END IF END IF * * get machine parameters for possible scaling to avoid overflow * * modified by Sven: parameters SMALL1, SMINV1, SMALL2, * SMINV2, EPS are not SAVEd anymore between one call to the * others but recomputed at each call * EPS = DLAMCH( 'Epsilon' ) BASE = DLAMCH( 'Base' ) SMALL1 = BASE**( INT( LOG( DLAMCH( 'SafMin' ) ) / LOG( BASE ) / \$ THREE ) ) SMINV1 = ONE / SMALL1 SMALL2 = SMALL1*SMALL1 SMINV2 = SMINV1*SMINV1 * * Determine if scaling of inputs necessary to avoid overflow * when computing 1/TEMP**3 * IF( ORGATI ) THEN TEMP = MIN( ABS( D( 2 )-TAU ), ABS( D( 3 )-TAU ) ) ELSE TEMP = MIN( ABS( D( 1 )-TAU ), ABS( D( 2 )-TAU ) ) END IF SCALE = .FALSE. IF( TEMP.LE.SMALL1 ) THEN SCALE = .TRUE. IF( TEMP.LE.SMALL2 ) THEN * * Scale up by power of radix nearest 1/SAFMIN**(2/3) * SCLFAC = SMINV2 SCLINV = SMALL2 ELSE * * Scale up by power of radix nearest 1/SAFMIN**(1/3) * SCLFAC = SMINV1 SCLINV = SMALL1 END IF * * Scaling up safe because D, Z, TAU scaled elsewhere to be O(1) * DO 10 I = 1, 3 DSCALE( I ) = D( I )*SCLFAC ZSCALE( I ) = Z( I )*SCLFAC 10 CONTINUE TAU = TAU*SCLFAC LBD = LBD*SCLFAC UBD = UBD*SCLFAC ELSE * * Copy D and Z to DSCALE and ZSCALE * DO 20 I = 1, 3 DSCALE( I ) = D( I ) ZSCALE( I ) = Z( I ) 20 CONTINUE END IF * FC = ZERO DF = ZERO DDF = ZERO DO 30 I = 1, 3 TEMP = ONE / ( DSCALE( I )-TAU ) TEMP1 = ZSCALE( I )*TEMP TEMP2 = TEMP1*TEMP TEMP3 = TEMP2*TEMP FC = FC + TEMP1 / DSCALE( I ) DF = DF + TEMP2 DDF = DDF + TEMP3 30 CONTINUE F = FINIT + TAU*FC * IF( ABS( F ).LE.ZERO ) \$ GO TO 60 IF( F .LE. ZERO )THEN LBD = TAU ELSE UBD = TAU END IF * * Iteration begins -- Use Gragg-Thornton-Warner cubic convergent * scheme * * It is not hard to see that * * 1) Iterations will go up monotonically * if FINIT < 0; * * 2) Iterations will go down monotonically * if FINIT > 0. * ITER = NITER + 1 * DO 50 NITER = ITER, MAXIT * IF( ORGATI ) THEN TEMP1 = DSCALE( 2 ) - TAU TEMP2 = DSCALE( 3 ) - TAU ELSE TEMP1 = DSCALE( 1 ) - TAU TEMP2 = DSCALE( 2 ) - TAU END IF A = ( TEMP1+TEMP2 )*F - TEMP1*TEMP2*DF B = TEMP1*TEMP2*F C = F - ( TEMP1+TEMP2 )*DF + TEMP1*TEMP2*DDF TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) ) A = A / TEMP B = B / TEMP C = C / TEMP IF( C.EQ.ZERO ) THEN ETA = B / A ELSE IF( A.LE.ZERO ) THEN ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) ELSE ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) END IF IF( F*ETA.GE.ZERO ) THEN ETA = -F / DF END IF * TAU = TAU + ETA IF( TAU .LT. LBD .OR. TAU .GT. UBD ) \$ TAU = ( LBD + UBD )/TWO * FC = ZERO ERRETM = ZERO DF = ZERO DDF = ZERO DO 40 I = 1, 3 IF ( ( DSCALE( I )-TAU ).NE.ZERO ) THEN TEMP = ONE / ( DSCALE( I )-TAU ) TEMP1 = ZSCALE( I )*TEMP TEMP2 = TEMP1*TEMP TEMP3 = TEMP2*TEMP TEMP4 = TEMP1 / DSCALE( I ) FC = FC + TEMP4 ERRETM = ERRETM + ABS( TEMP4 ) DF = DF + TEMP2 DDF = DDF + TEMP3 ELSE GO TO 60 END IF 40 CONTINUE F = FINIT + TAU*FC ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) + \$ ABS( TAU )*DF IF( ( ABS( F ).LE.FOUR*EPS*ERRETM ) .OR. \$ ( (UBD-LBD).LE.FOUR*EPS*ABS(TAU) ) ) \$ GO TO 60 IF( F .LE. ZERO )THEN LBD = TAU ELSE UBD = TAU END IF 50 CONTINUE INFO = 1 60 CONTINUE * * Undo scaling * IF( SCALE ) \$ TAU = TAU*SCLINV RETURN * * End of DLAED6 * END