[BACK]Return to dlasd4.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dlasd4.f between versions 1.4 and 1.18

version 1.4, 2010/08/06 15:32:29 version 1.18, 2017/06/17 10:53:57
Line 1 Line 1
   *> \brief \b DLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modification to a positive diagonal matrix. Used by dbdsdc.
   *
   *  =========== DOCUMENTATION ===========
   *
   * Online html documentation available at
   *            http://www.netlib.org/lapack/explore-html/
   *
   *> \htmlonly
   *> Download DLASD4 + dependencies
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd4.f">
   *> [TGZ]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd4.f">
   *> [ZIP]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd4.f">
   *> [TXT]</a>
   *> \endhtmlonly
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE DLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )
   *
   *       .. Scalar Arguments ..
   *       INTEGER            I, INFO, N
   *       DOUBLE PRECISION   RHO, SIGMA
   *       ..
   *       .. Array Arguments ..
   *       DOUBLE PRECISION   D( * ), DELTA( * ), WORK( * ), Z( * )
   *       ..
   *
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> This subroutine computes the square root of the I-th updated
   *> eigenvalue of a positive symmetric rank-one modification to
   *> a positive diagonal matrix whose entries are given as the squares
   *> of the corresponding entries in the array d, and that
   *>
   *>        0 <= 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 ) * 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, 0 <= 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 .ne. 1, DELTA contains (D(j) - sigma_I) in its  j-th
   *>         component.  If N = 1, then DELTA(1) = 1.  The vector DELTA
   *>         contains the information necessary to construct the
   *>         (singular) eigenvectors.
   *> \endverbatim
   *>
   *> \param[in] RHO
   *> \verbatim
   *>          RHO is DOUBLE PRECISION
   *>         The scalar in the symmetric updating formula.
   *> \endverbatim
   *>
   *> \param[out] SIGMA
   *> \verbatim
   *>          SIGMA is DOUBLE PRECISION
   *>         The computed sigma_I, the I-th updated eigenvalue.
   *> \endverbatim
   *>
   *> \param[out] WORK
   *> \verbatim
   *>          WORK is DOUBLE PRECISION array, dimension ( N )
   *>         If N .ne. 1, WORK contains (D(j) + sigma_I) in its  j-th
   *>         component.  If N = 1, then WORK( 1 ) = 1.
   *> \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 OTHERauxiliary
   *
   *> \par Contributors:
   *  ==================
   *>
   *>     Ren-Cang Li, Computer Science Division, University of California
   *>     at Berkeley, USA
   *>
   *  =====================================================================
       SUBROUTINE DLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )        SUBROUTINE DLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )
 *  *
 *  -- LAPACK auxiliary routine (version 3.2) --  *  -- LAPACK auxiliary routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 *     November 2006  *     December 2016
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       INTEGER            I, INFO, N        INTEGER            I, INFO, N
Line 13 Line 166
       DOUBLE PRECISION   D( * ), DELTA( * ), WORK( * ), Z( * )        DOUBLE PRECISION   D( * ), DELTA( * ), WORK( * ), Z( * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  This subroutine computes the square root of the I-th updated  
 *  eigenvalue of a positive symmetric rank-one modification to  
 *  a positive diagonal matrix whose entries are given as the squares  
 *  of the corresponding entries in the array d, and that  
 *  
 *         0 <= 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 ) * 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.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  N      (input) INTEGER  
 *         The length of all arrays.  
 *  
 *  I      (input) INTEGER  
 *         The index of the eigenvalue to be computed.  1 <= I <= N.  
 *  
 *  D      (input) DOUBLE PRECISION array, dimension ( N )  
 *         The original eigenvalues.  It is assumed that they are in  
 *         order, 0 <= D(I) < D(J)  for I < J.  
 *  
 *  Z      (input) DOUBLE PRECISION array, dimension ( N )  
 *         The components of the updating vector.  
 *  
 *  DELTA  (output) DOUBLE PRECISION array, dimension ( N )  
 *         If N .ne. 1, DELTA contains (D(j) - sigma_I) in its  j-th  
 *         component.  If N = 1, then DELTA(1) = 1.  The vector DELTA  
 *         contains the information necessary to construct the  
 *         (singular) eigenvectors.  
 *  
 *  RHO    (input) DOUBLE PRECISION  
 *         The scalar in the symmetric updating formula.  
 *  
 *  SIGMA  (output) DOUBLE PRECISION  
 *         The computed sigma_I, the I-th updated eigenvalue.  
 *  
 *  WORK   (workspace) DOUBLE PRECISION array, dimension ( N )  
 *         If N .ne. 1, WORK contains (D(j) + sigma_I) in its  j-th  
 *         component.  If N = 1, then WORK( 1 ) = 1.  
 *  
 *  INFO   (output) INTEGER  
 *         = 0:  successful exit  
 *         > 0:  if INFO = 1, the updating process failed.  
 *  
 *  Internal Parameters  
 *  ===================  
 *  
 *  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.  
 *  
 *  Further Details  
 *  ===============  
 *  
 *  Based on contributions by  
 *     Ren-Cang Li, Computer Science Division, University of California  
 *     at Berkeley, USA  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
       INTEGER            MAXIT        INTEGER            MAXIT
       PARAMETER          ( MAXIT = 20 )        PARAMETER          ( MAXIT = 400 )
       DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN        DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,        PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
      $                   THREE = 3.0D+0, FOUR = 4.0D+0, EIGHT = 8.0D+0,       $                   THREE = 3.0D+0, FOUR = 4.0D+0, EIGHT = 8.0D+0,
      $                   TEN = 10.0D+0 )       $                   TEN = 10.0D+0 )
 *     ..  *     ..
 *     .. Local Scalars ..  *     .. Local Scalars ..
       LOGICAL            ORGATI, SWTCH, SWTCH3        LOGICAL            ORGATI, SWTCH, SWTCH3, GEOMAVG
       INTEGER            II, IIM1, IIP1, IP1, ITER, J, NITER        INTEGER            II, IIM1, IIP1, IP1, ITER, J, NITER
       DOUBLE PRECISION   A, B, C, DELSQ, DELSQ2, DPHI, DPSI, DTIIM,        DOUBLE PRECISION   A, B, C, DELSQ, DELSQ2, SQ2, DPHI, DPSI, DTIIM,
      $                   DTIIP, DTIPSQ, DTISQ, DTNSQ, DTNSQ1, DW, EPS,       $                   DTIIP, DTIPSQ, DTISQ, DTNSQ, DTNSQ1, DW, EPS,
      $                   ERRETM, ETA, PHI, PREW, PSI, RHOINV, SG2LB,       $                   ERRETM, ETA, PHI, PREW, PSI, RHOINV, SGLB,
      $                   SG2UB, TAU, TEMP, TEMP1, TEMP2, W       $                   SGUB, TAU, TAU2, TEMP, TEMP1, TEMP2, W
 *     ..  *     ..
 *     .. Local Arrays ..  *     .. Local Arrays ..
       DOUBLE PRECISION   DD( 3 ), ZZ( 3 )        DOUBLE PRECISION   DD( 3 ), ZZ( 3 )
