diff options
author | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
---|---|---|
committer | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
commit | baba851215b44ac3b60b9248eb02bcce7eb76247 (patch) | |
tree | 8c0f5c006875532a30d4409f5e94b0f310ff00a7 /SRC/dlaneg.f | |
download | lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.gz lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.bz2 lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.zip |
Move LAPACK trunk into position.
Diffstat (limited to 'SRC/dlaneg.f')
-rw-r--r-- | SRC/dlaneg.f | 164 |
1 files changed, 164 insertions, 0 deletions
diff --git a/SRC/dlaneg.f b/SRC/dlaneg.f new file mode 100644 index 00000000..fead657c --- /dev/null +++ b/SRC/dlaneg.f @@ -0,0 +1,164 @@ + FUNCTION DLANEG( N, D, LLD, SIGMA, PIVMIN, R ) + IMPLICIT NONE + INTEGER DLANEG +* +* -- LAPACK auxiliary routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER N, R + DOUBLE PRECISION PIVMIN, SIGMA +* .. +* .. Array Arguments .. + DOUBLE PRECISION D( * ), LLD( * ) +* .. +* +* Purpose +* ======= +* +* DLANEG computes the Sturm count, the number of negative pivots +* encountered while factoring tridiagonal T - sigma I = L D L^T. +* This implementation works directly on the factors without forming +* the tridiagonal matrix T. The Sturm count is also the number of +* eigenvalues of T less than sigma. +* +* This routine is called from DLARRB. +* +* The current routine does not use the PIVMIN parameter but rather +* requires IEEE-754 propagation of Infinities and NaNs. This +* routine also has no input range restrictions but does require +* default exception handling such that x/0 produces Inf when x is +* non-zero, and Inf/Inf produces NaN. For more information, see: +* +* Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in +* Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on +* Scientific Computing, v28, n5, 2006. DOI 10.1137/050641624 +* (Tech report version in LAWN 172 with the same title.) +* +* Arguments +* ========= +* +* N (input) INTEGER +* The order of the matrix. +* +* D (input) DOUBLE PRECISION array, dimension (N) +* The N diagonal elements of the diagonal matrix D. +* +* LLD (input) DOUBLE PRECISION array, dimension (N-1) +* The (N-1) elements L(i)*L(i)*D(i). +* +* SIGMA (input) DOUBLE PRECISION +* Shift amount in T - sigma I = L D L^T. +* +* PIVMIN (input) DOUBLE PRECISION +* The minimum pivot in the Sturm sequence. May be used +* when zero pivots are encountered on non-IEEE-754 +* architectures. +* +* R (input) INTEGER +* The twist index for the twisted factorization that is used +* for the negcount. +* +* Further Details +* =============== +* +* Based on contributions by +* Osni Marques, LBNL/NERSC, USA +* Christof Voemel, University of California, Berkeley, USA +* Jason Riedy, University of California, Berkeley, USA +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +* Some architectures propagate Infinities and NaNs very slowly, so +* the code computes counts in BLKLEN chunks. Then a NaN can +* propagate at most BLKLEN columns before being detected. This is +* not a general tuning parameter; it needs only to be just large +* enough that the overhead is tiny in common cases. + INTEGER BLKLEN + PARAMETER ( BLKLEN = 128 ) +* .. +* .. Local Scalars .. + INTEGER BJ, J, NEG1, NEG2, NEGCNT + DOUBLE PRECISION BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP + LOGICAL SAWNAN +* .. +* .. Intrinsic Functions .. + INTRINSIC MIN, MAX +* .. +* .. External Functions .. + LOGICAL DISNAN + EXTERNAL DISNAN +* .. +* .. Executable Statements .. + + NEGCNT = 0 + +* I) upper part: L D L^T - SIGMA I = L+ D+ L+^T + 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 = DISNAN( T ) +* Run a slower version of the above loop if a NaN is detected. +* A NaN should occur only with a zero pivot after an infinite +* pivot. In that case, substituting 1 for T/DPLUS is the +* correct limit. + 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 (DISNAN(TMP)) TMP = ONE + T = TMP * LLD(J) - SIGMA + 22 CONTINUE + END IF + NEGCNT = NEGCNT + NEG1 + 210 CONTINUE +* +* II) lower part: L D L^T - SIGMA I = U- D- U-^T + 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 = DISNAN( P ) +* As above, run a slower version that substitutes 1 for Inf/Inf. +* + 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 (DISNAN(TMP)) TMP = ONE + P = TMP * D(J) - SIGMA + 24 CONTINUE + END IF + NEGCNT = NEGCNT + NEG2 + 230 CONTINUE +* +* III) Twist index +* T was shifted by SIGMA initially. + GAMMA = (T + SIGMA) + P + IF( GAMMA.LT.ZERO ) NEGCNT = NEGCNT+1 + + DLANEG = NEGCNT + END |