LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlaneg.f
Go to the documentation of this file.
00001 *> \brief \b DLANEG
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLANEG + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaneg.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaneg.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaneg.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       INTEGER FUNCTION DLANEG( N, D, LLD, SIGMA, PIVMIN, R )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            N, R
00025 *       DOUBLE PRECISION   PIVMIN, SIGMA
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       DOUBLE PRECISION   D( * ), LLD( * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> DLANEG computes the Sturm count, the number of negative pivots
00038 *> encountered while factoring tridiagonal T - sigma I = L D L^T.
00039 *> This implementation works directly on the factors without forming
00040 *> the tridiagonal matrix T.  The Sturm count is also the number of
00041 *> eigenvalues of T less than sigma.
00042 *>
00043 *> This routine is called from DLARRB.
00044 *>
00045 *> The current routine does not use the PIVMIN parameter but rather
00046 *> requires IEEE-754 propagation of Infinities and NaNs.  This
00047 *> routine also has no input range restrictions but does require
00048 *> default exception handling such that x/0 produces Inf when x is
00049 *> non-zero, and Inf/Inf produces NaN.  For more information, see:
00050 *>
00051 *>   Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in
00052 *>   Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on
00053 *>   Scientific Computing, v28, n5, 2006.  DOI 10.1137/050641624
00054 *>   (Tech report version in LAWN 172 with the same title.)
00055 *> \endverbatim
00056 *
00057 *  Arguments:
00058 *  ==========
00059 *
00060 *> \param[in] N
00061 *> \verbatim
00062 *>          N is INTEGER
00063 *>          The order of the matrix.
00064 *> \endverbatim
00065 *>
00066 *> \param[in] D
00067 *> \verbatim
00068 *>          D is DOUBLE PRECISION array, dimension (N)
00069 *>          The N diagonal elements of the diagonal matrix D.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] LLD
00073 *> \verbatim
00074 *>          LLD is DOUBLE PRECISION array, dimension (N-1)
00075 *>          The (N-1) elements L(i)*L(i)*D(i).
00076 *> \endverbatim
00077 *>
00078 *> \param[in] SIGMA
00079 *> \verbatim
00080 *>          SIGMA is DOUBLE PRECISION
00081 *>          Shift amount in T - sigma I = L D L^T.
00082 *> \endverbatim
00083 *>
00084 *> \param[in] PIVMIN
00085 *> \verbatim
00086 *>          PIVMIN is DOUBLE PRECISION
00087 *>          The minimum pivot in the Sturm sequence.  May be used
00088 *>          when zero pivots are encountered on non-IEEE-754
00089 *>          architectures.
00090 *> \endverbatim
00091 *>
00092 *> \param[in] R
00093 *> \verbatim
00094 *>          R is INTEGER
00095 *>          The twist index for the twisted factorization that is used
00096 *>          for the negcount.
00097 *> \endverbatim
00098 *
00099 *  Authors:
00100 *  ========
00101 *
00102 *> \author Univ. of Tennessee 
00103 *> \author Univ. of California Berkeley 
00104 *> \author Univ. of Colorado Denver 
00105 *> \author NAG Ltd. 
00106 *
00107 *> \date November 2011
00108 *
00109 *> \ingroup auxOTHERauxiliary
00110 *
00111 *> \par Contributors:
00112 *  ==================
00113 *>
00114 *>     Osni Marques, LBNL/NERSC, USA \n
00115 *>     Christof Voemel, University of California, Berkeley, USA \n
00116 *>     Jason Riedy, University of California, Berkeley, USA \n
00117 *>
00118 *  =====================================================================
00119       INTEGER FUNCTION DLANEG( N, D, LLD, SIGMA, PIVMIN, R )
00120 *
00121 *  -- LAPACK auxiliary routine (version 3.4.0) --
00122 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00123 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00124 *     November 2011
00125 *
00126 *     .. Scalar Arguments ..
00127       INTEGER            N, R
00128       DOUBLE PRECISION   PIVMIN, SIGMA
00129 *     ..
00130 *     .. Array Arguments ..
00131       DOUBLE PRECISION   D( * ), LLD( * )
00132 *     ..
00133 *
00134 *  =====================================================================
00135 *
00136 *     .. Parameters ..
00137       DOUBLE PRECISION   ZERO, ONE
00138       PARAMETER        ( ZERO = 0.0D0, ONE = 1.0D0 )
00139 *     Some architectures propagate Infinities and NaNs very slowly, so
00140 *     the code computes counts in BLKLEN chunks.  Then a NaN can
00141 *     propagate at most BLKLEN columns before being detected.  This is
00142 *     not a general tuning parameter; it needs only to be just large
00143 *     enough that the overhead is tiny in common cases.
00144       INTEGER BLKLEN
00145       PARAMETER ( BLKLEN = 128 )
00146 *     ..
00147 *     .. Local Scalars ..
00148       INTEGER            BJ, J, NEG1, NEG2, NEGCNT
00149       DOUBLE PRECISION   BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP
00150       LOGICAL SAWNAN
00151 *     ..
00152 *     .. Intrinsic Functions ..
00153       INTRINSIC MIN, MAX
00154 *     ..
00155 *     .. External Functions ..
00156       LOGICAL DISNAN
00157       EXTERNAL DISNAN
00158 *     ..
00159 *     .. Executable Statements ..
00160 
00161       NEGCNT = 0
00162 
00163 *     I) upper part: L D L^T - SIGMA I = L+ D+ L+^T
00164       T = -SIGMA
00165       DO 210 BJ = 1, R-1, BLKLEN
00166          NEG1 = 0
00167          BSAV = T
00168          DO 21 J = BJ, MIN(BJ+BLKLEN-1, R-1)
00169             DPLUS = D( J ) + T
00170             IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
00171             TMP = T / DPLUS
00172             T = TMP * LLD( J ) - SIGMA
00173  21      CONTINUE
00174          SAWNAN = DISNAN( T )
00175 *     Run a slower version of the above loop if a NaN is detected.
00176 *     A NaN should occur only with a zero pivot after an infinite
00177 *     pivot.  In that case, substituting 1 for T/DPLUS is the
00178 *     correct limit.
00179          IF( SAWNAN ) THEN
00180             NEG1 = 0
00181             T = BSAV
00182             DO 22 J = BJ, MIN(BJ+BLKLEN-1, R-1)
00183                DPLUS = D( J ) + T
00184                IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
00185                TMP = T / DPLUS
00186                IF (DISNAN(TMP)) TMP = ONE
00187                T = TMP * LLD(J) - SIGMA
00188  22         CONTINUE
00189          END IF
00190          NEGCNT = NEGCNT + NEG1
00191  210  CONTINUE
00192 *
00193 *     II) lower part: L D L^T - SIGMA I = U- D- U-^T
00194       P = D( N ) - SIGMA
00195       DO 230 BJ = N-1, R, -BLKLEN
00196          NEG2 = 0
00197          BSAV = P
00198          DO 23 J = BJ, MAX(BJ-BLKLEN+1, R), -1
00199             DMINUS = LLD( J ) + P
00200             IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
00201             TMP = P / DMINUS
00202             P = TMP * D( J ) - SIGMA
00203  23      CONTINUE
00204          SAWNAN = DISNAN( P )
00205 *     As above, run a slower version that substitutes 1 for Inf/Inf.
00206 *
00207          IF( SAWNAN ) THEN
00208             NEG2 = 0
00209             P = BSAV
00210             DO 24 J = BJ, MAX(BJ-BLKLEN+1, R), -1
00211                DMINUS = LLD( J ) + P
00212                IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
00213                TMP = P / DMINUS
00214                IF (DISNAN(TMP)) TMP = ONE
00215                P = TMP * D(J) - SIGMA
00216  24         CONTINUE
00217          END IF
00218          NEGCNT = NEGCNT + NEG2
00219  230  CONTINUE
00220 *
00221 *     III) Twist index
00222 *       T was shifted by SIGMA initially.
00223       GAMMA = (T + SIGMA) + P
00224       IF( GAMMA.LT.ZERO ) NEGCNT = NEGCNT+1
00225 
00226       DLANEG = NEGCNT
00227       END
 All Files Functions