Line 148 Line 223
 *  *
       EPS = DLAMCH( 'Epsilon' )        EPS = DLAMCH( 'Epsilon' )
       RHOINV = ONE / RHO        RHOINV = ONE / RHO
         TAU2= ZERO
 *  *
 *     The case I = N  *     The case I = N
 *  *
Line 186 Line 262
      $             ( D( N )-D( N-1 )+RHO / ( D( N )+TEMP1 ) ) ) +       $             ( D( N )-D( N-1 )+RHO / ( D( N )+TEMP1 ) ) ) +
      $             Z( N )*Z( N ) / RHO       $             Z( N )*Z( N ) / RHO
 *  *
 *           The following TAU is to approximate  *           The following TAU2 is to approximate
 *           SIGMA_n^2 - D( N )*D( N )  *           SIGMA_n^2 - D( N )*D( N )
 *  *
             IF( C.LE.TEMP ) THEN              IF( C.LE.TEMP ) THEN
Line 196 Line 272
                A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )                 A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
                B = Z( N )*Z( N )*DELSQ                 B = Z( N )*Z( N )*DELSQ
                IF( A.LT.ZERO ) THEN                 IF( A.LT.ZERO ) THEN
                   TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )                    TAU2 = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
                ELSE                 ELSE
                   TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )                    TAU2 = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
                END IF                 END IF
                  TAU = TAU2 / ( D( N )+SQRT( D( N )*D( N )+TAU2 ) )
             END IF              END IF
 *  *
 *           It can be proved that  *           It can be proved that
 *               D(N)^2+RHO/2 <= SIGMA_n^2 < D(N)^2+TAU <= D(N)^2+RHO  *               D(N)^2+RHO/2 <= SIGMA_n^2 < D(N)^2+TAU2 <= D(N)^2+RHO
 *  *
          ELSE           ELSE
             DELSQ = ( D( N )-D( N-1 ) )*( D( N )+D( N-1 ) )              DELSQ = ( D( N )-D( N-1 ) )*( D( N )+D( N-1 ) )
             A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )              A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
             B = Z( N )*Z( N )*DELSQ              B = Z( N )*Z( N )*DELSQ
 *  *
 *           The following TAU is to approximate  *           The following TAU2 is to approximate
 *           SIGMA_n^2 - D( N )*D( N )  *           SIGMA_n^2 - D( N )*D( N )
 *  *
             IF( A.LT.ZERO ) THEN              IF( A.LT.ZERO ) THEN
                TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )                 TAU2 = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
             ELSE              ELSE
                TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )                 TAU2 = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
             END IF              END IF
               TAU = TAU2 / ( D( N )+SQRT( D( N )*D( N )+TAU2 ) )
   
 *  *
 *           It can be proved that  *           It can be proved that
 *           D(N)^2 < D(N)^2+TAU < SIGMA(N)^2 < D(N)^2+RHO/2  *           D(N)^2 < D(N)^2+TAU2 < SIGMA(N)^2 < D(N)^2+RHO/2
 *  *
          END IF           END IF
 *  *
 *        The following ETA is to approximate SIGMA_n - D( N )  *        The following TAU is to approximate SIGMA_n - D( N )
 *  *
          ETA = TAU / ( D( N )+SQRT( D( N )*D( N )+TAU ) )  *         TAU = TAU2 / ( D( N )+SQRT( D( N )*D( N )+TAU2 ) )
 *  *
          SIGMA = D( N ) + ETA           SIGMA = D( N ) + TAU
          DO 30 J = 1, N           DO 30 J = 1, N
             DELTA( J ) = ( D( J )-D( I ) ) - ETA              DELTA( J ) = ( D( J )-D( N ) ) - TAU
             WORK( J ) = D( J ) + D( I ) + ETA              WORK( J ) = D( J ) + D( N ) + TAU
    30    CONTINUE     30    CONTINUE
 *  *
 *        Evaluate PSI and the derivative DPSI  *        Evaluate PSI and the derivative DPSI
