slarrd.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 738 行 · 第 1/4 页
HTML
738 行
</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"> FUDGE REAL , 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"> Based on contributions by
</span><span class="comment">*</span><span class="comment"> W. Kahan, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment"> Beresford Parlett, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment"> Jim Demmel, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment"> Inderjit Dhillon, University of Texas, Austin, USA
</span><span class="comment">*</span><span class="comment"> Osni Marques, LBNL/NERSC, USA
</span><span class="comment">*</span><span class="comment"> Christof Voemel, University of California, Berkeley, USA
</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> REAL ZERO, ONE, TWO, HALF, FUDGE
PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0,
$ TWO = 2.0E0, HALF = ONE/TWO,
$ FUDGE = TWO )
INTEGER ALLRNG, VALRNG, INDRNG
PARAMETER ( ALLRNG = 1, VALRNG = 2, INDRNG = 3 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Scalars ..
</span> LOGICAL NCNVRG, TOOFEW
INTEGER I, IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
$ IM, IN, IOFF, IOUT, IRANGE, ITMAX, ITMP1,
$ ITMP2, IW, IWOFF, J, JBLK, JDISC, JE, JEE, NB,
$ NWL, NWU
REAL ATOLI, EPS, GL, GU, RTOLI, SPDIAM, TMP1, TMP2,
$ TNORM, UFLOW, WKILL, WLU, 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.229"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
INTEGER <a name="ILAENV.230"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
REAL <a name="SLAMCH.231"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
EXTERNAL <a name="LSAME.232"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.232"></a><a href="hfy-index.html#ILAENV">ILAENV</a>, <a name="SLAMCH.232"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Subroutines ..
</span> EXTERNAL <a name="SLAEBZ.235"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</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
<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">
</span><span class="comment">*</span><span class="comment"> Decode RANGE
</span><span class="comment">*</span><span class="comment">
</span> IF( <a name="LSAME.246"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'A'</span> ) ) THEN
IRANGE = ALLRNG
ELSE IF( <a name="LSAME.248"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'V'</span> ) ) THEN
IRANGE = VALRNG
ELSE IF( <a name="LSAME.250"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( RANGE, <span class="string">'I'</span> ) ) THEN
IRANGE = INDRNG
ELSE
IRANGE = 0
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Check for Errors
</span><span class="comment">*</span><span class="comment">
</span> IF( IRANGE.LE.0 ) THEN
INFO = -1
ELSE IF( .NOT.(<a name="LSAME.260"></a><a href="lsame.f.html#LSAME.1">LSAME</a>(ORDER,<span class="string">'B'</span>).OR.<a name="LSAME.260"></a><a href="lsame.f.html#LSAME.1">LSAME</a>(ORDER,<span class="string">'E'</span>)) ) THEN
INFO = -2
ELSE IF( N.LT.0 ) THEN
INFO = -3
ELSE IF( IRANGE.EQ.VALRNG ) THEN
IF( VL.GE.VU )
$ INFO = -5
ELSE IF( IRANGE.EQ.INDRNG .AND.
$ ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) THEN
INFO = -6
ELSE IF( IRANGE.EQ.INDRNG .AND.
$ ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN
INFO = -7
END IF
<span class="comment">*</span><span class="comment">
</span> IF( INFO.NE.0 ) THEN
RETURN
END IF
<span class="comment">*</span><span class="comment"> Initialize error flags
</span> INFO = 0
NCNVRG = .FALSE.
TOOFEW = .FALSE.
<span class="comment">*</span><span class="comment"> Quick return if possible
</span> M = 0
IF( N.EQ.0 ) RETURN
<span class="comment">*</span><span class="comment"> Simplification:
</span> IF( IRANGE.EQ.INDRNG .AND. IL.EQ.1 .AND. IU.EQ.N ) IRANGE = 1
<span class="comment">*</span><span class="comment"> Get machine constants
</span> EPS = <a name="SLAMCH.292"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )
UFLOW = <a name="SLAMCH.293"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'U'</span> )
<span class="comment">*</span><span class="comment"> Special Case when N=1
</span><span class="comment">*</span><span class="comment"> Treat case of 1x1 matrix for quick return
</span> IF( N.EQ.1 ) THEN
IF( (IRANGE.EQ.ALLRNG).OR.
$ ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
$ ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
M = 1
W(1) = D(1)
<span class="comment">*</span><span class="comment"> The computation error of the eigenvalue is zero
</span> WERR(1) = ZERO
IBLOCK( 1 ) = 1
INDEXW( 1 ) = 1
ENDIF
RETURN
END IF
<span class="comment">*</span><span class="comment"> NB is the minimum vector length for vector bisection, or 0
</span><span class="comment">*</span><span class="comment"> if only scalar is to be done.
</span> NB = <a name="ILAENV.314"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="SSTEBZ.314"></a><a href="sstebz.f.html#SSTEBZ.1">SSTEBZ</a>'</span>, <span class="string">' '</span>, N, -1, -1, -1 )
IF( NB.LE.1 ) NB = 0
<span class="comment">*</span><span class="comment"> Find global spectral radius
</span> GL = D(1)
GU = D(1)
DO 5 I = 1,N
GL = MIN( GL, GERS( 2*I - 1))
GU = MAX( GU, GERS(2*I) )
5 CONTINUE
<span class="comment">*</span><span class="comment"> Compute global Gerschgorin bounds and spectral diameter
</span> TNORM = MAX( ABS( GL ), ABS( GU ) )
GL = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
GU = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
SPDIAM = GU - GL
<span class="comment">*</span><span class="comment"> Input arguments for <a name="SLAEBZ.329"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a>:
</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 < RELTOL*max(|a|,|b|),
</span> RTOLI = RELTOL
<span class="comment">*</span><span class="comment"> Set the absolute tolerance for interval convergence to zero to force
</span><span class="comment">*</span><span class="comment"> interval convergence based on relative size of the interval.
</span><span class="comment">*</span><span class="comment"> This is dangerous because intervals might not converge when RELTOL is
</span><span class="comment">*</span><span class="comment"> small. But at least a very small number should be selected so that for
</span><span class="comment">*</span><span class="comment"> strongly graded matrices, the code can get relatively accurate
</span><span class="comment">*</span><span class="comment"> eigenvalues.
</span> ATOLI = FUDGE*TWO*UFLOW + FUDGE*TWO*PIVMIN
IF( IRANGE.EQ.INDRNG ) THEN
<span class="comment">*</span><span class="comment"> RANGE='I': Compute an interval containing eigenvalues
</span><span class="comment">*</span><span class="comment"> IL through IU. The initial interval [GL,GU] from the global
</span><span class="comment">*</span><span class="comment"> Gerschgorin bounds GL and GU is refined by <a name="SLAEBZ.345"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a>.
</span> ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
$ LOG( TWO ) ) + 2
WORK( N+1 ) = GL
WORK( N+2 ) = GL
WORK( N+3 ) = GU
WORK( N+4 ) = GU
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?