version 1.6, 2010/08/07 13:22:19
|
version 1.24, 2023/08/07 08:38:57
|
Line 1
|
Line 1
|
|
*> \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 |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarre.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarre.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarre.f"> |
|
*> [TXT]</a> |
|
*> \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, |
SUBROUTINE DLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2, |
$ RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, |
$ RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, |
$ W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, |
$ W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, |
$ WORK, IWORK, INFO ) |
$ 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, -- |
* -- 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..-- |
* June 2010 |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER RANGE |
CHARACTER RANGE |
Line 21
|
Line 319
|
$ W( * ),WERR( * ), WGAP( * ), WORK( * ) |
$ 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 .. |
* .. Parameters .. |
Line 215
|
Line 357
|
* .. |
* .. |
* .. External Subroutines .. |
* .. External Subroutines .. |
EXTERNAL DCOPY, DLARNV, DLARRA, DLARRB, DLARRC, DLARRD, |
EXTERNAL DCOPY, DLARNV, DLARRA, DLARRB, DLARRC, DLARRD, |
$ DLASQ2 |
$ DLASQ2, DLARRK |
* .. |
* .. |
* .. Intrinsic Functions .. |
* .. Intrinsic Functions .. |
INTRINSIC ABS, MAX, MIN |
INTRINSIC ABS, MAX, MIN |
Line 225
|
Line 367
|
* |
* |
|
|
INFO = 0 |
INFO = 0 |
|
* |
|
* Quick return if possible |
|
* |
|
IF( N.LE.0 ) THEN |
|
RETURN |
|
END IF |
* |
* |
* Decode RANGE |
* Decode RANGE |
* |
* |
Line 532
|
Line 679
|
* The initial SIGMA was to the outer end of the spectrum |
* The initial SIGMA was to the outer end of the spectrum |
* the matrix is definite and we need not retreat. |
* the matrix is definite and we need not retreat. |
TAU = SPDIAM*EPS*N + TWO*PIVMIN |
TAU = SPDIAM*EPS*N + TWO*PIVMIN |
|
TAU = MAX( TAU,TWO*EPS*ABS(SIGMA) ) |
ELSE |
ELSE |
IF(MB.GT.1) THEN |
IF(MB.GT.1) THEN |
CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN) |
CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN) |
Line 748
|
Line 896
|
|
|
RETURN |
RETURN |
* |
* |
* end of DLARRE |
* End of DLARRE |
* |
* |
END |
END |