Line 252 Line 331
          TEMP = Z( N ) / ( DELTA( N )*WORK( N ) )           TEMP = Z( N ) / ( DELTA( N )*WORK( N ) )
          PHI = Z( N )*TEMP           PHI = Z( N )*TEMP
          DPHI = TEMP*TEMP           DPHI = TEMP*TEMP
          ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +           ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV
      $            ABS( TAU )*( DPSI+DPHI )  *    $          + ABS( TAU2 )*( DPSI+DPHI )
 *  *
          W = RHOINV + PHI + PSI           W = RHOINV + PHI + PSI
 *  *
Line 293 Line 372
          IF( TEMP.GT.RHO )           IF( TEMP.GT.RHO )
      $      ETA = RHO + DTNSQ       $      ETA = RHO + DTNSQ
 *  *
          TAU = TAU + ETA  
          ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )           ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
            TAU = TAU + ETA
            SIGMA = SIGMA + ETA
   *
          DO 50 J = 1, N           DO 50 J = 1, N
             DELTA( J ) = DELTA( J ) - ETA              DELTA( J ) = DELTA( J ) - ETA
             WORK( J ) = WORK( J ) + ETA              WORK( J ) = WORK( J ) + ETA
    50    CONTINUE     50    CONTINUE
 *  *
          SIGMA = SIGMA + ETA  
 *  
 *        Evaluate PSI and the derivative DPSI  *        Evaluate PSI and the derivative DPSI
 *  *
          DPSI = ZERO           DPSI = ZERO
