--- rpl/lapack/lapack/dlarre.f 2010/08/13 21:03:51 1.7
+++ 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.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..--
-* June 2010
*
* .. 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