forked from lijiext/lammps
918 lines
26 KiB
Fortran
918 lines
26 KiB
Fortran
*> \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
|