#include "ESMF_LapackBlas.inc" !> \brief \b DLAED6 ! ! =========== DOCUMENTATION =========== ! ! Online html documentation available at ! http://www.netlib.org/lapack/explore-html/ ! !> \htmlonly !> Download DLAED6 + dependencies !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed6.f"> !> [TGZ]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed6.f"> !> [ZIP]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed6.f"> !> [TXT]</a> !> \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. ! !> \date November 2011 ! !> \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 (version 3.4.0) -- ! -- LAPACK is a software package provided by Univ. of Tennessee, -- ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! November 2011 ! ! .. 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 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 40 CONTINUE F = FINIT + TAU*FC ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) + & & ABS( TAU )*DF IF( ABS( F ).LE.EPS*ERRETM ) & & 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