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

Diff for /rpl/lapack/lapack/dlasq2.f between versions 1.14 and 1.15

version 1.14, 2016/08/27 15:34:31 version 1.15, 2017/06/17 10:53:57
Line 2 Line 2
 *  *
 *  =========== 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 DLASQ2 + dependencies   *> Download DLASQ2 + dependencies
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq2.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq2.f">
 *> [TGZ]</a>   *> [TGZ]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq2.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq2.f">
 *> [ZIP]</a>   *> [ZIP]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq2.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq2.f">
 *> [TXT]</a>  *> [TXT]</a>
 *> \endhtmlonly   *> \endhtmlonly
 *  *
 *  Definition:  *  Definition:
 *  ===========  *  ===========
 *  *
 *       SUBROUTINE DLASQ2( N, Z, INFO )  *       SUBROUTINE DLASQ2( N, Z, INFO )
 *   *
 *       .. Scalar Arguments ..  *       .. Scalar Arguments ..
 *       INTEGER            INFO, N  *       INTEGER            INFO, N
 *       ..  *       ..
 *       .. Array Arguments ..  *       .. Array Arguments ..
 *       DOUBLE PRECISION   Z( * )  *       DOUBLE PRECISION   Z( * )
 *       ..  *       ..
 *    *
 *  *
 *> \par Purpose:  *> \par Purpose:
 *  =============  *  =============
 *>  *>
 *> \verbatim  *> \verbatim
 *>  *>
 *> DLASQ2 computes all the eigenvalues of the symmetric positive   *> DLASQ2 computes all the eigenvalues of the symmetric positive
 *> definite tridiagonal matrix associated with the qd array Z to high  *> definite tridiagonal matrix associated with the qd array Z to high
 *> relative accuracy are computed to high relative accuracy, in the  *> relative accuracy are computed to high relative accuracy, in the
 *> absence of denormalization, underflow and overflow.  *> absence of denormalization, underflow and overflow.
Line 83 Line 83
 *>              = 2, current block of Z not diagonalized after 100*N  *>              = 2, current block of Z not diagonalized after 100*N
 *>                   iterations (in inner while loop).  On exit Z holds  *>                   iterations (in inner while loop).  On exit Z holds
 *>                   a qd array with the same eigenvalues as the given Z.  *>                   a qd array with the same eigenvalues as the given Z.
 *>              = 3, termination criterion of outer while loop not met   *>              = 3, termination criterion of outer while loop not met
 *>                   (program created more than N unreduced blocks)  *>                   (program created more than N unreduced blocks)
 *> \endverbatim  *> \endverbatim
 *  *
 *  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 September 2012  *> \date December 2016
 *  *
 *> \ingroup auxOTHERcomputational  *> \ingroup auxOTHERcomputational
 *  *
Line 112 Line 112
 *  =====================================================================  *  =====================================================================
       SUBROUTINE DLASQ2( N, Z, INFO )        SUBROUTINE DLASQ2( N, Z, INFO )
 *  *
 *  -- LAPACK computational routine (version 3.4.2) --  *  -- LAPACK computational 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..--
 *     September 2012  *     December 2016
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       INTEGER            INFO, N        INTEGER            INFO, N
Line 136 Line 136
 *     .. Local Scalars ..  *     .. Local Scalars ..
       LOGICAL            IEEE        LOGICAL            IEEE
       INTEGER            I0, I1, I4, IINFO, IPN4, ITER, IWHILA, IWHILB,        INTEGER            I0, I1, I4, IINFO, IPN4, ITER, IWHILA, IWHILB,
      $                   K, KMIN, N0, N1, NBIG, NDIV, NFAIL, PP, SPLT,        $                   K, KMIN, N0, N1, NBIG, NDIV, NFAIL, PP, SPLT,
      $                   TTYPE       $                   TTYPE
       DOUBLE PRECISION   D, DEE, DEEMIN, DESIG, DMIN, DMIN1, DMIN2, DN,        DOUBLE PRECISION   D, DEE, DEEMIN, DESIG, DMIN, DMIN1, DMIN2, DN,
      $                   DN1, DN2, E, EMAX, EMIN, EPS, G, OLDEMN, QMAX,       $                   DN1, DN2, E, EMAX, EMIN, EPS, G, OLDEMN, QMAX,
Line 155 Line 155
       INTRINSIC          ABS, DBLE, MAX, MIN, SQRT        INTRINSIC          ABS, DBLE, MAX, MIN, SQRT
 *     ..  *     ..
 *     .. Executable Statements ..  *     .. Executable Statements ..
 *        *
 *     Test the input arguments.  *     Test the input arguments.
 *     (in case DLASQ2 is not called by DLASQ1)  *     (in case DLASQ2 is not called by DLASQ1)
 *  *
Line 195 Line 195
          END IF           END IF
          Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 )           Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 )
          IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN           IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN
             T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) )               T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) )
             S = Z( 3 )*( Z( 2 ) / T )              S = Z( 3 )*( Z( 2 ) / T )
             IF( S.LE.T ) THEN              IF( S.LE.T ) THEN
                S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) )                 S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) )
Line 264 Line 264
          Z( 2*N-1 ) = ZERO           Z( 2*N-1 ) = ZERO
          RETURN           RETURN
       END IF        END IF
 *           *
 *     Check whether the machine is IEEE conformable.  *     Check whether the machine is IEEE conformable.
 *           *
       IEEE = ILAENV( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND.        IEEE = ILAENV( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND.
      $       ILAENV( 11, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1             $       ILAENV( 11, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1
 *           *
 *     Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).  *     Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
 *  *
       DO 30 K = 2*N, 2, -2        DO 30 K = 2*N, 2, -2
          Z( 2*K ) = ZERO            Z( 2*K ) = ZERO
          Z( 2*K-1 ) = Z( K )            Z( 2*K-1 ) = Z( K )
          Z( 2*K-2 ) = ZERO            Z( 2*K-2 ) = ZERO
          Z( 2*K-3 ) = Z( K-1 )            Z( 2*K-3 ) = Z( K-1 )
    30 CONTINUE     30 CONTINUE
 *  *
       I0 = 1        I0 = 1
Line 333 Line 333
                D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) )                 D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) )
             END IF              END IF
             EMIN = MIN( EMIN, Z( I4-2*PP ) )              EMIN = MIN( EMIN, Z( I4-2*PP ) )
    60    CONTINUE      60    CONTINUE
          Z( 4*N0-PP-2 ) = D           Z( 4*N0-PP-2 ) = D
 *  *
 *        Now find qmax.  *        Now find qmax.
Line 364 Line 364
       NDIV = 2*( N0-I0 )        NDIV = 2*( N0-I0 )
 *  *
       DO 160 IWHILA = 1, N + 1        DO 160 IWHILA = 1, N + 1
          IF( N0.LT.1 )            IF( N0.LT.1 )
      $      GO TO 170       $      GO TO 170
 *  *
 *        While array unfinished do   *        While array unfinished do
 *  *
 *        E(N0) holds the value of SIGMA when submatrix in I0:N0  *        E(N0) holds the value of SIGMA when submatrix in I0:N0
 *        splits from the rest of the array, but is negated.  *        splits from the rest of the array, but is negated.
 *        *
          DESIG = ZERO           DESIG = ZERO
          IF( N0.EQ.N ) THEN           IF( N0.EQ.N ) THEN
             SIGMA = ZERO              SIGMA = ZERO
Line 386 Line 386
 *        Find last unreduced submatrix's top index I0, find QMAX and  *        Find last unreduced submatrix's top index I0, find QMAX and
 *        EMIN. Find Gershgorin-type bound if Q's much greater than E's.  *        EMIN. Find Gershgorin-type bound if Q's much greater than E's.
 *  *
          EMAX = ZERO            EMAX = ZERO
          IF( N0.GT.I0 ) THEN           IF( N0.GT.I0 ) THEN
             EMIN = ABS( Z( 4*N0-5 ) )              EMIN = ABS( Z( 4*N0-5 ) )
          ELSE           ELSE
