[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.12 and 1.13

version 1.12, 2012/08/22 09:48:20 version 1.13, 2012/12/14 12:30:25
Line 1 Line 1
 *> \brief \b DLASD4  *> \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 ===========  *  =========== DOCUMENTATION ===========
 *  *
 * Online html documentation available at   * Online html documentation available at
 *            http://www.netlib.org/lapack/explore-html/   *            http://www.netlib.org/lapack/explore-html/
 *  *
 *> \htmlonly  *> \htmlonly
 *> Download DLASD4 + dependencies   *> Download DLASD4 + dependencies
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd4.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd4.f">
 *> [TGZ]</a>   *> [TGZ]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd4.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd4.f">
 *> [ZIP]</a>   *> [ZIP]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd4.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd4.f">
 *> [TXT]</a>  *> [TXT]</a>
 *> \endhtmlonly   *> \endhtmlonly
 *  *
 *  Definition:  *  Definition:
 *  ===========  *  ===========
 *  *
 *       SUBROUTINE DLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )  *       SUBROUTINE DLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )
 *   *
 *       .. Scalar Arguments ..  *       .. Scalar Arguments ..
 *       INTEGER            I, INFO, N  *       INTEGER            I, INFO, N
 *       DOUBLE PRECISION   RHO, SIGMA  *       DOUBLE PRECISION   RHO, SIGMA
Line 27 Line 27
 *       .. Array Arguments ..  *       .. Array Arguments ..
 *       DOUBLE PRECISION   D( * ), DELTA( * ), WORK( * ), Z( * )  *       DOUBLE PRECISION   D( * ), DELTA( * ), WORK( * ), Z( * )
 *       ..  *       ..
 *    *
 *  *
 *> \par Purpose:  *> \par Purpose:
 *  =============  *  =============
Line 135 Line 135
 *  Authors:  *  Authors:
 *  ========  *  ========
 *  *
 *> \author Univ. of Tennessee   *> \author Univ. of Tennessee
 *> \author Univ. of California Berkeley   *> \author Univ. of California Berkeley
 *> \author Univ. of Colorado Denver   *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.   *> \author NAG Ltd.
 *  *
 *> \date November 2011  *> \date September 2012
 *  *
 *> \ingroup auxOTHERauxiliary  *> \ingroup auxOTHERauxiliary
 *  *
Line 153 Line 153
 *  =====================================================================  *  =====================================================================
       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.4.0) --  *  -- LAPACK auxiliary routine (version 3.4.2) --
 *  -- 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 2011  *     September 2012
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       INTEGER            I, INFO, N        INTEGER            I, INFO, N
Line 170 Line 170
 *  *
 *     .. Parameters ..  *     .. Parameters ..
       INTEGER            MAXIT        INTEGER            MAXIT
       PARAMETER          ( MAXIT = 64 )        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 261 Line 261
      $             ( 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 271 Line 271
                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
             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
 *  *
 *           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 327 Line 327
          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 368 Line 368
          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 392 Line 392
 *  *
 *        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 437 Line 438
             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 461 Line 462
 *  *
 *           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 488 Line 490
 *  *
          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 507 Line 510
          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 514 Line 518
 *           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 536 Line 547
 *           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 616 Line 618
          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 626 Line 629
          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 693 Line 696
             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 705 Line 738
 *  *
          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 755 Line 794
             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 786 Line 821
 *           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 873 Line 915
                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 885 Line 973
 *  *
             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 935 Line 1029
                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.12  
changed lines
  Added in v.1.13

CVSweb interface <joel.bertrand@systella.fr>