dlaneg.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 187 行
HTML
187 行
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<title>dlaneg.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.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>
FUNCTION <a name="DLANEG.1"></a><a href="dlaneg.f.html#DLANEG.1">DLANEG</a>( N, D, LLD, SIGMA, PIVMIN, R )
IMPLICIT NONE
INTEGER <a name="DLANEG.3"></a><a href="dlaneg.f.html#DLANEG.1">DLANEG</a>
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> -- LAPACK auxiliary 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">
</span><span class="comment">*</span><span class="comment"> .. Scalar Arguments ..
</span> INTEGER N, R
DOUBLE PRECISION PIVMIN, SIGMA
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Array Arguments ..
</span> DOUBLE PRECISION D( * ), LLD( * )
<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="DLANEG.20"></a><a href="dlaneg.f.html#DLANEG.1">DLANEG</a> computes the Sturm count, the number of negative pivots
</span><span class="comment">*</span><span class="comment"> encountered while factoring tridiagonal T - sigma I = L D L^T.
</span><span class="comment">*</span><span class="comment"> This implementation works directly on the factors without forming
</span><span class="comment">*</span><span class="comment"> the tridiagonal matrix T. The Sturm count is also the number of
</span><span class="comment">*</span><span class="comment"> eigenvalues of T less than sigma.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> This routine is called from <a name="DLARRB.26"></a><a href="dlarrb.f.html#DLARRB.1">DLARRB</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The current routine does not use the PIVMIN parameter but rather
</span><span class="comment">*</span><span class="comment"> requires IEEE-754 propagation of Infinities and NaNs. This
</span><span class="comment">*</span><span class="comment"> routine also has no input range restrictions but does require
</span><span class="comment">*</span><span class="comment"> default exception handling such that x/0 produces Inf when x is
</span><span class="comment">*</span><span class="comment"> non-zero, and Inf/Inf produces NaN. For more information, see:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in
</span><span class="comment">*</span><span class="comment"> Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on
</span><span class="comment">*</span><span class="comment"> Scientific Computing, v28, n5, 2006. DOI 10.1137/050641624
</span><span class="comment">*</span><span class="comment"> (Tech report version in LAWN 172 with the same title.)
</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"> N (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The order of the matrix.
</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 diagonal matrix D.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LLD (input) DOUBLE PRECISION array, dimension (N-1)
</span><span class="comment">*</span><span class="comment"> The (N-1) elements L(i)*L(i)*D(i).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> SIGMA (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment"> Shift amount in T - sigma I = L D L^T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> PIVMIN (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment"> The minimum pivot in the Sturm sequence. May be used
</span><span class="comment">*</span><span class="comment"> when zero pivots are encountered on non-IEEE-754
</span><span class="comment">*</span><span class="comment"> architectures.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> R (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The twist index for the twisted factorization that is used
</span><span class="comment">*</span><span class="comment"> for the negcount.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Further Details
</span><span class="comment">*</span><span class="comment"> ===============
</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"> 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"> Jason Riedy, 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> DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
<span class="comment">*</span><span class="comment"> Some architectures propagate Infinities and NaNs very slowly, so
</span><span class="comment">*</span><span class="comment"> the code computes counts in BLKLEN chunks. Then a NaN can
</span><span class="comment">*</span><span class="comment"> propagate at most BLKLEN columns before being detected. This is
</span><span class="comment">*</span><span class="comment"> not a general tuning parameter; it needs only to be just large
</span><span class="comment">*</span><span class="comment"> enough that the overhead is tiny in common cases.
</span> INTEGER BLKLEN
PARAMETER ( BLKLEN = 128 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Scalars ..
</span> INTEGER BJ, J, NEG1, NEG2, NEGCNT
DOUBLE PRECISION BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP
LOGICAL SAWNAN
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Intrinsic Functions ..
</span> INTRINSIC MIN, MAX
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Functions ..
</span> LOGICAL <a name="DISNAN.93"></a><a href="disnan.f.html#DISNAN.1">DISNAN</a>
EXTERNAL <a name="DISNAN.94"></a><a href="disnan.f.html#DISNAN.1">DISNAN</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Executable Statements ..
</span>
NEGCNT = 0
<span class="comment">*</span><span class="comment"> I) upper part: L D L^T - SIGMA I = L+ D+ L+^T
</span> T = -SIGMA
DO 210 BJ = 1, R-1, BLKLEN
NEG1 = 0
BSAV = T
DO 21 J = BJ, MIN(BJ+BLKLEN-1, R-1)
DPLUS = D( J ) + T
IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
TMP = T / DPLUS
T = TMP * LLD( J ) - SIGMA
21 CONTINUE
SAWNAN = <a name="DISNAN.111"></a><a href="disnan.f.html#DISNAN.1">DISNAN</a>( T )
<span class="comment">*</span><span class="comment"> Run a slower version of the above loop if a NaN is detected.
</span><span class="comment">*</span><span class="comment"> A NaN should occur only with a zero pivot after an infinite
</span><span class="comment">*</span><span class="comment"> pivot. In that case, substituting 1 for T/DPLUS is the
</span><span class="comment">*</span><span class="comment"> correct limit.
</span> IF( SAWNAN ) THEN
NEG1 = 0
T = BSAV
DO 22 J = BJ, MIN(BJ+BLKLEN-1, R-1)
DPLUS = D( J ) + T
IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
TMP = T / DPLUS
IF (<a name="DISNAN.123"></a><a href="disnan.f.html#DISNAN.1">DISNAN</a>(TMP)) TMP = ONE
T = TMP * LLD(J) - SIGMA
22 CONTINUE
END IF
NEGCNT = NEGCNT + NEG1
210 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> II) lower part: L D L^T - SIGMA I = U- D- U-^T
</span> P = D( N ) - SIGMA
DO 230 BJ = N-1, R, -BLKLEN
NEG2 = 0
BSAV = P
DO 23 J = BJ, MAX(BJ-BLKLEN+1, R), -1
DMINUS = LLD( J ) + P
IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
TMP = P / DMINUS
P = TMP * D( J ) - SIGMA
23 CONTINUE
SAWNAN = <a name="DISNAN.141"></a><a href="disnan.f.html#DISNAN.1">DISNAN</a>( P )
<span class="comment">*</span><span class="comment"> As above, run a slower version that substitutes 1 for Inf/Inf.
</span><span class="comment">*</span><span class="comment">
</span> IF( SAWNAN ) THEN
NEG2 = 0
P = BSAV
DO 24 J = BJ, MAX(BJ-BLKLEN+1, R), -1
DMINUS = LLD( J ) + P
IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
TMP = P / DMINUS
IF (<a name="DISNAN.151"></a><a href="disnan.f.html#DISNAN.1">DISNAN</a>(TMP)) TMP = ONE
P = TMP * D(J) - SIGMA
24 CONTINUE
END IF
NEGCNT = NEGCNT + NEG2
230 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> III) Twist index
</span><span class="comment">*</span><span class="comment"> T was shifted by SIGMA initially.
</span> GAMMA = (T + SIGMA) + P
IF( GAMMA.LT.ZERO ) NEGCNT = NEGCNT+1
<a name="DLANEG.163"></a><a href="dlaneg.f.html#DLANEG.1">DLANEG</a> = NEGCNT
END
</pre>
</body>
</html>
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?