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

Diff for /rpl/lapack/lapack/dlarrv.f between versions 1.1.1.1 and 1.24

version 1.1.1.1, 2010/01/26 15:22:46 version 1.24, 2023/08/07 08:38:58
Line 1 Line 1
   *> \brief \b DLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT.
   *
   *  =========== DOCUMENTATION ===========
   *
   * Online html documentation available at
   *            http://www.netlib.org/lapack/explore-html/
   *
   *> \htmlonly
   *> Download DLARRV + dependencies
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrv.f">
   *> [TGZ]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrv.f">
   *> [ZIP]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrv.f">
   *> [TXT]</a>
   *> \endhtmlonly
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE DLARRV( N, VL, VU, D, L, PIVMIN,
   *                          ISPLIT, M, DOL, DOU, MINRGP,
   *                          RTOL1, RTOL2, W, WERR, WGAP,
   *                          IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
   *                          WORK, IWORK, INFO )
   *
   *       .. Scalar Arguments ..
   *       INTEGER            DOL, DOU, INFO, LDZ, M, N
   *       DOUBLE PRECISION   MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
   *       ..
   *       .. Array Arguments ..
   *       INTEGER            IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
   *      $                   ISUPPZ( * ), IWORK( * )
   *       DOUBLE PRECISION   D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
   *      $                   WGAP( * ), WORK( * )
   *       DOUBLE PRECISION  Z( LDZ, * )
   *       ..
   *
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> DLARRV computes the eigenvectors of the tridiagonal matrix
   *> T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T.
   *> The input eigenvalues should have been computed by DLARRE.
   *> \endverbatim
   *
   *  Arguments:
   *  ==========
   *
   *> \param[in] N
   *> \verbatim
   *>          N is INTEGER
   *>          The order of the matrix.  N >= 0.
   *> \endverbatim
   *>
   *> \param[in] VL
   *> \verbatim
   *>          VL is DOUBLE PRECISION
   *>          Lower bound of the interval that contains the desired
   *>          eigenvalues. VL < VU. Needed to compute gaps on the left or right
   *>          end of the extremal eigenvalues in the desired RANGE.
   *> \endverbatim
   *>
   *> \param[in] VU
   *> \verbatim
   *>          VU is DOUBLE PRECISION
   *>          Upper bound of the interval that contains the desired
   *>          eigenvalues. VL < VU. 
   *>          Note: VU is currently not used by this implementation of DLARRV, VU is
   *>          passed to DLARRV because it could be used compute gaps on the right end
   *>          of the extremal eigenvalues. However, with not much initial accuracy in
   *>          LAMBDA and VU, the formula can lead to an overestimation of the right gap
   *>          and thus to inadequately early RQI 'convergence'. This is currently
   *>          prevented this by forcing a small right gap. And so it turns out that VU
   *>          is currently not used by this implementation of DLARRV.
   *> \endverbatim
   *>
   *> \param[in,out] D
   *> \verbatim
   *>          D is DOUBLE PRECISION array, dimension (N)
   *>          On entry, the N diagonal elements of the diagonal matrix D.
   *>          On exit, D may be overwritten.
   *> \endverbatim
   *>
   *> \param[in,out] L
   *> \verbatim
   *>          L is DOUBLE PRECISION array, dimension (N)
   *>          On entry, the (N-1) subdiagonal elements of the unit
   *>          bidiagonal matrix L are in elements 1 to N-1 of L
   *>          (if the matrix is not split.) At the end of each block
   *>          is stored the corresponding shift as given by DLARRE.
   *>          On exit, L is overwritten.
   *> \endverbatim
   *>
   *> \param[in] PIVMIN
   *> \verbatim
   *>          PIVMIN is DOUBLE PRECISION
   *>          The minimum pivot allowed in the Sturm sequence.
   *> \endverbatim
   *>
   *> \param[in] ISPLIT
   *> \verbatim
   *>          ISPLIT is INTEGER array, dimension (N)
   *>          The splitting points, at which T breaks up into blocks.
   *>          The first block consists of rows/columns 1 to
   *>          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
   *>          through ISPLIT( 2 ), etc.
   *> \endverbatim
   *>
   *> \param[in] M
   *> \verbatim
   *>          M is INTEGER
   *>          The total number of input eigenvalues.  0 <= M <= N.
   *> \endverbatim
   *>
   *> \param[in] DOL
   *> \verbatim
   *>          DOL is INTEGER
   *> \endverbatim
   *>
   *> \param[in] DOU
   *> \verbatim
   *>          DOU is INTEGER
   *>          If the user wants to compute only selected eigenvectors from all
   *>          the eigenvalues supplied, he can specify an index range DOL:DOU.
   *>          Or else the setting DOL=1, DOU=M should be applied.
   *>          Note that DOL and DOU refer to the order in which the eigenvalues
   *>          are stored in W.
   *>          If the user wants to compute only selected eigenpairs, then
   *>          the columns DOL-1 to DOU+1 of the eigenvector space Z contain the
   *>          computed eigenvectors. All other columns of Z are set to zero.
   *> \endverbatim
   *>
   *> \param[in] MINRGP
   *> \verbatim
   *>          MINRGP is DOUBLE PRECISION
   *> \endverbatim
   *>
   *> \param[in] RTOL1
   *> \verbatim
   *>          RTOL1 is DOUBLE PRECISION
   *> \endverbatim
   *>
   *> \param[in] RTOL2
   *> \verbatim
   *>          RTOL2 is DOUBLE PRECISION
   *>           Parameters for bisection.
   *>           An interval [LEFT,RIGHT] has converged if
   *>           RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
   *> \endverbatim
   *>
   *> \param[in,out] W
   *> \verbatim
   *>          W is DOUBLE PRECISION array, dimension (N)
   *>          The first M elements of W contain the APPROXIMATE eigenvalues for
   *>          which eigenvectors are to be computed.  The eigenvalues
   *>          should be grouped by split-off block and ordered from
   *>          smallest to largest within the block ( The output array
   *>          W from DLARRE is expected here ). Furthermore, they are with
   *>          respect to the shift of the corresponding root representation
   *>          for their block. On exit, W holds the eigenvalues of the
   *>          UNshifted matrix.
   *> \endverbatim
   *>
   *> \param[in,out] WERR
   *> \verbatim
   *>          WERR is DOUBLE PRECISION array, dimension (N)
   *>          The first M elements contain the semiwidth of the uncertainty
   *>          interval of the corresponding eigenvalue in W
   *> \endverbatim
   *>
   *> \param[in,out] WGAP
   *> \verbatim
   *>          WGAP is DOUBLE PRECISION array, dimension (N)
   *>          The separation from the right neighbor eigenvalue in W.
   *> \endverbatim
   *>
   *> \param[in] IBLOCK
   *> \verbatim
   *>          IBLOCK is INTEGER array, dimension (N)
   *>          The indices of the blocks (submatrices) associated with the
   *>          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
   *>          W(i) belongs to the first block from the top, =2 if W(i)
   *>          belongs to the second block, etc.
   *> \endverbatim
   *>
   *> \param[in] INDEXW
   *> \verbatim
   *>          INDEXW is INTEGER array, dimension (N)
   *>          The indices of the eigenvalues within each block (submatrix);
   *>          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
   *>          i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.
   *> \endverbatim
   *>
   *> \param[in] GERS
   *> \verbatim
   *>          GERS is DOUBLE PRECISION array, dimension (2*N)
   *>          The N Gerschgorin intervals (the i-th Gerschgorin interval
   *>          is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should
   *>          be computed from the original UNshifted matrix.
   *> \endverbatim
   *>
   *> \param[out] Z
   *> \verbatim
   *>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
   *>          If INFO = 0, the first M columns of Z contain the
   *>          orthonormal eigenvectors of the matrix T
   *>          corresponding to the input eigenvalues, with the i-th
   *>          column of Z holding the eigenvector associated with W(i).
   *>          Note: the user must ensure that at least max(1,M) columns are
   *>          supplied in the array Z.
   *> \endverbatim
   *>
   *> \param[in] LDZ
   *> \verbatim
   *>          LDZ is INTEGER
   *>          The leading dimension of the array Z.  LDZ >= 1, and if
   *>          JOBZ = 'V', LDZ >= max(1,N).
   *> \endverbatim
   *>
   *> \param[out] ISUPPZ
   *> \verbatim
   *>          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
   *>          The support of the eigenvectors in Z, i.e., the indices
   *>          indicating the nonzero elements in Z. The I-th eigenvector
   *>          is nonzero only in elements ISUPPZ( 2*I-1 ) through
   *>          ISUPPZ( 2*I ).
   *> \endverbatim
   *>
   *> \param[out] WORK
   *> \verbatim
   *>          WORK is DOUBLE PRECISION array, dimension (12*N)
   *> \endverbatim
   *>
   *> \param[out] IWORK
   *> \verbatim
   *>          IWORK is INTEGER array, dimension (7*N)
   *> \endverbatim
   *>
   *> \param[out] INFO
   *> \verbatim
   *>          INFO is INTEGER
   *>          = 0:  successful exit
   *>
   *>          > 0:  A problem occurred in DLARRV.
   *>          < 0:  One of the called subroutines signaled an internal problem.
   *>                Needs inspection of the corresponding parameter IINFO
   *>                for further information.
   *>
   *>          =-1:  Problem in DLARRB when refining a child's eigenvalues.
   *>          =-2:  Problem in DLARRF when computing the RRR of a child.
   *>                When a child is inside a tight cluster, it can be difficult
   *>                to find an RRR. A partial remedy from the user's point of
   *>                view is to make the parameter MINRGP smaller and recompile.
   *>                However, as the orthogonality of the computed vectors is
   *>                proportional to 1/MINRGP, the user should be aware that
   *>                he might be trading in precision when he decreases MINRGP.
   *>          =-3:  Problem in DLARRB when refining a single eigenvalue
   *>                after the Rayleigh correction was rejected.
   *>          = 5:  The Rayleigh Quotient Iteration failed to converge to
   *>                full accuracy in MAXITR steps.
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee
   *> \author Univ. of California Berkeley
   *> \author Univ. of Colorado Denver
   *> \author NAG Ltd.
   *
   *> \ingroup doubleOTHERauxiliary
   *
   *> \par Contributors:
   *  ==================
   *>
   *> Beresford Parlett, University of California, Berkeley, USA \n
   *> Jim Demmel, University of California, Berkeley, USA \n
   *> Inderjit Dhillon, University of Texas, Austin, USA \n
   *> Osni Marques, LBNL/NERSC, USA \n
   *> Christof Voemel, University of California, Berkeley, USA
   *
   *  =====================================================================
       SUBROUTINE DLARRV( N, VL, VU, D, L, PIVMIN,        SUBROUTINE DLARRV( N, VL, VU, D, L, PIVMIN,
      $                   ISPLIT, M, DOL, DOU, MINRGP,       $                   ISPLIT, M, DOL, DOU, MINRGP,
      $                   RTOL1, RTOL2, W, WERR, WGAP,       $                   RTOL1, RTOL2, W, WERR, WGAP,
      $                   IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,       $                   IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
      $                   WORK, IWORK, INFO )       $                   WORK, IWORK, INFO )
 *  *
 *  -- LAPACK auxiliary routine (version 3.2) --  *  -- LAPACK auxiliary routine --
 *  -- 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  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       INTEGER            DOL, DOU, INFO, LDZ, M, N        INTEGER            DOL, DOU, INFO, LDZ, M, N
Line 21 Line 306
       DOUBLE PRECISION  Z( LDZ, * )        DOUBLE PRECISION  Z( LDZ, * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  DLARRV computes the eigenvectors of the tridiagonal matrix  
 *  T = L D L^T given L, D and APPROXIMATIONS to the eigenvalues of L D L^T.  
 *  The input eigenvalues should have been computed by DLARRE.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  N       (input) INTEGER  
 *          The order of the matrix.  N >= 0.  
 *  
 *  VL      (input) DOUBLE PRECISION  
 *  VU      (input) DOUBLE PRECISION  
 *          Lower and upper bounds of the interval that contains the desired  
 *          eigenvalues. VL < VU. Needed to compute gaps on the left or right  
 *          end of the extremal eigenvalues in the desired RANGE.  
 *  
 *  D       (input/output) DOUBLE PRECISION array, dimension (N)  
 *          On entry, the N diagonal elements of the diagonal matrix D.  
 *          On exit, D may be overwritten.  
 *  
 *  L       (input/output) DOUBLE PRECISION array, dimension (N)  
 *          On entry, the (N-1) subdiagonal elements of the unit  
 *          bidiagonal matrix L are in elements 1 to N-1 of L  
 *          (if the matrix is not splitted.) At the end of each block  
 *          is stored the corresponding shift as given by DLARRE.  
 *          On exit, L is overwritten.  
 *  
 *  PIVMIN  (in) DOUBLE PRECISION  
 *          The minimum pivot allowed in the Sturm sequence.  
 *  
 *  ISPLIT  (input) INTEGER array, dimension (N)  
 *          The splitting points, at which T breaks up into blocks.  
 *          The first block consists of rows/columns 1 to  
 *          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1  
 *          through ISPLIT( 2 ), etc.  
 *  
 *  M       (input) INTEGER  
 *          The total number of input eigenvalues.  0 <= M <= N.  
 *  
 *  DOL     (input) INTEGER  
 *  DOU     (input) INTEGER  
 *          If the user wants to compute only selected eigenvectors from all  
 *          the eigenvalues supplied, he can specify an index range DOL:DOU.  
 *          Or else the setting DOL=1, DOU=M should be applied.  
 *          Note that DOL and DOU refer to the order in which the eigenvalues  
 *          are stored in W.  
 *          If the user wants to compute only selected eigenpairs, then  
 *          the columns DOL-1 to DOU+1 of the eigenvector space Z contain the  
 *          computed eigenvectors. All other columns of Z are set to zero.  
 *  
 *  MINRGP  (input) DOUBLE PRECISION  
 *  
 *  RTOL1   (input) DOUBLE PRECISION  
 *  RTOL2   (input) DOUBLE PRECISION  
 *           Parameters for bisection.  
 *           An interval [LEFT,RIGHT] has converged if  
 *           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )  
 *  
 *  W       (input/output) DOUBLE PRECISION array, dimension (N)  
 *          The first M elements of W contain the APPROXIMATE eigenvalues for  
 *          which eigenvectors are to be computed.  The eigenvalues  
 *          should be grouped by split-off block and ordered from  
 *          smallest to largest within the block ( The output array  
 *          W from DLARRE is expected here ). Furthermore, they are with  
 *          respect to the shift of the corresponding root representation  
 *          for their block. On exit, W holds the eigenvalues of the  
 *          UNshifted matrix.  
 *  
 *  WERR    (input/output) DOUBLE PRECISION array, dimension (N)  
 *          The first M elements contain the semiwidth of the uncertainty  
 *          interval of the corresponding eigenvalue in W  
 *  
 *  WGAP    (input/output) DOUBLE PRECISION array, dimension (N)  
 *          The separation from the right neighbor eigenvalue in W.  
 *  
 *  IBLOCK  (input) INTEGER array, dimension (N)  
 *          The indices of the blocks (submatrices) associated with the  
 *          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue  
 *          W(i) belongs to the first block from the top, =2 if W(i)  
 *          belongs to the second block, etc.  
 *  
 *  INDEXW  (input) INTEGER array, dimension (N)  
 *          The indices of the eigenvalues within each block (submatrix);  
 *          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the  
 *          i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.  
 *  
 *  GERS    (input) DOUBLE PRECISION array, dimension (2*N)  
 *          The N Gerschgorin intervals (the i-th Gerschgorin interval  
 *          is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should  
 *          be computed from the original UNshifted matrix.  
 *  
 *  Z       (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )  
 *          If INFO = 0, the first M columns of Z contain the  
 *          orthonormal eigenvectors of the matrix T  
 *          corresponding to the input eigenvalues, with the i-th  
 *          column of Z holding the eigenvector associated with W(i).  
 *          Note: the user must ensure that at least max(1,M) columns are  
 *          supplied in the array Z.  
 *  
 *  LDZ     (input) INTEGER  
 *          The leading dimension of the array Z.  LDZ >= 1, and if  
 *          JOBZ = 'V', LDZ >= max(1,N).  
 *  
 *  ISUPPZ  (output) INTEGER array, dimension ( 2*max(1,M) )  
 *          The support of the eigenvectors in Z, i.e., the indices  
 *          indicating the nonzero elements in Z. The I-th eigenvector  
 *          is nonzero only in elements ISUPPZ( 2*I-1 ) through  
 *          ISUPPZ( 2*I ).  
 *  
 *  WORK    (workspace) DOUBLE PRECISION array, dimension (12*N)  
 *  
 *  IWORK   (workspace) INTEGER array, dimension (7*N)  
 *  
 *  INFO    (output) INTEGER  
 *          = 0:  successful exit  
 *  
 *          > 0:  A problem occured in DLARRV.  
 *          < 0:  One of the called subroutines signaled an internal problem.  
 *                Needs inspection of the corresponding parameter IINFO  
 *                for further information.  
 *  
 *          =-1:  Problem in DLARRB when refining a child's eigenvalues.  
 *          =-2:  Problem in DLARRF when computing the RRR of a child.  
 *                When a child is inside a tight cluster, it can be difficult  
 *                to find an RRR. A partial remedy from the user's point of  
 *                view is to make the parameter MINRGP smaller and recompile.  
 *                However, as the orthogonality of the computed vectors is  
 *                proportional to 1/MINRGP, the user should be aware that  
 *                he might be trading in precision when he decreases MINRGP.  
 *          =-3:  Problem in DLARRB when refining a single eigenvalue  
 *                after the Rayleigh correction was rejected.  
 *          = 5:  The Rayleigh Quotient Iteration failed to converge to  
 *                full accuracy in MAXITR steps.  
 *  
 *  Further Details  
 *  ===============  
 *  
 *  Based on contributions by  
 *     Beresford Parlett, University of California, Berkeley, USA  
 *     Jim Demmel, University of California, Berkeley, USA  
 *     Inderjit Dhillon, University of Texas, Austin, USA  
 *     Osni Marques, LBNL/NERSC, USA  
 *     Christof Voemel, University of California, Berkeley, USA  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
Line 208 Line 346
 *     .. Executable Statements ..  *     .. Executable Statements ..
 *     ..  *     ..
   
         INFO = 0
   *
   *     Quick return if possible
   *
         IF( (N.LE.0).OR.(M.LE.0) ) THEN
            RETURN
         END IF
   *
 *     The first N entries of WORK are reserved for the eigenvalues  *     The first N entries of WORK are reserved for the eigenvalues
       INDLD = N+1        INDLD = N+1
       INDLLD= 2*N+1        INDLLD= 2*N+1
Line 332 Line 478
 *        high relative accuracy is required for the computation of the  *        high relative accuracy is required for the computation of the
 *        corresponding eigenvectors.  *        corresponding eigenvectors.
          CALL DCOPY( IM, W( WBEGIN ), 1,           CALL DCOPY( IM, W( WBEGIN ), 1,
      &                   WORK( WBEGIN ), 1 )       $                   WORK( WBEGIN ), 1 )
   
 *        We store in W the eigenvalue approximations w.r.t. the original  *        We store in W the eigenvalue approximations w.r.t. the original
 *        matrix T.  *        matrix T.
Line 429 Line 575
 *                 within the current block  *                 within the current block
                   P = INDEXW( WBEGIN-1+OLDFST )                    P = INDEXW( WBEGIN-1+OLDFST )
                   Q = INDEXW( WBEGIN-1+OLDLST )                    Q = INDEXW( WBEGIN-1+OLDLST )
 *                 Offset for the arrays WORK, WGAP and WERR, i.e., th P-OFFSET  *                 Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET
 *                 thru' Q-OFFSET elements of these arrays are to be used.  *                 through the Q-OFFSET elements of these arrays are to be used.
 C                  OFFSET = P-OLDFST  *                  OFFSET = P-OLDFST
                   OFFSET = INDEXW( WBEGIN ) - 1                    OFFSET = INDEXW( WBEGIN ) - 1
 *                 perform limited bisection (if necessary) to get approximate  *                 perform limited bisection (if necessary) to get approximate
 *                 eigenvalues to the precision needed.  *                 eigenvalues to the precision needed.
Line 570  C                  OFFSET = P-OLDFST Line 716  C                  OFFSET = P-OLDFST
 *                    Compute RRR of child cluster.  *                    Compute RRR of child cluster.
 *                    Note that the new RRR is stored in Z  *                    Note that the new RRR is stored in Z
 *  *
 C                    DLARRF needs LWORK = 2*N  *                    DLARRF needs LWORK = 2*N
                      CALL DLARRF( IN, D( IBEGIN ), L( IBEGIN ),                       CALL DLARRF( IN, D( IBEGIN ), L( IBEGIN ),
      $                         WORK(INDLD+IBEGIN-1),       $                         WORK(INDLD+IBEGIN-1),
      $                         NEWFST, NEWLST, WORK(WBEGIN),       $                         NEWFST, NEWLST, WORK(WBEGIN),
Removed from v.1.1.1.1  
changed lines
  Added in v.1.24

CVSweb interface <joel.bertrand@systella.fr>