Line 317 Line 396
 *  *
 *        Evaluate PHI and the derivative DPHI  *        Evaluate PHI and the derivative DPHI
 *  *
          TEMP = Z( N ) / ( WORK( N )*DELTA( N ) )           TAU2 = WORK( N )*DELTA( N )
            TEMP = Z( N ) / TAU2
          PHI = Z( N )*TEMP           PHI = Z( N )*TEMP
          DPHI = TEMP*TEMP           DPHI = TEMP*TEMP
          ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +           ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV
      $            ABS( TAU )*( DPSI+DPHI )  *    $          + ABS( TAU2 )*( DPSI+DPHI )
 *  *
          W = RHOINV + PHI + PSI           W = RHOINV + PHI + PSI
 *  *
Line 362 Line 442
             IF( TEMP.LE.ZERO )              IF( TEMP.LE.ZERO )
      $         ETA = ETA / TWO       $         ETA = ETA / TWO
 *  *
             TAU = TAU + ETA  
             ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )              ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
               TAU = TAU + ETA
               SIGMA = SIGMA + ETA
   *
             DO 70 J = 1, N              DO 70 J = 1, N
                DELTA( J ) = DELTA( J ) - ETA                 DELTA( J ) = DELTA( J ) - ETA
                WORK( J ) = WORK( J ) + ETA                 WORK( J ) = WORK( J ) + ETA
    70       CONTINUE     70       CONTINUE
 *  *
             SIGMA = SIGMA + ETA  
 *  
 *           Evaluate PSI and the derivative DPSI  *           Evaluate PSI and the derivative DPSI
 *  *
             DPSI = ZERO              DPSI = ZERO
Line 386 Line 466
 *  *
 *           Evaluate PHI and the derivative DPHI  *           Evaluate PHI and the derivative DPHI
 *  *
             TEMP = Z( N ) / ( WORK( N )*DELTA( N ) )              TAU2 = WORK( N )*DELTA( N )
               TEMP = Z( N ) / TAU2
             PHI = Z( N )*TEMP              PHI = Z( N )*TEMP
             DPHI = TEMP*TEMP              DPHI = TEMP*TEMP
             ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +              ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV
      $               ABS( TAU )*( DPSI+DPHI )  *    $             + ABS( TAU2 )*( DPSI+DPHI )
 *  *
             W = RHOINV + PHI + PSI              W = RHOINV + PHI + PSI
    90    CONTINUE     90    CONTINUE
Line 413 Line 494
 *  *
          DELSQ = ( D( IP1 )-D( I ) )*( D( IP1 )+D( I ) )           DELSQ = ( D( IP1 )-D( I ) )*( D( IP1 )+D( I ) )
          DELSQ2 = DELSQ / TWO           DELSQ2 = DELSQ / TWO
          TEMP = DELSQ2 / ( D( I )+SQRT( D( I )*D( I )+DELSQ2 ) )           SQ2=SQRT( ( D( I )*D( I )+D( IP1 )*D( IP1 ) ) / TWO )
            TEMP = DELSQ2 / ( D( I )+SQ2 )
          DO 100 J = 1, N           DO 100 J = 1, N
             WORK( J ) = D( J ) + D( I ) + TEMP              WORK( J ) = D( J ) + D( I ) + TEMP
             DELTA( J ) = ( D( J )-D( I ) ) - TEMP              DELTA( J ) = ( D( J )-D( I ) ) - TEMP
Line 432 Line 514
          W = C + Z( I )*Z( I ) / ( WORK( I )*DELTA( I ) ) +           W = C + Z( I )*Z( I ) / ( WORK( I )*DELTA( I ) ) +
      $       Z( IP1 )*Z( IP1 ) / ( WORK( IP1 )*DELTA( IP1 ) )       $       Z( IP1 )*Z( IP1 ) / ( WORK( IP1 )*DELTA( IP1 ) )
 *  *
            GEOMAVG = .FALSE.
          IF( W.GT.ZERO ) THEN           IF( W.GT.ZERO ) THEN
 *  *
 *           d(i)^2 < the ith sigma^2 < (d(i)^2+d(i+1)^2)/2  *           d(i)^2 < the ith sigma^2 < (d(i)^2+d(i+1)^2)/2
