![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.2.2.
1: INTEGER FUNCTION DLANEG( N, D, LLD, SIGMA, PIVMIN, R ) 2: IMPLICIT NONE 3: * 4: * -- LAPACK auxiliary routine (version 3.2.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * June 2010 8: * 9: * .. Scalar Arguments .. 10: INTEGER N, R 11: DOUBLE PRECISION PIVMIN, SIGMA 12: * .. 13: * .. Array Arguments .. 14: DOUBLE PRECISION D( * ), LLD( * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * DLANEG computes the Sturm count, the number of negative pivots 21: * encountered while factoring tridiagonal T - sigma I = L D L^T. 22: * This implementation works directly on the factors without forming 23: * the tridiagonal matrix T. The Sturm count is also the number of 24: * eigenvalues of T less than sigma. 25: * 26: * This routine is called from DLARRB. 27: * 28: * The current routine does not use the PIVMIN parameter but rather 29: * requires IEEE-754 propagation of Infinities and NaNs. This 30: * routine also has no input range restrictions but does require 31: * default exception handling such that x/0 produces Inf when x is 32: * non-zero, and Inf/Inf produces NaN. For more information, see: 33: * 34: * Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in 35: * Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on 36: * Scientific Computing, v28, n5, 2006. DOI 10.1137/050641624 37: * (Tech report version in LAWN 172 with the same title.) 38: * 39: * Arguments 40: * ========= 41: * 42: * N (input) INTEGER 43: * The order of the matrix. 44: * 45: * D (input) DOUBLE PRECISION array, dimension (N) 46: * The N diagonal elements of the diagonal matrix D. 47: * 48: * LLD (input) DOUBLE PRECISION array, dimension (N-1) 49: * The (N-1) elements L(i)*L(i)*D(i). 50: * 51: * SIGMA (input) DOUBLE PRECISION 52: * Shift amount in T - sigma I = L D L^T. 53: * 54: * PIVMIN (input) DOUBLE PRECISION 55: * The minimum pivot in the Sturm sequence. May be used 56: * when zero pivots are encountered on non-IEEE-754 57: * architectures. 58: * 59: * R (input) INTEGER 60: * The twist index for the twisted factorization that is used 61: * for the negcount. 62: * 63: * Further Details 64: * =============== 65: * 66: * Based on contributions by 67: * Osni Marques, LBNL/NERSC, USA 68: * Christof Voemel, University of California, Berkeley, USA 69: * Jason Riedy, University of California, Berkeley, USA 70: * 71: * ===================================================================== 72: * 73: * .. Parameters .. 74: DOUBLE PRECISION ZERO, ONE 75: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 76: * Some architectures propagate Infinities and NaNs very slowly, so 77: * the code computes counts in BLKLEN chunks. Then a NaN can 78: * propagate at most BLKLEN columns before being detected. This is 79: * not a general tuning parameter; it needs only to be just large 80: * enough that the overhead is tiny in common cases. 81: INTEGER BLKLEN 82: PARAMETER ( BLKLEN = 128 ) 83: * .. 84: * .. Local Scalars .. 85: INTEGER BJ, J, NEG1, NEG2, NEGCNT 86: DOUBLE PRECISION BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP 87: LOGICAL SAWNAN 88: * .. 89: * .. Intrinsic Functions .. 90: INTRINSIC MIN, MAX 91: * .. 92: * .. External Functions .. 93: LOGICAL DISNAN 94: EXTERNAL DISNAN 95: * .. 96: * .. Executable Statements .. 97: 98: NEGCNT = 0 99: 100: * I) upper part: L D L^T - SIGMA I = L+ D+ L+^T 101: T = -SIGMA 102: DO 210 BJ = 1, R-1, BLKLEN 103: NEG1 = 0 104: BSAV = T 105: DO 21 J = BJ, MIN(BJ+BLKLEN-1, R-1) 106: DPLUS = D( J ) + T 107: IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1 108: TMP = T / DPLUS 109: T = TMP * LLD( J ) - SIGMA 110: 21 CONTINUE 111: SAWNAN = DISNAN( T ) 112: * Run a slower version of the above loop if a NaN is detected. 113: * A NaN should occur only with a zero pivot after an infinite 114: * pivot. In that case, substituting 1 for T/DPLUS is the 115: * correct limit. 116: IF( SAWNAN ) THEN 117: NEG1 = 0 118: T = BSAV 119: DO 22 J = BJ, MIN(BJ+BLKLEN-1, R-1) 120: DPLUS = D( J ) + T 121: IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1 122: TMP = T / DPLUS 123: IF (DISNAN(TMP)) TMP = ONE 124: T = TMP * LLD(J) - SIGMA 125: 22 CONTINUE 126: END IF 127: NEGCNT = NEGCNT + NEG1 128: 210 CONTINUE 129: * 130: * II) lower part: L D L^T - SIGMA I = U- D- U-^T 131: P = D( N ) - SIGMA 132: DO 230 BJ = N-1, R, -BLKLEN 133: NEG2 = 0 134: BSAV = P 135: DO 23 J = BJ, MAX(BJ-BLKLEN+1, R), -1 136: DMINUS = LLD( J ) + P 137: IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1 138: TMP = P / DMINUS 139: P = TMP * D( J ) - SIGMA 140: 23 CONTINUE 141: SAWNAN = DISNAN( P ) 142: * As above, run a slower version that substitutes 1 for Inf/Inf. 143: * 144: IF( SAWNAN ) THEN 145: NEG2 = 0 146: P = BSAV 147: DO 24 J = BJ, MAX(BJ-BLKLEN+1, R), -1 148: DMINUS = LLD( J ) + P 149: IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1 150: TMP = P / DMINUS 151: IF (DISNAN(TMP)) TMP = ONE 152: P = TMP * D(J) - SIGMA 153: 24 CONTINUE 154: END IF 155: NEGCNT = NEGCNT + NEG2 156: 230 CONTINUE 157: * 158: * III) Twist index 159: * T was shifted by SIGMA initially. 160: GAMMA = (T + SIGMA) + P 161: IF( GAMMA.LT.ZERO ) NEGCNT = NEGCNT+1 162: 163: DLANEG = NEGCNT 164: END