Annotation of rpl/lapack/lapack/dlaneg.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DLANEG
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLANEG + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaneg.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaneg.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaneg.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * INTEGER FUNCTION DLANEG( N, D, LLD, SIGMA, PIVMIN, R )
! 22: *
! 23: * .. Scalar Arguments ..
! 24: * INTEGER N, R
! 25: * DOUBLE PRECISION PIVMIN, SIGMA
! 26: * ..
! 27: * .. Array Arguments ..
! 28: * DOUBLE PRECISION D( * ), LLD( * )
! 29: * ..
! 30: *
! 31: *
! 32: *> \par Purpose:
! 33: * =============
! 34: *>
! 35: *> \verbatim
! 36: *>
! 37: *> DLANEG computes the Sturm count, the number of negative pivots
! 38: *> encountered while factoring tridiagonal T - sigma I = L D L^T.
! 39: *> This implementation works directly on the factors without forming
! 40: *> the tridiagonal matrix T. The Sturm count is also the number of
! 41: *> eigenvalues of T less than sigma.
! 42: *>
! 43: *> This routine is called from DLARRB.
! 44: *>
! 45: *> The current routine does not use the PIVMIN parameter but rather
! 46: *> requires IEEE-754 propagation of Infinities and NaNs. This
! 47: *> routine also has no input range restrictions but does require
! 48: *> default exception handling such that x/0 produces Inf when x is
! 49: *> non-zero, and Inf/Inf produces NaN. For more information, see:
! 50: *>
! 51: *> Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in
! 52: *> Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on
! 53: *> Scientific Computing, v28, n5, 2006. DOI 10.1137/050641624
! 54: *> (Tech report version in LAWN 172 with the same title.)
! 55: *> \endverbatim
! 56: *
! 57: * Arguments:
! 58: * ==========
! 59: *
! 60: *> \param[in] N
! 61: *> \verbatim
! 62: *> N is INTEGER
! 63: *> The order of the matrix.
! 64: *> \endverbatim
! 65: *>
! 66: *> \param[in] D
! 67: *> \verbatim
! 68: *> D is DOUBLE PRECISION array, dimension (N)
! 69: *> The N diagonal elements of the diagonal matrix D.
! 70: *> \endverbatim
! 71: *>
! 72: *> \param[in] LLD
! 73: *> \verbatim
! 74: *> LLD is DOUBLE PRECISION array, dimension (N-1)
! 75: *> The (N-1) elements L(i)*L(i)*D(i).
! 76: *> \endverbatim
! 77: *>
! 78: *> \param[in] SIGMA
! 79: *> \verbatim
! 80: *> SIGMA is DOUBLE PRECISION
! 81: *> Shift amount in T - sigma I = L D L^T.
! 82: *> \endverbatim
! 83: *>
! 84: *> \param[in] PIVMIN
! 85: *> \verbatim
! 86: *> PIVMIN is DOUBLE PRECISION
! 87: *> The minimum pivot in the Sturm sequence. May be used
! 88: *> when zero pivots are encountered on non-IEEE-754
! 89: *> architectures.
! 90: *> \endverbatim
! 91: *>
! 92: *> \param[in] R
! 93: *> \verbatim
! 94: *> R is INTEGER
! 95: *> The twist index for the twisted factorization that is used
! 96: *> for the negcount.
! 97: *> \endverbatim
! 98: *
! 99: * Authors:
! 100: * ========
! 101: *
! 102: *> \author Univ. of Tennessee
! 103: *> \author Univ. of California Berkeley
! 104: *> \author Univ. of Colorado Denver
! 105: *> \author NAG Ltd.
! 106: *
! 107: *> \date November 2011
! 108: *
! 109: *> \ingroup auxOTHERauxiliary
! 110: *
! 111: *> \par Contributors:
! 112: * ==================
! 113: *>
! 114: *> Osni Marques, LBNL/NERSC, USA \n
! 115: *> Christof Voemel, University of California, Berkeley, USA \n
! 116: *> Jason Riedy, University of California, Berkeley, USA \n
! 117: *>
! 118: * =====================================================================
1.5 bertrand 119: INTEGER FUNCTION DLANEG( N, D, LLD, SIGMA, PIVMIN, R )
1.1 bertrand 120: *
1.9 ! bertrand 121: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 122: * -- LAPACK is a software package provided by Univ. of Tennessee, --
123: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 124: * November 2011
1.1 bertrand 125: *
126: * .. Scalar Arguments ..
127: INTEGER N, R
128: DOUBLE PRECISION PIVMIN, SIGMA
129: * ..
130: * .. Array Arguments ..
131: DOUBLE PRECISION D( * ), LLD( * )
132: * ..
133: *
134: * =====================================================================
135: *
136: * .. Parameters ..
137: DOUBLE PRECISION ZERO, ONE
138: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
139: * Some architectures propagate Infinities and NaNs very slowly, so
140: * the code computes counts in BLKLEN chunks. Then a NaN can
141: * propagate at most BLKLEN columns before being detected. This is
142: * not a general tuning parameter; it needs only to be just large
143: * enough that the overhead is tiny in common cases.
144: INTEGER BLKLEN
145: PARAMETER ( BLKLEN = 128 )
146: * ..
147: * .. Local Scalars ..
148: INTEGER BJ, J, NEG1, NEG2, NEGCNT
149: DOUBLE PRECISION BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP
150: LOGICAL SAWNAN
151: * ..
152: * .. Intrinsic Functions ..
153: INTRINSIC MIN, MAX
154: * ..
155: * .. External Functions ..
156: LOGICAL DISNAN
157: EXTERNAL DISNAN
158: * ..
159: * .. Executable Statements ..
160:
161: NEGCNT = 0
162:
163: * I) upper part: L D L^T - SIGMA I = L+ D+ L+^T
164: T = -SIGMA
165: DO 210 BJ = 1, R-1, BLKLEN
166: NEG1 = 0
167: BSAV = T
168: DO 21 J = BJ, MIN(BJ+BLKLEN-1, R-1)
169: DPLUS = D( J ) + T
170: IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
171: TMP = T / DPLUS
172: T = TMP * LLD( J ) - SIGMA
173: 21 CONTINUE
174: SAWNAN = DISNAN( T )
175: * Run a slower version of the above loop if a NaN is detected.
176: * A NaN should occur only with a zero pivot after an infinite
177: * pivot. In that case, substituting 1 for T/DPLUS is the
178: * correct limit.
179: IF( SAWNAN ) THEN
180: NEG1 = 0
181: T = BSAV
182: DO 22 J = BJ, MIN(BJ+BLKLEN-1, R-1)
183: DPLUS = D( J ) + T
184: IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
185: TMP = T / DPLUS
186: IF (DISNAN(TMP)) TMP = ONE
187: T = TMP * LLD(J) - SIGMA
188: 22 CONTINUE
189: END IF
190: NEGCNT = NEGCNT + NEG1
191: 210 CONTINUE
192: *
193: * II) lower part: L D L^T - SIGMA I = U- D- U-^T
194: P = D( N ) - SIGMA
195: DO 230 BJ = N-1, R, -BLKLEN
196: NEG2 = 0
197: BSAV = P
198: DO 23 J = BJ, MAX(BJ-BLKLEN+1, R), -1
199: DMINUS = LLD( J ) + P
200: IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
201: TMP = P / DMINUS
202: P = TMP * D( J ) - SIGMA
203: 23 CONTINUE
204: SAWNAN = DISNAN( P )
205: * As above, run a slower version that substitutes 1 for Inf/Inf.
206: *
207: IF( SAWNAN ) THEN
208: NEG2 = 0
209: P = BSAV
210: DO 24 J = BJ, MAX(BJ-BLKLEN+1, R), -1
211: DMINUS = LLD( J ) + P
212: IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
213: TMP = P / DMINUS
214: IF (DISNAN(TMP)) TMP = ONE
215: P = TMP * D(J) - SIGMA
216: 24 CONTINUE
217: END IF
218: NEGCNT = NEGCNT + NEG2
219: 230 CONTINUE
220: *
221: * III) Twist index
222: * T was shifted by SIGMA initially.
223: GAMMA = (T + SIGMA) + P
224: IF( GAMMA.LT.ZERO ) NEGCNT = NEGCNT+1
225:
226: DLANEG = NEGCNT
227: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaneg.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack