Diff for /rpl/lapack/lapack/dlasq2.f between versions 1.7 and 1.8

version 1.7, 2010/12/21 13:53:33 version 1.8, 2011/11/21 20:42:59
Line 1 Line 1
       SUBROUTINE DLASQ2( N, Z, INFO )  *> \brief \b DLASQ2
 *  *
 *  -- LAPACK routine (version 3.2)                                    --  *  =========== DOCUMENTATION ===========
 *  *
 *  -- Contributed by Osni Marques of the Lawrence Berkeley National   --  * Online html documentation available at 
 *  -- Laboratory and Beresford Parlett of the Univ. of California at  --  *            http://www.netlib.org/lapack/explore-html/ 
 *  -- Berkeley                                                        --  
 *  -- November 2008                                                   --  
 *  *
   *> \htmlonly
   *> Download DLASQ2 + dependencies 
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq2.f"> 
   *> [TGZ]</a> 
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq2.f"> 
   *> [ZIP]</a> 
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq2.f"> 
   *> [TXT]</a>
   *> \endhtmlonly 
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE DLASQ2( N, Z, INFO )
   * 
   *       .. Scalar Arguments ..
   *       INTEGER            INFO, N
   *       ..
   *       .. Array Arguments ..
   *       DOUBLE PRECISION   Z( * )
   *       ..
   *  
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> DLASQ2 computes all the eigenvalues of the symmetric positive 
   *> definite tridiagonal matrix associated with the qd array Z to high
   *> relative accuracy are computed to high relative accuracy, in the
   *> absence of denormalization, underflow and overflow.
   *>
   *> To see the relation of Z to the tridiagonal matrix, let L be a
   *> unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
   *> let U be an upper bidiagonal matrix with 1's above and diagonal
   *> Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
   *> symmetric tridiagonal to which it is similar.
   *>
   *> Note : DLASQ2 defines a logical variable, IEEE, which is true
   *> on machines which follow ieee-754 floating-point standard in their
   *> handling of infinities and NaNs, and false otherwise. This variable
   *> is passed to DLASQ3.
   *> \endverbatim
   *
   *  Arguments:
   *  ==========
   *
   *> \param[in] N
   *> \verbatim
   *>          N is INTEGER
   *>        The number of rows and columns in the matrix. N >= 0.
   *> \endverbatim
   *>
   *> \param[in,out] Z
   *> \verbatim
   *>          Z is DOUBLE PRECISION array, dimension ( 4*N )
   *>        On entry Z holds the qd array. On exit, entries 1 to N hold
   *>        the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
   *>        trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
   *>        N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
   *>        holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
   *>        shifts that failed.
   *> \endverbatim
   *>
   *> \param[out] INFO
   *> \verbatim
   *>          INFO is INTEGER
   *>        = 0: successful exit
   *>        < 0: if the i-th argument is a scalar and had an illegal
   *>             value, then INFO = -i, if the i-th argument is an
   *>             array and the j-entry had an illegal value, then
   *>             INFO = -(i*100+j)
   *>        > 0: the algorithm failed
   *>              = 1, a split was marked by a positive value in E
   *>              = 2, current block of Z not diagonalized after 100*N
   *>                   iterations (in inner while loop).  On exit Z holds
   *>                   a qd array with the same eigenvalues as the given Z.
   *>              = 3, termination criterion of outer while loop not met 
   *>                   (program created more than N unreduced blocks)
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee 
   *> \author Univ. of California Berkeley 
   *> \author Univ. of Colorado Denver 
   *> \author NAG Ltd. 
   *
   *> \date November 2011
   *
   *> \ingroup auxOTHERcomputational
   *
   *> \par Further Details:
   *  =====================
   *>
   *> \verbatim
   *>
   *>  Local Variables: I0:N0 defines a current unreduced segment of Z.
   *>  The shifts are accumulated in SIGMA. Iteration count is in ITER.
   *>  Ping-pong is controlled by PP (alternates between 0 and 1).
   *> \endverbatim
   *>
   *  =====================================================================
         SUBROUTINE DLASQ2( N, Z, INFO )
   *
   *  -- LAPACK computational routine (version 3.4.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 2011
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       INTEGER            INFO, N        INTEGER            INFO, N
Line 17 Line 124
       DOUBLE PRECISION   Z( * )        DOUBLE PRECISION   Z( * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  DLASQ2 computes all the eigenvalues of the symmetric positive   
 *  definite tridiagonal matrix associated with the qd array Z to high  
 *  relative accuracy are computed to high relative accuracy, in the  
 *  absence of denormalization, underflow and overflow.  
 *  
 *  To see the relation of Z to the tridiagonal matrix, let L be a  
 *  unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and  
 *  let U be an upper bidiagonal matrix with 1's above and diagonal  
 *  Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the  
 *  symmetric tridiagonal to which it is similar.  
 *  
 *  Note : DLASQ2 defines a logical variable, IEEE, which is true  
 *  on machines which follow ieee-754 floating-point standard in their  
 *  handling of infinities and NaNs, and false otherwise. This variable  
 *  is passed to DLASQ3.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  N     (input) INTEGER  
 *        The number of rows and columns in the matrix. N >= 0.  
 *  
 *  Z     (input/output) DOUBLE PRECISION array, dimension ( 4*N )  
 *        On entry Z holds the qd array. On exit, entries 1 to N hold  
 *        the eigenvalues in decreasing order, Z( 2*N+1 ) holds the  
 *        trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If  
 *        N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )  
 *        holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of  
 *        shifts that failed.  
 *  
 *  INFO  (output) INTEGER  
 *        = 0: successful exit  
 *        < 0: if the i-th argument is a scalar and had an illegal  
 *             value, then INFO = -i, if the i-th argument is an  
 *             array and the j-entry had an illegal value, then  
 *             INFO = -(i*100+j)  
 *        > 0: the algorithm failed  
 *              = 1, a split was marked by a positive value in E  
 *              = 2, current block of Z not diagonalized after 30*N  
 *                   iterations (in inner while loop)  
 *              = 3, termination criterion of outer while loop not met   
 *                   (program created more than N unreduced blocks)  
 *  
 *  Further Details  
 *  ===============  
 *  Local Variables: I0:N0 defines a current unreduced segment of Z.  
 *  The shifts are accumulated in SIGMA. Iteration count is in ITER.  
 *  Ping-pong is controlled by PP (alternates between 0 and 1).  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
Line 80 Line 135
 *     ..  *     ..
 *     .. Local Scalars ..  *     .. Local Scalars ..
       LOGICAL            IEEE        LOGICAL            IEEE
       INTEGER            I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K,        INTEGER            I0, I1, I4, IINFO, IPN4, ITER, IWHILA, IWHILB,
      $                   KMIN, N0, NBIG, NDIV, NFAIL, PP, SPLT, TTYPE       $                   K, KMIN, N0, N1, NBIG, NDIV, NFAIL, PP, SPLT, 
        $                   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,
      $                   QMIN, S, SAFMIN, SIGMA, T, TAU, TEMP, TOL,       $                   QMIN, S, SAFMIN, SIGMA, T, TAU, TEMP, TOL,
      $                   TOL2, TRACE, ZMAX       $                   TOL2, TRACE, ZMAX, TEMPE, TEMPQ
 *     ..  *     ..
 *     .. External Subroutines ..  *     .. External Subroutines ..
       EXTERNAL           DLASQ3, DLASRT, XERBLA        EXTERNAL           DLASQ3, DLASRT, XERBLA
Line 396 Line 452
 *               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 = 30*( 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
Line 441 Line 497
   140    CONTINUE    140    CONTINUE
 *  *
          INFO = 2           INFO = 2
   *       
   *        Maximum number of iterations exceeded, restore the shift 
   *        SIGMA and place the new d's and e's in a qd array.
   *        This might need to be done for several blocks
   *
            I1 = I0
            N1 = N0
    145     CONTINUE
            TEMPQ = Z( 4*I0-3 )
            Z( 4*I0-3 ) = Z( 4*I0-3 ) + SIGMA
            DO K = I0+1, N0
               TEMPE = Z( 4*K-5 )
               Z( 4*K-5 ) = Z( 4*K-5 ) * (TEMPQ / Z( 4*K-7 ))
               TEMPQ = Z( 4*K-3 )
               Z( 4*K-3 ) = Z( 4*K-3 ) + SIGMA + TEMPE - Z( 4*K-5 )
            END DO
   *
   *        Prepare to do this on the previous block if there is one
   *
            IF( I1.GT.1 ) THEN
               N1 = I1-1
               DO WHILE( ( I1.GE.2 ) .AND. ( Z(4*I1-5).GE.ZERO ) )
                  I1 = I1 - 1
               END DO
               SIGMA = -Z(4*N1-1)
               GO TO 145
            END IF
   
            DO K = 1, N
               Z( 2*K-1 ) = Z( 4*K-3 )
   *
   *        Only the block 1..N0 is unfinished.  The rest of the e's
   *        must be essentially zero, although sometimes other data
   *        has been stored in them.
   *
               IF( K.LT.N0 ) THEN
                  Z( 2*K ) = Z( 4*K-1 )
               ELSE
                  Z( 2*K ) = 0
               END IF
            END DO
          RETURN           RETURN
 *  *
 *        end IWHILB  *        end IWHILB

Removed from v.1.7  
changed lines
  Added in v.1.8


CVSweb interface <joel.bertrand@systella.fr>