Line 439 Line 522
 *           We choose d(i) as origin.  *           We choose d(i) as origin.
 *  *
             ORGATI = .TRUE.              ORGATI = .TRUE.
             SG2LB = ZERO              II = I
             SG2UB = DELSQ2              SGLB = ZERO
               SGUB = DELSQ2  / ( D( I )+SQ2 )
             A = C*DELSQ + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 )              A = C*DELSQ + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 )
             B = Z( I )*Z( I )*DELSQ              B = Z( I )*Z( I )*DELSQ
             IF( A.GT.ZERO ) THEN              IF( A.GT.ZERO ) THEN
                TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )                 TAU2 = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
             ELSE              ELSE
                TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )                 TAU2 = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
             END IF              END IF
 *  *
 *           TAU now is an estimation of SIGMA^2 - D( I )^2. The  *           TAU2 now is an estimation of SIGMA^2 - D( I )^2. The
 *           following, however, is the corresponding estimation of  *           following, however, is the corresponding estimation of
 *           SIGMA - D( I ).  *           SIGMA - D( I ).
 *  *
             ETA = TAU / ( D( I )+SQRT( D( I )*D( I )+TAU ) )              TAU = TAU2 / ( D( I )+SQRT( D( I )*D( I )+TAU2 ) )
               TEMP = SQRT(EPS)
               IF( (D(I).LE.TEMP*D(IP1)).AND.(ABS(Z(I)).LE.TEMP)
        $                               .AND.(D(I).GT.ZERO) ) THEN
                  TAU = MIN( TEN*D(I), SGUB )
                  GEOMAVG = .TRUE.
               END IF
          ELSE           ELSE
 *  *
 *           (d(i)^2+d(i+1)^2)/2 <= the ith sigma^2 < d(i+1)^2/2  *           (d(i)^2+d(i+1)^2)/2 <= the ith sigma^2 < d(i+1)^2/2
Line 461 Line 551
 *           We choose d(i+1) as origin.  *           We choose d(i+1) as origin.
 *  *
             ORGATI = .FALSE.              ORGATI = .FALSE.
             SG2LB = -DELSQ2              II = IP1
             SG2UB = ZERO              SGLB = -DELSQ2  / ( D( II )+SQ2 )
               SGUB = ZERO
             A = C*DELSQ - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 )              A = C*DELSQ - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 )
             B = Z( IP1 )*Z( IP1 )*DELSQ              B = Z( IP1 )*Z( IP1 )*DELSQ
             IF( A.LT.ZERO ) THEN              IF( A.LT.ZERO ) THEN
                TAU = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) )                 TAU2 = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) )
             ELSE              ELSE
                TAU = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C )                 TAU2 = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C )
             END IF              END IF
 *  *
 *           TAU now is an estimation of SIGMA^2 - D( IP1 )^2. The  *           TAU2 now is an estimation of SIGMA^2 - D( IP1 )^2. The
 *           following, however, is the corresponding estimation of  *           following, however, is the corresponding estimation of
 *           SIGMA - D( IP1 ).  *           SIGMA - D( IP1 ).
 *  *
             ETA = TAU / ( D( IP1 )+SQRT( ABS( D( IP1 )*D( IP1 )+              TAU = TAU2 / ( D( IP1 )+SQRT( ABS( D( IP1 )*D( IP1 )+
      $            TAU ) ) )       $            TAU2 ) ) )
          END IF           END IF
 *  *
          IF( ORGATI ) THEN           SIGMA = D( II ) + TAU
             II = I           DO 130 J = 1, N
             SIGMA = D( I ) + ETA              WORK( J ) = D( J ) + D( II ) + TAU
             DO 130 J = 1, N              DELTA( J ) = ( D( J )-D( II ) ) - TAU
                WORK( J ) = D( J ) + D( I ) + ETA    130    CONTINUE
                DELTA( J ) = ( D( J )-D( I ) ) - ETA  
   130       CONTINUE  
          ELSE  
             II = I + 1  
             SIGMA = D( IP1 ) + ETA  
             DO 140 J = 1, N  
                WORK( J ) = D( J ) + D( IP1 ) + ETA  
                DELTA( J ) = ( D( J )-D( IP1 ) ) - ETA  
   140       CONTINUE  
          END IF  
          IIM1 = II - 1           IIM1 = II - 1
          IIP1 = II + 1           IIP1 = II + 1
 *  *
