*> \brief \b SLAED4 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 SLAED4 + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * SUBROUTINE SLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO ) * * .. Scalar Arguments .. * INTEGER I, INFO, N * REAL DLAM, RHO * .. * .. Array Arguments .. * REAL 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 REAL 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 REAL array, dimension (N) *> The components of the updating vector. *> \endverbatim *> *> \param[out] DELTA *> \verbatim *> DELTA is REAL array, dimension (N) *> If N > 2, DELTA contains (D(j) - lambda_I) in its j-th *> component. If N = 1, then DELTA(1) = 1. If N = 2, see SLAED5 *> for detail. The vector DELTA contains the information necessary *> to construct the eigenvectors by SLAED3 and SLAED9. *> \endverbatim *> *> \param[in] RHO *> \verbatim *> RHO is REAL *> The scalar in the symmetric updating formula. *> \endverbatim *> *> \param[out] DLAM *> \verbatim *> DLAM is REAL *> 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. * *> \ingroup auxOTHERcomputational * *> \par Contributors: * ================== *> *> Ren-Cang Li, Computer Science Division, University of California *> at Berkeley, USA *> * ===================================================================== SUBROUTINE SLAED4( N, I, D, Z, DELTA, RHO, DLAM, 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 .. INTEGER I, INFO, N REAL DLAM, RHO * .. * .. Array Arguments .. REAL D( * ), DELTA( * ), Z( * ) * .. * * ===================================================================== * * .. Parameters .. INTEGER MAXIT PARAMETER ( MAXIT = 30 ) REAL ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, $ THREE = 3.0E0, FOUR = 4.0E0, EIGHT = 8.0E0, $ TEN = 10.0E0 ) * .. * .. Local Scalars .. LOGICAL ORGATI, SWTCH, SWTCH3 INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER REAL A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW, $ EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI, $ RHOINV, TAU, TEMP, TEMP1, W * .. * .. Local Arrays .. REAL ZZ( 3 ) * .. * .. External Functions .. REAL SLAMCH EXTERNAL SLAMCH * .. * .. External Subroutines .. EXTERNAL SLAED5, SLAED6 * .. * .. 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 SLAED5( I, D, Z, DELTA, RHO, DLAM ) RETURN END IF * * Compute machine epsilon * EPS = SLAMCH( '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 SLAED6( 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 SLAED6( 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 SLAED4 * END