dstebz.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 677 行 · 第 1/3 页
HTML
677 行
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<title>dstebz.f</title>
<meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!--
body { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.default { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.default a { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
span.string { color: rgb(188, 143, 143); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.string a { color: rgb(188, 143, 143); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
span.comment { color: rgb(178, 34, 34); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.comment a { color: rgb(178, 34, 34); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
--></style>
</head>
<body>
<pre>
SUBROUTINE <a name="DSTEBZ.1"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a>( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
$ M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
$ INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment"> Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment"> November 2006
</span><span class="comment">*</span><span class="comment"> 8-18-00: Increase FUDGE factor for T3E (eca)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Scalar Arguments ..
</span> CHARACTER ORDER, RANGE
INTEGER IL, INFO, IU, M, N, NSPLIT
DOUBLE PRECISION ABSTOL, VL, VU
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Array Arguments ..
</span> INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * )
DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Purpose
</span><span class="comment">*</span><span class="comment"> =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> <a name="DSTEBZ.23"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> computes the eigenvalues of a symmetric tridiagonal
</span><span class="comment">*</span><span class="comment"> matrix T. The user may ask for all eigenvalues, all eigenvalues
</span><span class="comment">*</span><span class="comment"> in the half-open interval (VL, VU], or the IL-th through IU-th
</span><span class="comment">*</span><span class="comment"> eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> To avoid overflow, the matrix must be scaled so that its
</span><span class="comment">*</span><span class="comment"> largest element is no greater than overflow**(1/2) *
</span><span class="comment">*</span><span class="comment"> underflow**(1/4) in absolute value, and for greatest
</span><span class="comment">*</span><span class="comment"> accuracy, it should not be much smaller than that.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
</span><span class="comment">*</span><span class="comment"> Matrix", Report CS41, Computer Science Dept., Stanford
</span><span class="comment">*</span><span class="comment"> University, July 21, 1966.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Arguments
</span><span class="comment">*</span><span class="comment"> =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> RANGE (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment"> = 'A': ("All") all eigenvalues will be found.
</span><span class="comment">*</span><span class="comment"> = 'V': ("Value") all eigenvalues in the half-open interval
</span><span class="comment">*</span><span class="comment"> (VL, VU] will be found.
</span><span class="comment">*</span><span class="comment"> = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
</span><span class="comment">*</span><span class="comment"> entire matrix) will be found.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ORDER (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment"> = 'B': ("By Block") the eigenvalues will be grouped by
</span><span class="comment">*</span><span class="comment"> split-off block (see IBLOCK, ISPLIT) and
</span><span class="comment">*</span><span class="comment"> ordered from smallest to largest within
</span><span class="comment">*</span><span class="comment"> the block.
</span><span class="comment">*</span><span class="comment"> = 'E': ("Entire matrix")
</span><span class="comment">*</span><span class="comment"> the eigenvalues for the entire matrix
</span><span class="comment">*</span><span class="comment"> will be ordered from smallest to
</span><span class="comment">*</span><span class="comment"> largest.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> N (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The order of the tridiagonal matrix T. N >= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> VL (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment"> VU (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment"> If RANGE='V', the lower and upper bounds of the interval to
</span><span class="comment">*</span><span class="comment"> be searched for eigenvalues. Eigenvalues less than or equal
</span><span class="comment">*</span><span class="comment"> to VL, or greater than VU, will not be returned. VL < VU.
</span><span class="comment">*</span><span class="comment"> Not referenced if RANGE = 'A' or 'I'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> IL (input) INTEGER
</span><span class="comment">*</span><span class="comment"> IU (input) INTEGER
</span><span class="comment">*</span><span class="comment"> If RANGE='I', the indices (in ascending order) of the
</span><span class="comment">*</span><span class="comment"> smallest and largest eigenvalues to be returned.
</span><span class="comment">*</span><span class="comment"> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
</span><span class="comment">*</span><span class="comment"> Not referenced if RANGE = 'A' or 'V'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ABSTOL (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment"> The absolute tolerance for the eigenvalues. An eigenvalue
</span><span class="comment">*</span><span class="comment"> (or cluster) is considered to be located if it has been
</span><span class="comment">*</span><span class="comment"> determined to lie in an interval whose width is ABSTOL or
</span><span class="comment">*</span><span class="comment"> less. If ABSTOL is less than or equal to zero, then ULP*|T|
</span><span class="comment">*</span><span class="comment"> will be used, where |T| means the 1-norm of T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Eigenvalues will be computed most accurately when ABSTOL is
</span><span class="comment">*</span><span class="comment"> set to twice the underflow threshold 2*<a name="DLAMCH.82"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>('S'), not zero.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> D (input) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment"> The n diagonal elements of the tridiagonal matrix T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> E (input) DOUBLE PRECISION array, dimension (N-1)
</span><span class="comment">*</span><span class="comment"> The (n-1) off-diagonal elements of the tridiagonal matrix T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> M (output) INTEGER
</span><span class="comment">*</span><span class="comment"> The actual number of eigenvalues found. 0 <= M <= N.
</span><span class="comment">*</span><span class="comment"> (See also the description of INFO=2,3.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> NSPLIT (output) INTEGER
</span><span class="comment">*</span><span class="comment"> The number of diagonal blocks in the matrix T.
</span><span class="comment">*</span><span class="comment"> 1 <= NSPLIT <= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> W (output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment"> On exit, the first M elements of W will contain the
</span><span class="comment">*</span><span class="comment"> eigenvalues. (<a name="DSTEBZ.100"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> may use the remaining N-M elements as
</span><span class="comment">*</span><span class="comment"> workspace.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> IBLOCK (output) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment"> At each row/column j where E(j) is zero or small, the
</span><span class="comment">*</span><span class="comment"> matrix T is considered to split into a block diagonal
</span><span class="comment">*</span><span class="comment"> matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which
</span><span class="comment">*</span><span class="comment"> block (from 1 to the number of blocks) the eigenvalue W(i)
</span><span class="comment">*</span><span class="comment"> belongs. (<a name="DSTEBZ.108"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> may use the remaining N-M elements as
</span><span class="comment">*</span><span class="comment"> workspace.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ISPLIT (output) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment"> The splitting points, at which T breaks up into submatrices.
</span><span class="comment">*</span><span class="comment"> The first submatrix consists of rows/columns 1 to ISPLIT(1),
</span><span class="comment">*</span><span class="comment"> the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
</span><span class="comment">*</span><span class="comment"> etc., and the NSPLIT-th consists of rows/columns
</span><span class="comment">*</span><span class="comment"> ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
</span><span class="comment">*</span><span class="comment"> (Only the first NSPLIT elements will actually be used, but
</span><span class="comment">*</span><span class="comment"> since the user cannot know a priori what value NSPLIT will
</span><span class="comment">*</span><span class="comment"> have, N words must be reserved for ISPLIT.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> IWORK (workspace) INTEGER array, dimension (3*N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> INFO (output) INTEGER
</span><span class="comment">*</span><span class="comment"> = 0: successful exit
</span><span class="comment">*</span><span class="comment"> < 0: if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment"> > 0: some or all of the eigenvalues failed to converge or
</span><span class="comment">*</span><span class="comment"> were not computed:
</span><span class="comment">*</span><span class="comment"> =1 or 3: Bisection failed to converge for some
</span><span class="comment">*</span><span class="comment"> eigenvalues; these eigenvalues are flagged by a
</span><span class="comment">*</span><span class="comment"> negative block number. The effect is that the
</span><span class="comment">*</span><span class="comment"> eigenvalues may not be as accurate as the
</span><span class="comment">*</span><span class="comment"> absolute and relative tolerances. This is
</span><span class="comment">*</span><span class="comment"> generally caused by unexpectedly inaccurate
</span><span class="comment">*</span><span class="comment"> arithmetic.
</span><span class="comment">*</span><span class="comment"> =2 or 3: RANGE='I' only: Not all of the eigenvalues
</span><span class="comment">*</span><span class="comment"> IL:IU were found.
</span><span class="comment">*</span><span class="comment"> Effect: M < IU+1-IL
</span><span class="comment">*</span><span class="comment"> Cause: non-monotonic arithmetic, causing the
</span><span class="comment">*</span><span class="comment"> Sturm sequence to be non-monotonic.
</span><span class="comment">*</span><span class="comment"> Cure: recalculate, using RANGE='A', and pick
</span><span class="comment">*</span><span class="comment"> out eigenvalues IL:IU. In some cases,
</span><span class="comment">*</span><span class="comment"> increasing the PARAMETER "FUDGE" may
</span><span class="comment">*</span><span class="comment"> make things work.
</span><span class="comment">*</span><span class="comment"> = 4: RANGE='I', and the Gershgorin interval
</span><span class="comment">*</span><span class="comment"> initially used was too small. No eigenvalues
</span><span class="comment">*</span><span class="comment"> were computed.
</span><span class="comment">*</span><span class="comment"> Probable cause: your machine has sloppy
</span><span class="comment">*</span><span class="comment"> floating-point arithmetic.
</span><span class="comment">*</span><span class="comment"> Cure: Increase the PARAMETER "FUDGE",
</span><span class="comment">*</span><span class="comment"> recompile, and try again.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Internal Parameters
</span><span class="comment">*</span><span class="comment"> ===================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> RELFAC DOUBLE PRECISION, default = 2.0e0
</span><span class="comment">*</span><span class="comment"> The relative tolerance. An interval (a,b] lies within
</span><span class="comment">*</span><span class="comment"> "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|),
</span><span class="comment">*</span><span class="comment"> where "ulp" is the machine precision (distance from 1 to
</span><span class="comment">*</span><span class="comment"> the next larger floating point number.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> FUDGE DOUBLE PRECISION, default = 2
</span><span class="comment">*</span><span class="comment"> A "fudge factor" to widen the Gershgorin intervals. Ideally,
</span><span class="comment">*</span><span class="comment"> a value of 1 should work, but on machines with sloppy
</span><span class="comment">*</span><span class="comment"> arithmetic, this needs to be larger. The default for
</span><span class="comment">*</span><span class="comment"> publicly released versions should be large enough to handle
</span><span class="comment">*</span><span class="comment"> the worst machine around. Note that this has no effect
</span><span class="comment">*</span><span class="comment"> on accuracy of the solution.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Parameters ..
</span> DOUBLE PRECISION ZERO, ONE, TWO, HALF
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
$ HALF = 1.0D0 / TWO )
DOUBLE PRECISION FUDGE, RELFAC
PARAMETER ( FUDGE = 2.1D0, RELFAC = 2.0D0 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Scalars ..
</span> LOGICAL NCNVRG, TOOFEW
INTEGER IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
$ IM, IN, IOFF, IORDER, IOUT, IRANGE, ITMAX,
$ ITMP1, IW, IWOFF, J, JB, JDISC, JE, NB, NWL,
$ NWU
DOUBLE PRECISION ATOLI, BNORM, GL, GU, PIVMIN, RTOLI, SAFEMN,
$ TMP1, TMP2, TNORM, ULP, WKILL, WL, WLU, WU, WUL
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Arrays ..
</span> INTEGER IDUMMA( 1 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Functions ..
</span> LOGICAL <a name="LSAME.193"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
INTEGER <a name="ILAENV.194"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
DOUBLE PRECISION <a name="DLAMCH.195"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
EXTERNAL <a name="LSAME.196"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.196"></a><a href="hfy-index.html#ILAENV">ILAENV</a>, <a name="DLAMCH.196"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Subroutines ..
</span> EXTERNAL <a name="DLAEBZ.199"></a><a href="dlaebz.f.html#DLAEBZ.1">DLAEBZ</a>, <a name="XERBLA.199"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Intrinsic Functions ..
</span> INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span> INFO = 0
<span class="comment">*</span><span class="comment">
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?