#include "ESMF_LapackBlas.inc" !> \brief \b DLAED4 used by sstedc. Finds a single root of the secular equation. ! ! =========== DOCUMENTATION =========== ! ! Online html documentation available at ! http://www.netlib.org/lapack/explore-html/ ! !> \htmlonly !> Download DLAED4 + dependencies !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed4.f"> !> [TGZ]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed4.f"> !> [ZIP]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed4.f"> !> [TXT]</a> !> \endhtmlonly ! ! Definition: ! =========== ! ! SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO ) ! ! .. Scalar Arguments .. ! INTEGER I, INFO, N ! DOUBLE PRECISION DLAM, RHO ! .. ! .. Array Arguments .. ! DOUBLE PRECISION D( * ), DELTA( * ), Z( * ) ! .. ! ! !> \par Purpose: ! ============= !> !> \verbatim !> !> This subroutine computes the I-th updated eigenvalue of a symmetric !> rank-one modification to a diagonal matrix whose elements are !> given in the array d, and that !> !> D(i) < D(j) for i < j !> !> and that RHO > 0. This is arranged by the calling routine, and is !> no loss in generality. The rank-one modified system is thus !> !> diag( D ) + RHO * Z * Z_transpose. !> !> where we assume the Euclidean norm of Z is 1. !> !> The method consists of approximating the rational functions in the !> secular equation by simpler interpolating rational functions. !> \endverbatim ! ! Arguments: ! ========== ! !> \param[in] N !> \verbatim !> N is INTEGER !> The length of all arrays. !> \endverbatim !> !> \param[in] I !> \verbatim !> I is INTEGER !> The index of the eigenvalue to be computed. 1 <= I <= N. !> \endverbatim !> !> \param[in] D !> \verbatim !> D is DOUBLE PRECISION array, dimension (N) !> The original eigenvalues. It is assumed that they are in !> order, D(I) < D(J) for I < J. !> \endverbatim !> !> \param[in] Z !> \verbatim !> Z is DOUBLE PRECISION array, dimension (N) !> The components of the updating vector. !> \endverbatim !> !> \param[out] DELTA !> \verbatim !> DELTA is DOUBLE PRECISION array, dimension (N) !> If N .GT. 2, DELTA contains (D(j) - lambda_I) in its j-th !> component. If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5 !> for detail. The vector DELTA contains the information necessary !> to construct the eigenvectors by DLAED3 and DLAED9. !> \endverbatim !> !> \param[in] RHO !> \verbatim !> RHO is DOUBLE PRECISION !> The scalar in the symmetric updating formula. !> \endverbatim !> !> \param[out] DLAM !> \verbatim !> DLAM is DOUBLE PRECISION !> The computed lambda_I, the I-th updated eigenvalue. !> \endverbatim !> !> \param[out] INFO !> \verbatim !> INFO is INTEGER !> = 0: successful exit !> > 0: if INFO = 1, the updating process failed. !> \endverbatim ! !> \par Internal Parameters: ! ========================= !> !> \verbatim !> Logical variable ORGATI (origin-at-i?) is used for distinguishing !> whether D(i) or D(i+1) is treated as the origin. !> !> ORGATI = .true. origin at i !> ORGATI = .false. origin at i+1 !> !> Logical variable SWTCH3 (switch-for-3-poles?) is for noting !> if we are working with THREE poles! !> !> MAXIT is the maximum number of iterations allowed for each !> eigenvalue. !> \endverbatim ! ! Authors: ! ======== ! !> \author Univ. of Tennessee !> \author Univ. of California Berkeley !> \author Univ. of Colorado Denver !> \author NAG Ltd. ! !> \date December 2016 ! !> \ingroup auxOTHERcomputational ! !> \par Contributors: ! ================== !> !> Ren-Cang Li, Computer Science Division, University of California !> at Berkeley, USA !> ! ===================================================================== SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO ) ! ! -- LAPACK computational routine (version 3.7.0) -- ! -- LAPACK is a software package provided by Univ. of Tennessee, -- ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! December 2016 ! ! .. Scalar Arguments .. INTEGER I, INFO, N DOUBLE PRECISION DLAM, RHO ! .. ! .. Array Arguments .. DOUBLE PRECISION D( * ), DELTA( * ), Z( * ) ! .. ! ! ===================================================================== ! ! .. Parameters .. INTEGER MAXIT PARAMETER ( MAXIT = 30 ) DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, & THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0, & TEN = 10.0D0 ) ! .. ! .. Local Scalars .. LOGICAL ORGATI, SWTCH, SWTCH3 INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER DOUBLE PRECISION A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW, & EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI, & RHOINV, TAU, TEMP, TEMP1, W ! .. ! .. Local Arrays .. DOUBLE PRECISION ZZ( 3 ) ! .. ! .. External Functions .. DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH ! .. ! .. External Subroutines .. EXTERNAL DLAED5, DLAED6 ! .. ! .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT ! .. ! .. Executable Statements .. ! ! Since this routine is called in an inner loop, we do no argument ! checking. ! ! Quick return for N=1 and 2. ! INFO = 0 IF( N.EQ.1 ) THEN ! ! Presumably, I=1 upon entry ! DLAM = D( 1 ) + RHO*Z( 1 )*Z( 1 ) DELTA( 1 ) = ONE RETURN END IF IF( N.EQ.2 ) THEN CALL DLAED5( I, D, Z, DELTA, RHO, DLAM ) RETURN END IF ! ! Compute machine epsilon ! EPS = DLAMCH( 'Epsilon' ) RHOINV = ONE / RHO ! ! The case I = N ! IF( I.EQ.N ) THEN ! ! Initialize some basic variables ! II = N - 1 NITER = 1 ! ! Calculate initial guess ! MIDPT = RHO / TWO ! ! If ||Z||_2 is not one, then TEMP should be set to ! RHO * ||Z||_2^2 / TWO ! DO 10 J = 1, N DELTA( J ) = ( D( J )-D( I ) ) - MIDPT 10 CONTINUE ! PSI = ZERO DO 20 J = 1, N - 2 PSI = PSI + Z( J )*Z( J ) / DELTA( J ) 20 CONTINUE ! C = RHOINV + PSI W = C + Z( II )*Z( II ) / DELTA( II ) + & Z( N )*Z( N ) / DELTA( N ) ! IF( W.LE.ZERO ) THEN TEMP = Z( N-1 )*Z( N-1 ) / ( D( N )-D( N-1 )+RHO ) + & Z( N )*Z( N ) / RHO IF( C.LE.TEMP ) THEN TAU = RHO ELSE DEL = D( N ) - D( N-1 ) A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N ) B = Z( N )*Z( N )*DEL IF( A.LT.ZERO ) THEN TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A ) ELSE TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C ) END IF END IF ! ! It can be proved that ! D(N)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO ! DLTLB = MIDPT DLTUB = RHO ELSE DEL = D( N ) - D( N-1 ) A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N ) B = Z( N )*Z( N )*DEL IF( A.LT.ZERO ) THEN TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A ) ELSE TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C ) END IF ! ! It can be proved that ! D(N) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2 ! DLTLB = ZERO DLTUB = MIDPT END IF ! DO 30 J = 1, N DELTA( J ) = ( D( J )-D( I ) ) - TAU 30 CONTINUE ! ! Evaluate PSI and the derivative DPSI ! DPSI = ZERO PSI = ZERO ERRETM = ZERO DO 40 J = 1, II TEMP = Z( J ) / DELTA( J ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 40 CONTINUE ERRETM = ABS( ERRETM ) ! ! Evaluate PHI and the derivative DPHI ! TEMP = Z( N ) / DELTA( N ) PHI = Z( N )*TEMP DPHI = TEMP*TEMP ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV + & ABS( TAU )*( DPSI+DPHI ) ! W = RHOINV + PHI + PSI ! ! Test for convergence ! IF( ABS( W ).LE.EPS*ERRETM ) THEN DLAM = D( I ) + TAU GO TO 250 END IF ! IF( W.LE.ZERO ) THEN DLTLB = MAX( DLTLB, TAU ) ELSE DLTUB = MIN( DLTUB, TAU ) END IF ! ! Calculate the new step ! NITER = NITER + 1 C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI A = ( DELTA( N-1 )+DELTA( N ) )*W - & DELTA( N-1 )*DELTA( N )*( DPSI+DPHI ) B = DELTA( N-1 )*DELTA( N )*W IF( C.LT.ZERO ) & C = ABS( C ) IF( C.EQ.ZERO ) THEN ! ETA = B/A ! ETA = RHO - TAU ETA = DLTUB - TAU ELSE IF( A.GE.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 ! ! Note, eta should be positive if w is negative, and ! eta should be negative otherwise. However, ! if for some reason caused by roundoff, eta*w > 0, ! we simply use one Newton step instead. This way ! will guarantee eta*w < 0. ! IF( W*ETA.GT.ZERO ) & ETA = -W / ( DPSI+DPHI ) TEMP = TAU + ETA IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN IF( W.LT.ZERO ) THEN ETA = ( DLTUB-TAU ) / TWO ELSE ETA = ( DLTLB-TAU ) / TWO END IF END IF DO 50 J = 1, N DELTA( J ) = DELTA( J ) - ETA 50 CONTINUE ! TAU = TAU + ETA ! ! Evaluate PSI and the derivative DPSI ! DPSI = ZERO PSI = ZERO ERRETM = ZERO DO 60 J = 1, II TEMP = Z( J ) / DELTA( J ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 60 CONTINUE ERRETM = ABS( ERRETM ) ! ! Evaluate PHI and the derivative DPHI ! TEMP = Z( N ) / DELTA( N ) PHI = Z( N )*TEMP DPHI = TEMP*TEMP ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV + & ABS( TAU )*( DPSI+DPHI ) ! W = RHOINV + PHI + PSI ! ! Main loop to update the values of the array DELTA ! ITER = NITER + 1 ! DO 90 NITER = ITER, MAXIT ! ! Test for convergence ! IF( ABS( W ).LE.EPS*ERRETM ) THEN DLAM = D( I ) + TAU GO TO 250 END IF ! IF( W.LE.ZERO ) THEN DLTLB = MAX( DLTLB, TAU ) ELSE DLTUB = MIN( DLTUB, TAU ) END IF ! ! Calculate the new step ! C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI A = ( DELTA( N-1 )+DELTA( N ) )*W - & DELTA( N-1 )*DELTA( N )*( DPSI+DPHI ) B = DELTA( N-1 )*DELTA( N )*W IF( A.GE.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 ! ! Note, eta should be positive if w is negative, and ! eta should be negative otherwise. However, ! if for some reason caused by roundoff, eta*w > 0, ! we simply use one Newton step instead. This way ! will guarantee eta*w < 0. ! IF( W*ETA.GT.ZERO ) & ETA = -W / ( DPSI+DPHI ) TEMP = TAU + ETA IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN IF( W.LT.ZERO ) THEN ETA = ( DLTUB-TAU ) / TWO ELSE ETA = ( DLTLB-TAU ) / TWO END IF END IF DO 70 J = 1, N DELTA( J ) = DELTA( J ) - ETA 70 CONTINUE ! TAU = TAU + ETA ! ! Evaluate PSI and the derivative DPSI ! DPSI = ZERO PSI = ZERO ERRETM = ZERO DO 80 J = 1, II TEMP = Z( J ) / DELTA( J ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 80 CONTINUE ERRETM = ABS( ERRETM ) ! ! Evaluate PHI and the derivative DPHI ! TEMP = Z( N ) / DELTA( N ) PHI = Z( N )*TEMP DPHI = TEMP*TEMP ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV + & ABS( TAU )*( DPSI+DPHI ) ! W = RHOINV + PHI + PSI 90 CONTINUE ! ! Return with INFO = 1, NITER = MAXIT and not converged ! INFO = 1 DLAM = D( I ) + TAU GO TO 250 ! ! End for the case I = N ! ELSE ! ! The case for I < N ! NITER = 1 IP1 = I + 1 ! ! Calculate initial guess ! DEL = D( IP1 ) - D( I ) MIDPT = DEL / TWO DO 100 J = 1, N DELTA( J ) = ( D( J )-D( I ) ) - MIDPT 100 CONTINUE ! PSI = ZERO DO 110 J = 1, I - 1 PSI = PSI + Z( J )*Z( J ) / DELTA( J ) 110 CONTINUE ! PHI = ZERO DO 120 J = N, I + 2, -1 PHI = PHI + Z( J )*Z( J ) / DELTA( J ) 120 CONTINUE C = RHOINV + PSI + PHI W = C + Z( I )*Z( I ) / DELTA( I ) + & Z( IP1 )*Z( IP1 ) / DELTA( IP1 ) ! IF( W.GT.ZERO ) THEN ! ! d(i)< the ith eigenvalue < (d(i)+d(i+1))/2 ! ! We choose d(i) as origin. ! ORGATI = .TRUE. A = C*DEL + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 ) B = Z( I )*Z( I )*DEL IF( A.GT.ZERO ) THEN TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) ELSE TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) END IF DLTLB = ZERO DLTUB = MIDPT ELSE ! ! (d(i)+d(i+1))/2 <= the ith eigenvalue < d(i+1) ! ! We choose d(i+1) as origin. ! ORGATI = .FALSE. A = C*DEL - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 ) B = Z( IP1 )*Z( IP1 )*DEL IF( A.LT.ZERO ) THEN TAU = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) ) ELSE TAU = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C ) END IF DLTLB = -MIDPT DLTUB = ZERO END IF ! IF( ORGATI ) THEN DO 130 J = 1, N DELTA( J ) = ( D( J )-D( I ) ) - TAU 130 CONTINUE ELSE DO 140 J = 1, N DELTA( J ) = ( D( J )-D( IP1 ) ) - TAU 140 CONTINUE END IF IF( ORGATI ) THEN II = I ELSE II = I + 1 END IF IIM1 = II - 1 IIP1 = II + 1 ! ! Evaluate PSI and the derivative DPSI ! DPSI = ZERO PSI = ZERO ERRETM = ZERO DO 150 J = 1, IIM1 TEMP = Z( J ) / DELTA( J ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 150 CONTINUE ERRETM = ABS( ERRETM ) ! ! Evaluate PHI and the derivative DPHI ! DPHI = ZERO PHI = ZERO DO 160 J = N, IIP1, -1 TEMP = Z( J ) / DELTA( J ) PHI = PHI + Z( J )*TEMP DPHI = DPHI + TEMP*TEMP ERRETM = ERRETM + PHI 160 CONTINUE ! W = RHOINV + PHI + PSI ! ! W is the value of the secular function with ! its ii-th element removed. ! SWTCH3 = .FALSE. IF( ORGATI ) THEN IF( W.LT.ZERO ) & SWTCH3 = .TRUE. ELSE IF( W.GT.ZERO ) & SWTCH3 = .TRUE. END IF IF( II.EQ.1 .OR. II.EQ.N ) & SWTCH3 = .FALSE. ! TEMP = Z( II ) / DELTA( II ) DW = DPSI + DPHI + TEMP*TEMP TEMP = Z( II )*TEMP W = W + TEMP ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + & THREE*ABS( TEMP ) + ABS( TAU )*DW ! ! Test for convergence ! IF( ABS( W ).LE.EPS*ERRETM ) THEN IF( ORGATI ) THEN DLAM = D( I ) + TAU ELSE DLAM = D( IP1 ) + TAU END IF GO TO 250 END IF ! IF( W.LE.ZERO ) THEN DLTLB = MAX( DLTLB, TAU ) ELSE DLTUB = MIN( DLTUB, TAU ) END IF ! ! Calculate the new step ! NITER = NITER + 1 IF( .NOT.SWTCH3 ) THEN IF( ORGATI ) THEN C = W - DELTA( IP1 )*DW - ( D( I )-D( IP1 ) )* & ( Z( I ) / DELTA( I ) )**2 ELSE C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )* & ( Z( IP1 ) / DELTA( IP1 ) )**2 END IF A = ( DELTA( I )+DELTA( IP1 ) )*W - & DELTA( I )*DELTA( IP1 )*DW B = DELTA( I )*DELTA( IP1 )*W IF( C.EQ.ZERO ) THEN IF( A.EQ.ZERO ) THEN IF( ORGATI ) THEN A = Z( I )*Z( I ) + DELTA( IP1 )*DELTA( IP1 )* & ( DPSI+DPHI ) ELSE A = Z( IP1 )*Z( IP1 ) + DELTA( I )*DELTA( I )* & ( DPSI+DPHI ) END IF END IF 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 ELSE ! ! Interpolation using THREE most relevant poles ! TEMP = RHOINV + PSI + PHI IF( ORGATI ) THEN TEMP1 = Z( IIM1 ) / DELTA( IIM1 ) TEMP1 = TEMP1*TEMP1 C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) - & ( D( IIM1 )-D( IIP1 ) )*TEMP1 ZZ( 1 ) = Z( IIM1 )*Z( IIM1 ) ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )* & ( ( DPSI-TEMP1 )+DPHI ) ELSE TEMP1 = Z( IIP1 ) / DELTA( IIP1 ) TEMP1 = TEMP1*TEMP1 C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) - & ( D( IIP1 )-D( IIM1 ) )*TEMP1 ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )* & ( DPSI+( DPHI-TEMP1 ) ) ZZ( 3 ) = Z( IIP1 )*Z( IIP1 ) END IF ZZ( 2 ) = Z( II )*Z( II ) CALL DLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA, & INFO ) IF( INFO.NE.0 ) & GO TO 250 END IF ! ! Note, eta should be positive if w is negative, and ! eta should be negative otherwise. However, ! if for some reason caused by roundoff, eta*w > 0, ! we simply use one Newton step instead. This way ! will guarantee eta*w < 0. ! IF( W*ETA.GE.ZERO ) & ETA = -W / DW TEMP = TAU + ETA IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN IF( W.LT.ZERO ) THEN ETA = ( DLTUB-TAU ) / TWO ELSE ETA = ( DLTLB-TAU ) / TWO END IF END IF ! PREW = W ! DO 180 J = 1, N DELTA( J ) = DELTA( J ) - ETA 180 CONTINUE ! ! Evaluate PSI and the derivative DPSI ! DPSI = ZERO PSI = ZERO ERRETM = ZERO DO 190 J = 1, IIM1 TEMP = Z( J ) / DELTA( J ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 190 CONTINUE ERRETM = ABS( ERRETM ) ! ! Evaluate PHI and the derivative DPHI ! DPHI = ZERO PHI = ZERO DO 200 J = N, IIP1, -1 TEMP = Z( J ) / DELTA( J ) PHI = PHI + Z( J )*TEMP DPHI = DPHI + TEMP*TEMP ERRETM = ERRETM + PHI 200 CONTINUE ! TEMP = Z( II ) / DELTA( II ) DW = DPSI + DPHI + TEMP*TEMP TEMP = Z( II )*TEMP W = RHOINV + PHI + PSI + TEMP ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + & THREE*ABS( TEMP ) + ABS( TAU+ETA )*DW ! SWTCH = .FALSE. IF( ORGATI ) THEN IF( -W.GT.ABS( PREW ) / TEN ) & SWTCH = .TRUE. ELSE IF( W.GT.ABS( PREW ) / TEN ) & SWTCH = .TRUE. END IF ! TAU = TAU + ETA ! ! Main loop to update the values of the array DELTA ! ITER = NITER + 1 ! DO 240 NITER = ITER, MAXIT ! ! Test for convergence ! IF( ABS( W ).LE.EPS*ERRETM ) THEN IF( ORGATI ) THEN DLAM = D( I ) + TAU ELSE DLAM = D( IP1 ) + TAU END IF GO TO 250 END IF ! IF( W.LE.ZERO ) THEN DLTLB = MAX( DLTLB, TAU ) ELSE DLTUB = MIN( DLTUB, TAU ) END IF ! ! Calculate the new step ! IF( .NOT.SWTCH3 ) THEN IF( .NOT.SWTCH ) THEN IF( ORGATI ) THEN C = W - DELTA( IP1 )*DW - & ( D( I )-D( IP1 ) )*( Z( I ) / DELTA( I ) )**2 ELSE C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )* & ( Z( IP1 ) / DELTA( IP1 ) )**2 END IF ELSE TEMP = Z( II ) / DELTA( II ) IF( ORGATI ) THEN DPSI = DPSI + TEMP*TEMP ELSE DPHI = DPHI + TEMP*TEMP END IF C = W - DELTA( I )*DPSI - DELTA( IP1 )*DPHI END IF A = ( DELTA( I )+DELTA( IP1 ) )*W - & DELTA( I )*DELTA( IP1 )*DW B = DELTA( I )*DELTA( IP1 )*W IF( C.EQ.ZERO ) THEN IF( A.EQ.ZERO ) THEN IF( .NOT.SWTCH ) THEN IF( ORGATI ) THEN A = Z( I )*Z( I ) + DELTA( IP1 )* & DELTA( IP1 )*( DPSI+DPHI ) ELSE A = Z( IP1 )*Z( IP1 ) + & DELTA( I )*DELTA( I )*( DPSI+DPHI ) END IF ELSE A = DELTA( I )*DELTA( I )*DPSI + & DELTA( IP1 )*DELTA( IP1 )*DPHI END IF END IF 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 ELSE ! ! Interpolation using THREE most relevant poles ! TEMP = RHOINV + PSI + PHI IF( SWTCH ) THEN C = TEMP - DELTA( IIM1 )*DPSI - DELTA( IIP1 )*DPHI ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )*DPSI ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )*DPHI ELSE IF( ORGATI ) THEN TEMP1 = Z( IIM1 ) / DELTA( IIM1 ) TEMP1 = TEMP1*TEMP1 C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) - & ( D( IIM1 )-D( IIP1 ) )*TEMP1 ZZ( 1 ) = Z( IIM1 )*Z( IIM1 ) ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )* & ( ( DPSI-TEMP1 )+DPHI ) ELSE TEMP1 = Z( IIP1 ) / DELTA( IIP1 ) TEMP1 = TEMP1*TEMP1 C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) - & ( D( IIP1 )-D( IIM1 ) )*TEMP1 ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )* & ( DPSI+( DPHI-TEMP1 ) ) ZZ( 3 ) = Z( IIP1 )*Z( IIP1 ) END IF END IF CALL DLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA, & INFO ) IF( INFO.NE.0 ) & GO TO 250 END IF ! ! Note, eta should be positive if w is negative, and ! eta should be negative otherwise. However, ! if for some reason caused by roundoff, eta*w > 0, ! we simply use one Newton step instead. This way ! will guarantee eta*w < 0. ! IF( W*ETA.GE.ZERO ) & ETA = -W / DW TEMP = TAU + ETA IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN IF( W.LT.ZERO ) THEN ETA = ( DLTUB-TAU ) / TWO ELSE ETA = ( DLTLB-TAU ) / TWO END IF END IF ! DO 210 J = 1, N DELTA( J ) = DELTA( J ) - ETA 210 CONTINUE ! TAU = TAU + ETA PREW = W ! ! Evaluate PSI and the derivative DPSI ! DPSI = ZERO PSI = ZERO ERRETM = ZERO DO 220 J = 1, IIM1 TEMP = Z( J ) / DELTA( J ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 220 CONTINUE ERRETM = ABS( ERRETM ) ! ! Evaluate PHI and the derivative DPHI ! DPHI = ZERO PHI = ZERO DO 230 J = N, IIP1, -1 TEMP = Z( J ) / DELTA( J ) PHI = PHI + Z( J )*TEMP DPHI = DPHI + TEMP*TEMP ERRETM = ERRETM + PHI 230 CONTINUE ! TEMP = Z( II ) / DELTA( II ) DW = DPSI + DPHI + TEMP*TEMP TEMP = Z( II )*TEMP W = RHOINV + PHI + PSI + TEMP ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + & THREE*ABS( TEMP ) + ABS( TAU )*DW IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN ) & SWTCH = .NOT.SWTCH ! 240 CONTINUE ! ! Return with INFO = 1, NITER = MAXIT and not converged ! INFO = 1 IF( ORGATI ) THEN DLAM = D( I ) + TAU ELSE DLAM = D( IP1 ) + TAU END IF ! END IF ! 250 CONTINUE ! RETURN ! ! End of DLAED4 ! END