*> \brief \b DLANEG * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download DLANEG + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * INTEGER FUNCTION DLANEG( N, D, LLD, SIGMA, PIVMIN, R ) * * .. Scalar Arguments .. * INTEGER N, R * DOUBLE PRECISION PIVMIN, SIGMA * .. * .. Array Arguments .. * DOUBLE PRECISION D( * ), LLD( * ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> 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.) *> \endverbatim * * Arguments: * ========== * *> \param[in] N *> \verbatim *> N is INTEGER *> The order of the matrix. *> \endverbatim *> *> \param[in] D *> \verbatim *> D is DOUBLE PRECISION array, dimension (N) *> The N diagonal elements of the diagonal matrix D. *> \endverbatim *> *> \param[in] LLD *> \verbatim *> LLD is DOUBLE PRECISION array, dimension (N-1) *> The (N-1) elements L(i)*L(i)*D(i). *> \endverbatim *> *> \param[in] SIGMA *> \verbatim *> SIGMA is DOUBLE PRECISION *> Shift amount in T - sigma I = L D L^T. *> \endverbatim *> *> \param[in] PIVMIN *> \verbatim *> PIVMIN is DOUBLE PRECISION *> The minimum pivot in the Sturm sequence. May be used *> when zero pivots are encountered on non-IEEE-754 *> architectures. *> \endverbatim *> *> \param[in] R *> \verbatim *> R is INTEGER *> The twist index for the twisted factorization that is used *> for the negcount. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2011 * *> \ingroup auxOTHERauxiliary * *> \par Contributors: * ================== *> *> Osni Marques, LBNL/NERSC, USA \n *> Christof Voemel, University of California, Berkeley, USA \n *> Jason Riedy, University of California, Berkeley, USA \n *> * ===================================================================== INTEGER FUNCTION DLANEG( N, D, LLD, SIGMA, PIVMIN, R ) * * -- LAPACK auxiliary routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2011 * * .. Scalar Arguments .. INTEGER N, R DOUBLE PRECISION PIVMIN, SIGMA * .. * .. Array Arguments .. DOUBLE PRECISION D( * ), LLD( * ) * .. * * ===================================================================== * * .. 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