--- rpl/lapack/lapack/dlarrv.f 2010/04/21 13:45:19 1.2 +++ rpl/lapack/lapack/dlarrv.f 2023/08/07 08:38:58 1.24 @@ -1,13 +1,298 @@ +*> \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 +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \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, $ ISPLIT, M, DOL, DOU, MINRGP, $ RTOL1, RTOL2, W, WERR, WGAP, $ IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ, $ WORK, IWORK, INFO ) * -* -- LAPACK auxiliary routine (version 3.2) -- +* -- LAPACK auxiliary routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 * * .. Scalar Arguments .. INTEGER DOL, DOU, INFO, LDZ, M, N @@ -21,153 +306,6 @@ 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 .. @@ -208,6 +346,14 @@ * .. 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 INDLD = N+1 INDLLD= 2*N+1 @@ -332,7 +478,7 @@ * high relative accuracy is required for the computation of the * corresponding eigenvectors. CALL DCOPY( IM, W( WBEGIN ), 1, - & WORK( WBEGIN ), 1 ) + $ WORK( WBEGIN ), 1 ) * We store in W the eigenvalue approximations w.r.t. the original * matrix T. @@ -429,9 +575,9 @@ * within the current block P = INDEXW( WBEGIN-1+OLDFST ) Q = INDEXW( WBEGIN-1+OLDLST ) -* Offset for the arrays WORK, WGAP and WERR, i.e., th P-OFFSET -* thru' Q-OFFSET elements of these arrays are to be used. -C OFFSET = P-OLDFST +* Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET +* through the Q-OFFSET elements of these arrays are to be used. +* OFFSET = P-OLDFST OFFSET = INDEXW( WBEGIN ) - 1 * perform limited bisection (if necessary) to get approximate * eigenvalues to the precision needed. @@ -570,7 +716,7 @@ C OFFSET = P-OLDFST * Compute RRR of child cluster. * 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 ), $ WORK(INDLD+IBEGIN-1), $ NEWFST, NEWLST, WORK(WBEGIN),