Line 541 Line 622
          DW = DPSI + DPHI + TEMP*TEMP           DW = DPSI + DPHI + TEMP*TEMP
          TEMP = Z( II )*TEMP           TEMP = Z( II )*TEMP
          W = W + TEMP           W = W + TEMP
          ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +           ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV
      $            THREE*ABS( TEMP ) + ABS( TAU )*DW       $          + THREE*ABS( TEMP )
   *    $          + ABS( TAU2 )*DW
 *  *
 *        Test for convergence  *        Test for convergence
 *  *
Line 551 Line 633
          END IF           END IF
 *  *
          IF( W.LE.ZERO ) THEN           IF( W.LE.ZERO ) THEN
             SG2LB = MAX( SG2LB, TAU )              SGLB = MAX( SGLB, TAU )
          ELSE           ELSE
             SG2UB = MIN( SG2UB, TAU )              SGUB = MIN( SGUB, TAU )
          END IF           END IF
 *  *
 *        Calculate the new step  *        Calculate the new step
Line 618 Line 700
             DD( 2 ) = DELTA( II )*WORK( II )              DD( 2 ) = DELTA( II )*WORK( II )
             DD( 3 ) = DTIIP              DD( 3 ) = DTIIP
             CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )              CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )
             IF( INFO.NE.0 )  *
      $         GO TO 240              IF( INFO.NE.0 ) THEN
   *
   *              If INFO is not 0, i.e., DLAED6 failed, switch back
   *              to 2 pole interpolation.
   *
                  SWTCH3 = .FALSE.
                  INFO = 0
                  DTIPSQ = WORK( IP1 )*DELTA( IP1 )
                  DTISQ = WORK( I )*DELTA( I )
                  IF( ORGATI ) THEN
                     C = W - DTIPSQ*DW + DELSQ*( Z( I ) / DTISQ )**2
                  ELSE
                     C = W - DTISQ*DW - DELSQ*( Z( IP1 ) / DTIPSQ )**2
                  END IF
                  A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW
                  B = DTIPSQ*DTISQ*W
                  IF( C.EQ.ZERO ) THEN
                     IF( A.EQ.ZERO ) THEN
                        IF( ORGATI ) THEN
                           A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ*( DPSI+DPHI )
                        ELSE
                           A = Z( IP1 )*Z( IP1 ) + DTISQ*DTISQ*( 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
               END IF
          END IF           END IF
 *  *
 *        Note, eta should be positive if w is negative, and  *        Note, eta should be positive if w is negative, and
Line 630 Line 742
 *  *
          IF( W*ETA.GE.ZERO )           IF( W*ETA.GE.ZERO )
      $      ETA = -W / DW       $      ETA = -W / DW
          IF( ORGATI ) THEN  *
             TEMP1 = WORK( I )*DELTA( I )           ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )
             TEMP = ETA - TEMP1           TEMP = TAU + ETA
          ELSE           IF( TEMP.GT.SGUB .OR. TEMP.LT.SGLB ) THEN
             TEMP1 = WORK( IP1 )*DELTA( IP1 )  
             TEMP = ETA - TEMP1  
          END IF  
          IF( TEMP.GT.SG2UB .OR. TEMP.LT.SG2LB ) THEN  
             IF( W.LT.ZERO ) THEN              IF( W.LT.ZERO ) THEN
                ETA = ( SG2UB-TAU ) / TWO                 ETA = ( SGUB-TAU ) / TWO
             ELSE              ELSE
                ETA = ( SG2LB-TAU ) / TWO                 ETA = ( SGLB-TAU ) / TWO
               END IF
               IF( GEOMAVG ) THEN
                  IF( W .LT. ZERO ) THEN
                     IF( TAU .GT. ZERO ) THEN
                        ETA = SQRT(SGUB*TAU)-TAU
                     END IF
                  ELSE
                     IF( SGLB .GT. ZERO ) THEN
                        ETA = SQRT(SGLB*TAU)-TAU
                     END IF
                  END IF
             END IF              END IF
          END IF           END IF
 *  *
          TAU = TAU + ETA  
          ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )  
 *  
          PREW = W           PREW = W
 *  *
            TAU = TAU + ETA
          SIGMA = SIGMA + ETA           SIGMA = SIGMA + ETA
   *
          DO 170 J = 1, N           DO 170 J = 1, N
             WORK( J ) = WORK( J ) + ETA              WORK( J ) = WORK( J ) + ETA
             DELTA( J ) = DELTA( J ) - ETA              DELTA( J ) = DELTA( J ) - ETA
