|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 |
INTEGER FUNCTION SLANEG( N, D, LLD, SIGMA, PIVMIN, R )
IMPLICIT NONE * * -- LAPACK auxiliary routine (version 3.2.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * June 2010 * * .. Scalar Arguments .. INTEGER N, R REAL PIVMIN, SIGMA * .. * .. Array Arguments .. REAL D( * ), LLD( * ) * .. * * Purpose * ======= * * SLANEG 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 SLARRB. * * 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) REAL array, dimension (N) * The N diagonal elements of the diagonal matrix D. * * LLD (input) REAL array, dimension (N-1) * The (N-1) elements L(i)*L(i)*D(i). * * SIGMA (input) REAL * Shift amount in T - sigma I = L D L^T. * * PIVMIN (input) REAL * 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 .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) * 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 REAL BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP LOGICAL SAWNAN * .. * .. Intrinsic Functions .. INTRINSIC MIN, MAX * .. * .. External Functions .. LOGICAL SISNAN EXTERNAL SISNAN * .. * .. 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 = SISNAN( 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 (SISNAN(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 = SISNAN( 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 (SISNAN(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 SLANEG = NEGCNT END |