--- rpl/lapack/lapack/dlarre.f 2010/01/26 15:22:45 1.1 +++ rpl/lapack/lapack/dlarre.f 2023/08/07 08:38:57 1.24 @@ -1,13 +1,311 @@ +*> \brief \b DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues. +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLARRE + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2, +* RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, +* W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, +* WORK, IWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER RANGE +* INTEGER IL, INFO, IU, M, N, NSPLIT +* DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU +* .. +* .. Array Arguments .. +* INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ), +* $ INDEXW( * ) +* DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ), +* $ W( * ),WERR( * ), WGAP( * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> To find the desired eigenvalues of a given real symmetric +*> tridiagonal matrix T, DLARRE sets any "small" off-diagonal +*> elements to zero, and for each unreduced block T_i, it finds +*> (a) a suitable shift at one end of the block's spectrum, +*> (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and +*> (c) eigenvalues of each L_i D_i L_i^T. +*> The representations and eigenvalues found are then used by +*> DSTEMR to compute the eigenvectors of T. +*> The accuracy varies depending on whether bisection is used to +*> find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to +*> conpute all and then discard any unwanted one. +*> As an added benefit, DLARRE also outputs the n +*> Gerschgorin intervals for the matrices L_i D_i L_i^T. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] RANGE +*> \verbatim +*> RANGE is CHARACTER*1 +*> = 'A': ("All") all eigenvalues will be found. +*> = 'V': ("Value") all eigenvalues in the half-open interval +*> (VL, VU] will be found. +*> = 'I': ("Index") the IL-th through IU-th eigenvalues (of the +*> entire matrix) will be found. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix. N > 0. +*> \endverbatim +*> +*> \param[in,out] VL +*> \verbatim +*> VL is DOUBLE PRECISION +*> If RANGE='V', the lower bound for the eigenvalues. +*> Eigenvalues less than or equal to VL, or greater than VU, +*> will not be returned. VL < VU. +*> If RANGE='I' or ='A', DLARRE computes bounds on the desired +*> part of the spectrum. +*> \endverbatim +*> +*> \param[in,out] VU +*> \verbatim +*> VU is DOUBLE PRECISION +*> If RANGE='V', the upper bound for the eigenvalues. +*> Eigenvalues less than or equal to VL, or greater than VU, +*> will not be returned. VL < VU. +*> If RANGE='I' or ='A', DLARRE computes bounds on the desired +*> part of the spectrum. +*> \endverbatim +*> +*> \param[in] IL +*> \verbatim +*> IL is INTEGER +*> If RANGE='I', the index of the +*> smallest eigenvalue to be returned. +*> 1 <= IL <= IU <= N. +*> \endverbatim +*> +*> \param[in] IU +*> \verbatim +*> IU is INTEGER +*> If RANGE='I', the index of the +*> largest eigenvalue to be returned. +*> 1 <= IL <= IU <= N. +*> \endverbatim +*> +*> \param[in,out] D +*> \verbatim +*> D is DOUBLE PRECISION array, dimension (N) +*> On entry, the N diagonal elements of the tridiagonal +*> matrix T. +*> On exit, the N diagonal elements of the diagonal +*> matrices D_i. +*> \endverbatim +*> +*> \param[in,out] E +*> \verbatim +*> E is DOUBLE PRECISION array, dimension (N) +*> On entry, the first (N-1) entries contain the subdiagonal +*> elements of the tridiagonal matrix T; E(N) need not be set. +*> On exit, E contains the subdiagonal elements of the unit +*> bidiagonal matrices L_i. The entries E( ISPLIT( I ) ), +*> 1 <= I <= NSPLIT, contain the base points sigma_i on output. +*> \endverbatim +*> +*> \param[in,out] E2 +*> \verbatim +*> E2 is DOUBLE PRECISION array, dimension (N) +*> On entry, the first (N-1) entries contain the SQUARES of the +*> subdiagonal elements of the tridiagonal matrix T; +*> E2(N) need not be set. +*> On exit, the entries E2( ISPLIT( I ) ), +*> 1 <= I <= NSPLIT, have been set to zero +*> \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] SPLTOL +*> \verbatim +*> SPLTOL is DOUBLE PRECISION +*> The threshold for splitting. +*> \endverbatim +*> +*> \param[out] NSPLIT +*> \verbatim +*> NSPLIT is INTEGER +*> The number of blocks T splits into. 1 <= NSPLIT <= N. +*> \endverbatim +*> +*> \param[out] 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., and the NSPLIT-th consists of rows/columns +*> ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. +*> \endverbatim +*> +*> \param[out] M +*> \verbatim +*> M is INTEGER +*> The total number of eigenvalues (of all L_i D_i L_i^T) +*> found. +*> \endverbatim +*> +*> \param[out] W +*> \verbatim +*> W is DOUBLE PRECISION array, dimension (N) +*> The first M elements contain the eigenvalues. The +*> eigenvalues of each of the blocks, L_i D_i L_i^T, are +*> sorted in ascending order ( DLARRE may use the +*> remaining N-M elements as workspace). +*> \endverbatim +*> +*> \param[out] WERR +*> \verbatim +*> WERR is DOUBLE PRECISION array, dimension (N) +*> The error bound on the corresponding eigenvalue in W. +*> \endverbatim +*> +*> \param[out] WGAP +*> \verbatim +*> WGAP is DOUBLE PRECISION array, dimension (N) +*> The separation from the right neighbor eigenvalue in W. +*> The gap is only with respect to the eigenvalues of the same block +*> as each block has its own representation tree. +*> Exception: at the right end of a block we store the left gap +*> \endverbatim +*> +*> \param[out] 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[out] 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 block 2 +*> \endverbatim +*> +*> \param[out] 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)). +*> \endverbatim +*> +*> \param[out] PIVMIN +*> \verbatim +*> PIVMIN is DOUBLE PRECISION +*> The minimum pivot in the Sturm sequence for T. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (6*N) +*> Workspace. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (5*N) +*> Workspace. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> > 0: A problem occurred in DLARRE. +*> < 0: One of the called subroutines signaled an internal problem. +*> Needs inspection of the corresponding parameter IINFO +*> for further information. +*> +*> =-1: Problem in DLARRD. +*> = 2: No base representation could be found in MAXTRY iterations. +*> Increasing MAXTRY and recompilation might be a remedy. +*> =-3: Problem in DLARRB when computing the refined root +*> representation for DLASQ2. +*> =-4: Problem in DLARRB when preforming bisection on the +*> desired part of the spectrum. +*> =-5: Problem in DLASQ2. +*> =-6: Problem in DLASQ2. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \ingroup OTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The base representations are required to suffer very little +*> element growth and consequently define all their eigenvalues to +*> high relative accuracy. +*> \endverbatim +* +*> \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 \n +*> +* ===================================================================== SUBROUTINE DLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2, $ RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, $ W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, $ WORK, IWORK, INFO ) - IMPLICIT NONE * -* -- 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 .. CHARACTER RANGE @@ -21,162 +319,6 @@ $ W( * ),WERR( * ), WGAP( * ), WORK( * ) * .. * -* Purpose -* ======= -* -* To find the desired eigenvalues of a given real symmetric -* tridiagonal matrix T, DLARRE sets any "small" off-diagonal -* elements to zero, and for each unreduced block T_i, it finds -* (a) a suitable shift at one end of the block's spectrum, -* (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and -* (c) eigenvalues of each L_i D_i L_i^T. -* The representations and eigenvalues found are then used by -* DSTEMR to compute the eigenvectors of T. -* The accuracy varies depending on whether bisection is used to -* find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to -* conpute all and then discard any unwanted one. -* As an added benefit, DLARRE also outputs the n -* Gerschgorin intervals for the matrices L_i D_i L_i^T. -* -* Arguments -* ========= -* -* RANGE (input) CHARACTER -* = 'A': ("All") all eigenvalues will be found. -* = 'V': ("Value") all eigenvalues in the half-open interval -* (VL, VU] will be found. -* = 'I': ("Index") the IL-th through IU-th eigenvalues (of the -* entire matrix) will be found. -* -* N (input) INTEGER -* The order of the matrix. N > 0. -* -* VL (input/output) DOUBLE PRECISION -* VU (input/output) DOUBLE PRECISION -* If RANGE='V', the lower and upper bounds for the eigenvalues. -* Eigenvalues less than or equal to VL, or greater than VU, -* will not be returned. VL < VU. -* If RANGE='I' or ='A', DLARRE computes bounds on the desired -* part of the spectrum. -* -* IL (input) INTEGER -* IU (input) INTEGER -* If RANGE='I', the indices (in ascending order) of the -* smallest and largest eigenvalues to be returned. -* 1 <= IL <= IU <= N. -* -* D (input/output) DOUBLE PRECISION array, dimension (N) -* On entry, the N diagonal elements of the tridiagonal -* matrix T. -* On exit, the N diagonal elements of the diagonal -* matrices D_i. -* -* E (input/output) DOUBLE PRECISION array, dimension (N) -* On entry, the first (N-1) entries contain the subdiagonal -* elements of the tridiagonal matrix T; E(N) need not be set. -* On exit, E contains the subdiagonal elements of the unit -* bidiagonal matrices L_i. The entries E( ISPLIT( I ) ), -* 1 <= I <= NSPLIT, contain the base points sigma_i on output. -* -* E2 (input/output) DOUBLE PRECISION array, dimension (N) -* On entry, the first (N-1) entries contain the SQUARES of the -* subdiagonal elements of the tridiagonal matrix T; -* E2(N) need not be set. -* On exit, the entries E2( ISPLIT( I ) ), -* 1 <= I <= NSPLIT, have been set to zero -* -* 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|) ) -* -* SPLTOL (input) DOUBLE PRECISION -* The threshold for splitting. -* -* NSPLIT (output) INTEGER -* The number of blocks T splits into. 1 <= NSPLIT <= N. -* -* ISPLIT (output) 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., and the NSPLIT-th consists of rows/columns -* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. -* -* M (output) INTEGER -* The total number of eigenvalues (of all L_i D_i L_i^T) -* found. -* -* W (output) DOUBLE PRECISION array, dimension (N) -* The first M elements contain the eigenvalues. The -* eigenvalues of each of the blocks, L_i D_i L_i^T, are -* sorted in ascending order ( DLARRE may use the -* remaining N-M elements as workspace). -* -* WERR (output) DOUBLE PRECISION array, dimension (N) -* The error bound on the corresponding eigenvalue in W. -* -* WGAP (output) DOUBLE PRECISION array, dimension (N) -* The separation from the right neighbor eigenvalue in W. -* The gap is only with respect to the eigenvalues of the same block -* as each block has its own representation tree. -* Exception: at the right end of a block we store the left gap -* -* IBLOCK (output) 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 (output) 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 block 2 -* -* GERS (output) DOUBLE PRECISION array, dimension (2*N) -* The N Gerschgorin intervals (the i-th Gerschgorin interval -* is (GERS(2*i-1), GERS(2*i)). -* -* PIVMIN (output) DOUBLE PRECISION -* The minimum pivot in the Sturm sequence for T. -* -* WORK (workspace) DOUBLE PRECISION array, dimension (6*N) -* Workspace. -* -* IWORK (workspace) INTEGER array, dimension (5*N) -* Workspace. -* -* INFO (output) INTEGER -* = 0: successful exit -* > 0: A problem occured in DLARRE. -* < 0: One of the called subroutines signaled an internal problem. -* Needs inspection of the corresponding parameter IINFO -* for further information. -* -* =-1: Problem in DLARRD. -* = 2: No base representation could be found in MAXTRY iterations. -* Increasing MAXTRY and recompilation might be a remedy. -* =-3: Problem in DLARRB when computing the refined root -* representation for DLASQ2. -* =-4: Problem in DLARRB when preforming bisection on the -* desired part of the spectrum. -* =-5: Problem in DLASQ2. -* =-6: Problem in DLASQ2. -* -* Further Details -* The base representations are required to suffer very little -* element growth and consequently define all their eigenvalues to -* high relative accuracy. -* =============== -* -* 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 .. @@ -215,7 +357,7 @@ * .. * .. External Subroutines .. EXTERNAL DCOPY, DLARNV, DLARRA, DLARRB, DLARRC, DLARRD, - $ DLASQ2 + $ DLASQ2, DLARRK * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN @@ -225,7 +367,12 @@ * INFO = 0 - +* +* Quick return if possible +* + IF( N.LE.0 ) THEN + RETURN + END IF * * Decode RANGE * @@ -532,6 +679,7 @@ * The initial SIGMA was to the outer end of the spectrum * the matrix is definite and we need not retreat. TAU = SPDIAM*EPS*N + TWO*PIVMIN + TAU = MAX( TAU,TWO*EPS*ABS(SIGMA) ) ELSE IF(MB.GT.1) THEN CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN) @@ -748,6 +896,6 @@ RETURN * -* end of DLARRE +* End of DLARRE * END