Line 680 Line 798
             ERRETM = ERRETM + PHI              ERRETM = ERRETM + PHI
   190    CONTINUE    190    CONTINUE
 *  *
          TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )           TAU2 = WORK( II )*DELTA( II )
            TEMP = Z( II ) / TAU2
          DW = DPSI + DPHI + TEMP*TEMP           DW = DPSI + DPHI + TEMP*TEMP
          TEMP = Z( II )*TEMP           TEMP = Z( II )*TEMP
          W = RHOINV + PHI + PSI + TEMP           W = RHOINV + PHI + PSI + TEMP
          ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +           ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV
      $            THREE*ABS( TEMP ) + ABS( TAU )*DW       $          + THREE*ABS( TEMP )
 *  *    $          + ABS( TAU2 )*DW
          IF( W.LE.ZERO ) THEN  
             SG2LB = MAX( SG2LB, TAU )  
          ELSE  
             SG2UB = MIN( SG2UB, TAU )  
          END IF  
 *  *
          SWTCH = .FALSE.           SWTCH = .FALSE.
          IF( ORGATI ) THEN           IF( ORGATI ) THEN
Line 711 Line 825
 *           Test for convergence  *           Test for convergence
 *  *
             IF( ABS( W ).LE.EPS*ERRETM ) THEN              IF( ABS( W ).LE.EPS*ERRETM ) THEN
   *     $          .OR. (SGUB-SGLB).LE.EIGHT*ABS(SGUB+SGLB) ) THEN
                GO TO 240                 GO TO 240
             END IF              END IF
 *  *
               IF( W.LE.ZERO ) THEN
                  SGLB = MAX( SGLB, TAU )
               ELSE
                  SGUB = MIN( SGUB, TAU )
               END IF
   *
 *           Calculate the new step  *           Calculate the new step
 *  *
             IF( .NOT.SWTCH3 ) THEN              IF( .NOT.SWTCH3 ) THEN
Line 798 Line 919
                DD( 2 ) = DELTA( II )*WORK( II )                 DD( 2 ) = DELTA( II )*WORK( II )
                DD( 3 ) = DTIIP                 DD( 3 ) = DTIIP
                CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )                 CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )
                IF( INFO.NE.0 )  *
      $            GO TO 240                 IF( INFO.NE.0 ) THEN
   *
   *                 If INFO is not 0, i.e., DLAED6 failed, switch
   *                 back to two pole interpolation
   *
                     SWTCH3 = .FALSE.
                     INFO = 0
                     DTIPSQ = WORK( IP1 )*DELTA( IP1 )
                     DTISQ = WORK( I )*DELTA( I )
                     IF( .NOT.SWTCH ) THEN
                        IF( ORGATI ) THEN
                           C = W - DTIPSQ*DW + DELSQ*( Z( I )/DTISQ )**2
                        ELSE
                           C = W - DTISQ*DW - DELSQ*( Z( IP1 )/DTIPSQ )**2
                        END IF
                     ELSE
                        TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
                        IF( ORGATI ) THEN
                           DPSI = DPSI + TEMP*TEMP
                        ELSE
                           DPHI = DPHI + TEMP*TEMP
                        END IF
                        C = W - DTISQ*DPSI - DTIPSQ*DPHI
                     END IF
                     A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW
                     B = DTIPSQ*DTISQ*W
                     IF( C.EQ.ZERO ) THEN
                        IF( A.EQ.ZERO ) THEN
                           IF( .NOT.SWTCH ) THEN
                              IF( ORGATI ) THEN
                                 A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ*
        $                            ( DPSI+DPHI )
                              ELSE
                                 A = Z( IP1 )*Z( IP1 ) +
        $                            DTISQ*DTISQ*( DPSI+DPHI )
                              END IF
                           ELSE
                              A = DTISQ*DTISQ*DPSI + DTIPSQ*DTIPSQ*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
                  END IF
             END IF              END IF
 *  *
 *           Note, eta should be positive if w is negative, and  *           Note, eta should be positive if w is negative, and