Line 404 Line 404
             QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) )              QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) )
             EMIN = MIN( EMIN, Z( I4-5 ) )              EMIN = MIN( EMIN, Z( I4-5 ) )
    90    CONTINUE     90    CONTINUE
          I4 = 4            I4 = 4
 *  *
   100    CONTINUE    100    CONTINUE
          I0 = I4 / 4           I0 = I4 / 4
Line 421 Line 421
                   KMIN = ( I4+3 )/4                    KMIN = ( I4+3 )/4
                END IF                 END IF
   110       CONTINUE    110       CONTINUE
             IF( (KMIN-I0)*2.LT.N0-KMIN .AND.               IF( (KMIN-I0)*2.LT.N0-KMIN .AND.
      $         DEEMIN.LE.HALF*Z(4*N0-3) ) THEN       $         DEEMIN.LE.HALF*Z(4*N0-3) ) THEN
                IPN4 = 4*( I0+N0 )                 IPN4 = 4*( I0+N0 )
                PP = 2                 PP = 2
Line 446 Line 446
 *  *
          DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) )           DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) )
 *  *
 *        Now I0:N0 is unreduced.   *        Now I0:N0 is unreduced.
 *        PP = 0 for ping, PP = 1 for pong.  *        PP = 0 for ping, PP = 1 for pong.
 *        PP = 2 indicates that flipping was applied to the Z array and  *        PP = 2 indicates that flipping was applied to the Z array and
 *               and that the tests for deflation upon entry in DLASQ3   *               and that the tests for deflation upon entry in DLASQ3
 *               should not be performed.  *               should not be performed.
 *  *
          NBIG = 100*( N0-I0+1 )           NBIG = 100*( N0-I0+1 )
          DO 140 IWHILB = 1, NBIG           DO 140 IWHILB = 1, NBIG
             IF( I0.GT.N0 )               IF( I0.GT.N0 )
      $         GO TO 150       $         GO TO 150
 *  *
 *           While submatrix unfinished take a good dqds step.  *           While submatrix unfinished take a good dqds step.
Line 497 Line 497
   140    CONTINUE    140    CONTINUE
 *  *
          INFO = 2           INFO = 2
 *         *
 *        Maximum number of iterations exceeded, restore the shift   *        Maximum number of iterations exceeded, restore the shift
 *        SIGMA and place the new d's and e's in a qd array.  *        SIGMA and place the new d's and e's in a qd array.
 *        This might need to be done for several blocks  *        This might need to be done for several blocks
 *  *
Line 549 Line 549
       INFO = 3        INFO = 3
       RETURN        RETURN
 *  *
 *     end IWHILA     *     end IWHILA
 *  *
   170 CONTINUE    170 CONTINUE
 *        *
 *     Move q's to the front.  *     Move q's to the front.
 *        *
       DO 180 K = 2, N        DO 180 K = 2, N
          Z( K ) = Z( 4*K-3 )           Z( K ) = Z( 4*K-3 )
   180 CONTINUE    180 CONTINUE
 *        *
 *     Sort and compute sum of eigenvalues.  *     Sort and compute sum of eigenvalues.
 *  *
       CALL DLASRT( 'D', N, Z, IINFO )        CALL DLASRT( 'D', N, Z, IINFO )
Line 570 Line 570
 *  *
 *     Store trace, sum(eigenvalues) and information on performance.  *     Store trace, sum(eigenvalues) and information on performance.
 *  *
       Z( 2*N+1 ) = TRACE         Z( 2*N+1 ) = TRACE
       Z( 2*N+2 ) = E        Z( 2*N+2 ) = E
       Z( 2*N+3 ) = DBLE( ITER )        Z( 2*N+3 ) = DBLE( ITER )
       Z( 2*N+4 ) = DBLE( NDIV ) / DBLE( N**2 )        Z( 2*N+4 ) = DBLE( NDIV ) / DBLE( N**2 )
Removed from v.1.14  
changed lines
  Added in v.1.15

CVSweb interface <joel.bertrand@systella.fr>