SRC/slaneg.f

   1 *> \brief \b SLANEG computes the Sturm count.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SLANEG + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaneg.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaneg.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaneg.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       INTEGER FUNCTION SLANEG( N, D, LLD, SIGMA, PIVMIN, R )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            N, R
  25 *       REAL               PIVMIN, SIGMA
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       REAL               D( * ), LLD( * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> SLANEG 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 SLARRB.
  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 REAL 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 REAL 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 REAL
  81 *>          Shift amount in T - sigma I = L D L^T.
  82 *> \endverbatim
  83 *>
  84 *> \param[in] PIVMIN
  85 *> \verbatim
  86 *>          PIVMIN is REAL
  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 September 2012
 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 *  =====================================================================
 119       INTEGER FUNCTION SLANEG( N, D, LLD, SIGMA, PIVMIN, R )
 120 *
 121 *  -- LAPACK auxiliary routine (version 3.4.2) --
 122 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 123 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 124 *     September 2012
 125 *
 126 *     .. Scalar Arguments ..
 127       INTEGER            N, R
 128       REAL               PIVMIN, SIGMA
 129 *     ..
 130 *     .. Array Arguments ..
 131       REAL               D( * ), LLD( * )
 132 *     ..
 133 *
 134 *  =====================================================================
 135 *
 136 *     .. Parameters ..
 137       REAL               ZERO, ONE
 138       PARAMETER        ( ZERO = 0.0E0, ONE = 1.0E0 )
 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       REAL               BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP
 150       LOGICAL SAWNAN
 151 *     ..
 152 *     .. Intrinsic Functions ..
 153       INTRINSIC MIN, MAX
 154 *     ..
 155 *     .. External Functions ..
 156       LOGICAL SISNAN
 157       EXTERNAL SISNAN
 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 = SISNAN( 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 (SISNAN(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 = SISNAN( 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 (SISNAN(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       SLANEG = NEGCNT
 227       END