Line 810 Line 977
 *  *
             IF( W*ETA.GE.ZERO )              IF( W*ETA.GE.ZERO )
      $         ETA = -W / DW       $         ETA = -W / DW
             IF( ORGATI ) THEN  *
                TEMP1 = WORK( I )*DELTA( I )              ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )
                TEMP = ETA - TEMP1              TEMP=TAU+ETA
             ELSE              IF( TEMP.GT.SGUB .OR. TEMP.LT.SGLB ) THEN
                TEMP1 = WORK( IP1 )*DELTA( IP1 )  
                TEMP = ETA - TEMP1  
             END IF  
             IF( TEMP.GT.SG2UB .OR. TEMP.LT.SG2LB ) THEN  
                IF( W.LT.ZERO ) THEN                 IF( W.LT.ZERO ) THEN
                   ETA = ( SG2UB-TAU ) / TWO                    ETA = ( SGUB-TAU ) / TWO
                ELSE                 ELSE
                   ETA = ( SG2LB-TAU ) / TWO                    ETA = ( SGLB-TAU ) / TWO
                  END IF
                  IF( GEOMAVG ) THEN
                     IF( W .LT. ZERO ) THEN
                        IF( TAU .GT. ZERO ) THEN
                           ETA = SQRT(SGUB*TAU)-TAU
                        END IF
                     ELSE
                        IF( SGLB .GT. ZERO ) THEN
                           ETA = SQRT(SGLB*TAU)-TAU
                        END IF
                     END IF
                END IF                 END IF
             END IF              END IF
 *  *
             TAU = TAU + ETA              PREW = W
             ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )  
 *  *
               TAU = TAU + ETA
             SIGMA = SIGMA + ETA              SIGMA = SIGMA + ETA
   *
             DO 200 J = 1, N              DO 200 J = 1, N
                WORK( J ) = WORK( J ) + ETA                 WORK( J ) = WORK( J ) + ETA
                DELTA( J ) = DELTA( J ) - ETA                 DELTA( J ) = DELTA( J ) - ETA
   200       CONTINUE    200       CONTINUE
 *  *
             PREW = W  
 *  
 *           Evaluate PSI and the derivative DPSI  *           Evaluate PSI and the derivative DPSI
 *  *
             DPSI = ZERO              DPSI = ZERO
Line 860 Line 1033
                ERRETM = ERRETM + PHI                 ERRETM = ERRETM + PHI
   220       CONTINUE    220       CONTINUE
 *  *
             TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )              TAU2 = WORK( II )*DELTA( II )
               TEMP = Z( II ) / TAU2
             DW = DPSI + DPHI + TEMP*TEMP              DW = DPSI + DPHI + TEMP*TEMP
             TEMP = Z( II )*TEMP              TEMP = Z( II )*TEMP
             W = RHOINV + PHI + PSI + TEMP              W = RHOINV + PHI + PSI + TEMP
             ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +              ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV
      $               THREE*ABS( TEMP ) + ABS( TAU )*DW       $             + THREE*ABS( TEMP )
   *    $             + ABS( TAU2 )*DW
   *
             IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN )              IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN )
      $         SWTCH = .NOT.SWTCH       $         SWTCH = .NOT.SWTCH
 *  *
             IF( W.LE.ZERO ) THEN  
                SG2LB = MAX( SG2LB, TAU )  
             ELSE  
                SG2UB = MIN( SG2UB, TAU )  
             END IF  
 *  
   230    CONTINUE    230    CONTINUE
 *  *
 *        Return with INFO = 1, NITER = MAXIT and not converged  *        Return with INFO = 1, NITER = MAXIT and not converged
Removed from v.1.4  
changed lines
  Added in v.1.18

CVSweb interface <joel.bertrand@systella.fr>