version 1.6, 2010/08/13 21:03:58
|
version 1.13, 2014/01/27 09:28:27
|
Line 1
|
Line 1
|
|
*> \brief \b DSTEBZ |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download DSTEBZ + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstebz.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstebz.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstebz.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, |
|
* M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, |
|
* INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* CHARACTER ORDER, RANGE |
|
* INTEGER IL, INFO, IU, M, N, NSPLIT |
|
* DOUBLE PRECISION ABSTOL, VL, VU |
|
* .. |
|
* .. Array Arguments .. |
|
* INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ) |
|
* DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> DSTEBZ computes the eigenvalues of a symmetric tridiagonal |
|
*> matrix T. The user may ask for all eigenvalues, all eigenvalues |
|
*> in the half-open interval (VL, VU], or the IL-th through IU-th |
|
*> eigenvalues. |
|
*> |
|
*> To avoid overflow, the matrix must be scaled so that its |
|
*> largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest |
|
*> accuracy, it should not be much smaller than that. |
|
*> |
|
*> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal |
|
*> Matrix", Report CS41, Computer Science Dept., Stanford |
|
*> University, July 21, 1966. |
|
*> \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] ORDER |
|
*> \verbatim |
|
*> ORDER is CHARACTER*1 |
|
*> = 'B': ("By Block") the eigenvalues will be grouped by |
|
*> split-off block (see IBLOCK, ISPLIT) and |
|
*> ordered from smallest to largest within |
|
*> the block. |
|
*> = 'E': ("Entire matrix") |
|
*> the eigenvalues for the entire matrix |
|
*> will be ordered from smallest to |
|
*> largest. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] N |
|
*> \verbatim |
|
*> N is INTEGER |
|
*> The order of the tridiagonal matrix T. N >= 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] VL |
|
*> \verbatim |
|
*> VL is DOUBLE PRECISION |
|
*> \endverbatim |
|
*> |
|
*> \param[in] VU |
|
*> \verbatim |
|
*> VU is DOUBLE PRECISION |
|
*> |
|
*> If RANGE='V', the lower and upper bounds of the interval to |
|
*> be searched for eigenvalues. Eigenvalues less than or equal |
|
*> to VL, or greater than VU, will not be returned. VL < VU. |
|
*> Not referenced if RANGE = 'A' or 'I'. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] IL |
|
*> \verbatim |
|
*> IL is INTEGER |
|
*> \endverbatim |
|
*> |
|
*> \param[in] IU |
|
*> \verbatim |
|
*> IU is INTEGER |
|
*> |
|
*> If RANGE='I', the indices (in ascending order) of the |
|
*> smallest and largest eigenvalues to be returned. |
|
*> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. |
|
*> Not referenced if RANGE = 'A' or 'V'. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] ABSTOL |
|
*> \verbatim |
|
*> ABSTOL is DOUBLE PRECISION |
|
*> The absolute tolerance for the eigenvalues. An eigenvalue |
|
*> (or cluster) is considered to be located if it has been |
|
*> determined to lie in an interval whose width is ABSTOL or |
|
*> less. If ABSTOL is less than or equal to zero, then ULP*|T| |
|
*> will be used, where |T| means the 1-norm of T. |
|
*> |
|
*> Eigenvalues will be computed most accurately when ABSTOL is |
|
*> set to twice the underflow threshold 2*DLAMCH('S'), not zero. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] D |
|
*> \verbatim |
|
*> D is DOUBLE PRECISION array, dimension (N) |
|
*> The n diagonal elements of the tridiagonal matrix T. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] E |
|
*> \verbatim |
|
*> E is DOUBLE PRECISION array, dimension (N-1) |
|
*> The (n-1) off-diagonal elements of the tridiagonal matrix T. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] M |
|
*> \verbatim |
|
*> M is INTEGER |
|
*> The actual number of eigenvalues found. 0 <= M <= N. |
|
*> (See also the description of INFO=2,3.) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] NSPLIT |
|
*> \verbatim |
|
*> NSPLIT is INTEGER |
|
*> The number of diagonal blocks in the matrix T. |
|
*> 1 <= NSPLIT <= N. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] W |
|
*> \verbatim |
|
*> W is DOUBLE PRECISION array, dimension (N) |
|
*> On exit, the first M elements of W will contain the |
|
*> eigenvalues. (DSTEBZ may use the remaining N-M elements as |
|
*> workspace.) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] IBLOCK |
|
*> \verbatim |
|
*> IBLOCK is INTEGER array, dimension (N) |
|
*> At each row/column j where E(j) is zero or small, the |
|
*> matrix T is considered to split into a block diagonal |
|
*> matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which |
|
*> block (from 1 to the number of blocks) the eigenvalue W(i) |
|
*> belongs. (DSTEBZ may use the remaining N-M elements as |
|
*> workspace.) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] ISPLIT |
|
*> \verbatim |
|
*> ISPLIT is INTEGER array, dimension (N) |
|
*> The splitting points, at which T breaks up into submatrices. |
|
*> The first submatrix 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. |
|
*> (Only the first NSPLIT elements will actually be used, but |
|
*> since the user cannot know a priori what value NSPLIT will |
|
*> have, N words must be reserved for ISPLIT.) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] WORK |
|
*> \verbatim |
|
*> WORK is DOUBLE PRECISION array, dimension (4*N) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] IWORK |
|
*> \verbatim |
|
*> IWORK is INTEGER array, dimension (3*N) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] INFO |
|
*> \verbatim |
|
*> INFO is INTEGER |
|
*> = 0: successful exit |
|
*> < 0: if INFO = -i, the i-th argument had an illegal value |
|
*> > 0: some or all of the eigenvalues failed to converge or |
|
*> were not computed: |
|
*> =1 or 3: Bisection failed to converge for some |
|
*> eigenvalues; these eigenvalues are flagged by a |
|
*> negative block number. The effect is that the |
|
*> eigenvalues may not be as accurate as the |
|
*> absolute and relative tolerances. This is |
|
*> generally caused by unexpectedly inaccurate |
|
*> arithmetic. |
|
*> =2 or 3: RANGE='I' only: Not all of the eigenvalues |
|
*> IL:IU were found. |
|
*> Effect: M < IU+1-IL |
|
*> Cause: non-monotonic arithmetic, causing the |
|
*> Sturm sequence to be non-monotonic. |
|
*> Cure: recalculate, using RANGE='A', and pick |
|
*> out eigenvalues IL:IU. In some cases, |
|
*> increasing the PARAMETER "FUDGE" may |
|
*> make things work. |
|
*> = 4: RANGE='I', and the Gershgorin interval |
|
*> initially used was too small. No eigenvalues |
|
*> were computed. |
|
*> Probable cause: your machine has sloppy |
|
*> floating-point arithmetic. |
|
*> Cure: Increase the PARAMETER "FUDGE", |
|
*> recompile, and try again. |
|
*> \endverbatim |
|
* |
|
*> \par Internal Parameters: |
|
* ========================= |
|
*> |
|
*> \verbatim |
|
*> RELFAC DOUBLE PRECISION, default = 2.0e0 |
|
*> The relative tolerance. An interval (a,b] lies within |
|
*> "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|), |
|
*> where "ulp" is the machine precision (distance from 1 to |
|
*> the next larger floating point number.) |
|
*> |
|
*> FUDGE DOUBLE PRECISION, default = 2 |
|
*> A "fudge factor" to widen the Gershgorin intervals. Ideally, |
|
*> a value of 1 should work, but on machines with sloppy |
|
*> arithmetic, this needs to be larger. The default for |
|
*> publicly released versions should be large enough to handle |
|
*> the worst machine around. Note that this has no effect |
|
*> on accuracy of the solution. |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \date November 2011 |
|
* |
|
*> \ingroup auxOTHERcomputational |
|
* |
|
* ===================================================================== |
SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, |
SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, |
$ M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, |
$ M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, |
$ INFO ) |
$ INFO ) |
* |
* |
* -- LAPACK routine (version 3.2) -- |
* -- LAPACK computational routine (version 3.4.0) -- |
* -- 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 |
* November 2011 |
* 8-18-00: Increase FUDGE factor for T3E (eca) |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER ORDER, RANGE |
CHARACTER ORDER, RANGE |
Line 18
|
Line 278
|
DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) |
DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* DSTEBZ computes the eigenvalues of a symmetric tridiagonal |
|
* matrix T. The user may ask for all eigenvalues, all eigenvalues |
|
* in the half-open interval (VL, VU], or the IL-th through IU-th |
|
* eigenvalues. |
|
* |
|
* To avoid overflow, the matrix must be scaled so that its |
|
* largest element is no greater than overflow**(1/2) * |
|
* underflow**(1/4) in absolute value, and for greatest |
|
* accuracy, it should not be much smaller than that. |
|
* |
|
* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal |
|
* Matrix", Report CS41, Computer Science Dept., Stanford |
|
* University, July 21, 1966. |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* RANGE (input) 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. |
|
* |
|
* ORDER (input) CHARACTER*1 |
|
* = 'B': ("By Block") the eigenvalues will be grouped by |
|
* split-off block (see IBLOCK, ISPLIT) and |
|
* ordered from smallest to largest within |
|
* the block. |
|
* = 'E': ("Entire matrix") |
|
* the eigenvalues for the entire matrix |
|
* will be ordered from smallest to |
|
* largest. |
|
* |
|
* N (input) INTEGER |
|
* The order of the tridiagonal matrix T. N >= 0. |
|
* |
|
* VL (input) DOUBLE PRECISION |
|
* VU (input) DOUBLE PRECISION |
|
* If RANGE='V', the lower and upper bounds of the interval to |
|
* be searched for eigenvalues. Eigenvalues less than or equal |
|
* to VL, or greater than VU, will not be returned. VL < VU. |
|
* Not referenced if RANGE = 'A' or 'I'. |
|
* |
|
* 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, if N > 0; IL = 1 and IU = 0 if N = 0. |
|
* Not referenced if RANGE = 'A' or 'V'. |
|
* |
|
* ABSTOL (input) DOUBLE PRECISION |
|
* The absolute tolerance for the eigenvalues. An eigenvalue |
|
* (or cluster) is considered to be located if it has been |
|
* determined to lie in an interval whose width is ABSTOL or |
|
* less. If ABSTOL is less than or equal to zero, then ULP*|T| |
|
* will be used, where |T| means the 1-norm of T. |
|
* |
|
* Eigenvalues will be computed most accurately when ABSTOL is |
|
* set to twice the underflow threshold 2*DLAMCH('S'), not zero. |
|
* |
|
* D (input) DOUBLE PRECISION array, dimension (N) |
|
* The n diagonal elements of the tridiagonal matrix T. |
|
* |
|
* E (input) DOUBLE PRECISION array, dimension (N-1) |
|
* The (n-1) off-diagonal elements of the tridiagonal matrix T. |
|
* |
|
* M (output) INTEGER |
|
* The actual number of eigenvalues found. 0 <= M <= N. |
|
* (See also the description of INFO=2,3.) |
|
* |
|
* NSPLIT (output) INTEGER |
|
* The number of diagonal blocks in the matrix T. |
|
* 1 <= NSPLIT <= N. |
|
* |
|
* W (output) DOUBLE PRECISION array, dimension (N) |
|
* On exit, the first M elements of W will contain the |
|
* eigenvalues. (DSTEBZ may use the remaining N-M elements as |
|
* workspace.) |
|
* |
|
* IBLOCK (output) INTEGER array, dimension (N) |
|
* At each row/column j where E(j) is zero or small, the |
|
* matrix T is considered to split into a block diagonal |
|
* matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which |
|
* block (from 1 to the number of blocks) the eigenvalue W(i) |
|
* belongs. (DSTEBZ may use the remaining N-M elements as |
|
* workspace.) |
|
* |
|
* ISPLIT (output) INTEGER array, dimension (N) |
|
* The splitting points, at which T breaks up into submatrices. |
|
* The first submatrix 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. |
|
* (Only the first NSPLIT elements will actually be used, but |
|
* since the user cannot know a priori what value NSPLIT will |
|
* have, N words must be reserved for ISPLIT.) |
|
* |
|
* WORK (workspace) DOUBLE PRECISION array, dimension (4*N) |
|
* |
|
* IWORK (workspace) INTEGER array, dimension (3*N) |
|
* |
|
* INFO (output) INTEGER |
|
* = 0: successful exit |
|
* < 0: if INFO = -i, the i-th argument had an illegal value |
|
* > 0: some or all of the eigenvalues failed to converge or |
|
* were not computed: |
|
* =1 or 3: Bisection failed to converge for some |
|
* eigenvalues; these eigenvalues are flagged by a |
|
* negative block number. The effect is that the |
|
* eigenvalues may not be as accurate as the |
|
* absolute and relative tolerances. This is |
|
* generally caused by unexpectedly inaccurate |
|
* arithmetic. |
|
* =2 or 3: RANGE='I' only: Not all of the eigenvalues |
|
* IL:IU were found. |
|
* Effect: M < IU+1-IL |
|
* Cause: non-monotonic arithmetic, causing the |
|
* Sturm sequence to be non-monotonic. |
|
* Cure: recalculate, using RANGE='A', and pick |
|
* out eigenvalues IL:IU. In some cases, |
|
* increasing the PARAMETER "FUDGE" may |
|
* make things work. |
|
* = 4: RANGE='I', and the Gershgorin interval |
|
* initially used was too small. No eigenvalues |
|
* were computed. |
|
* Probable cause: your machine has sloppy |
|
* floating-point arithmetic. |
|
* Cure: Increase the PARAMETER "FUDGE", |
|
* recompile, and try again. |
|
* |
|
* Internal Parameters |
|
* =================== |
|
* |
|
* RELFAC DOUBLE PRECISION, default = 2.0e0 |
|
* The relative tolerance. An interval (a,b] lies within |
|
* "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|), |
|
* where "ulp" is the machine precision (distance from 1 to |
|
* the next larger floating point number.) |
|
* |
|
* FUDGE DOUBLE PRECISION, default = 2 |
|
* A "fudge factor" to widen the Gershgorin intervals. Ideally, |
|
* a value of 1 should work, but on machines with sloppy |
|
* arithmetic, this needs to be larger. The default for |
|
* publicly released versions should be large enough to handle |
|
* the worst machine around. Note that this has no effect |
|
* on accuracy of the solution. |
|
* |
|
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. Parameters .. |
* .. Parameters .. |
Line 301
|
Line 410
|
WORK( N ) = ZERO |
WORK( N ) = ZERO |
PIVMIN = ONE |
PIVMIN = ONE |
* |
* |
*DIR$ NOVECTOR |
|
DO 10 J = 2, N |
DO 10 J = 2, N |
TMP1 = E( J-1 )**2 |
TMP1 = E( J-1 )**2 |
IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN |
IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN |