LAPACK 3.3.0

dlaneg.f

Go to the documentation of this file.
00001       INTEGER FUNCTION DLANEG( N, D, LLD, SIGMA, PIVMIN, R )
00002       IMPLICIT NONE
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.2.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     June 2010
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            N, R
00011       DOUBLE PRECISION   PIVMIN, SIGMA
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   D( * ), LLD( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DLANEG computes the Sturm count, the number of negative pivots
00021 *  encountered while factoring tridiagonal T - sigma I = L D L^T.
00022 *  This implementation works directly on the factors without forming
00023 *  the tridiagonal matrix T.  The Sturm count is also the number of
00024 *  eigenvalues of T less than sigma.
00025 *
00026 *  This routine is called from DLARRB.
00027 *
00028 *  The current routine does not use the PIVMIN parameter but rather
00029 *  requires IEEE-754 propagation of Infinities and NaNs.  This
00030 *  routine also has no input range restrictions but does require
00031 *  default exception handling such that x/0 produces Inf when x is
00032 *  non-zero, and Inf/Inf produces NaN.  For more information, see:
00033 *
00034 *    Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in
00035 *    Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on
00036 *    Scientific Computing, v28, n5, 2006.  DOI 10.1137/050641624
00037 *    (Tech report version in LAWN 172 with the same title.)
00038 *
00039 *  Arguments
00040 *  =========
00041 *
00042 *  N       (input) INTEGER
00043 *          The order of the matrix.
00044 *
00045 *  D       (input) DOUBLE PRECISION array, dimension (N)
00046 *          The N diagonal elements of the diagonal matrix D.
00047 *
00048 *  LLD     (input) DOUBLE PRECISION array, dimension (N-1)
00049 *          The (N-1) elements L(i)*L(i)*D(i).
00050 *
00051 *  SIGMA   (input) DOUBLE PRECISION
00052 *          Shift amount in T - sigma I = L D L^T.
00053 *
00054 *  PIVMIN  (input) DOUBLE PRECISION
00055 *          The minimum pivot in the Sturm sequence.  May be used
00056 *          when zero pivots are encountered on non-IEEE-754
00057 *          architectures.
00058 *
00059 *  R       (input) INTEGER
00060 *          The twist index for the twisted factorization that is used
00061 *          for the negcount.
00062 *
00063 *  Further Details
00064 *  ===============
00065 *
00066 *  Based on contributions by
00067 *     Osni Marques, LBNL/NERSC, USA
00068 *     Christof Voemel, University of California, Berkeley, USA
00069 *     Jason Riedy, University of California, Berkeley, USA
00070 *
00071 *  =====================================================================
00072 *
00073 *     .. Parameters ..
00074       DOUBLE PRECISION   ZERO, ONE
00075       PARAMETER        ( ZERO = 0.0D0, ONE = 1.0D0 )
00076 *     Some architectures propagate Infinities and NaNs very slowly, so
00077 *     the code computes counts in BLKLEN chunks.  Then a NaN can
00078 *     propagate at most BLKLEN columns before being detected.  This is
00079 *     not a general tuning parameter; it needs only to be just large
00080 *     enough that the overhead is tiny in common cases.
00081       INTEGER BLKLEN
00082       PARAMETER ( BLKLEN = 128 )
00083 *     ..
00084 *     .. Local Scalars ..
00085       INTEGER            BJ, J, NEG1, NEG2, NEGCNT
00086       DOUBLE PRECISION   BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP
00087       LOGICAL SAWNAN
00088 *     ..
00089 *     .. Intrinsic Functions ..
00090       INTRINSIC MIN, MAX
00091 *     ..
00092 *     .. External Functions ..
00093       LOGICAL DISNAN
00094       EXTERNAL DISNAN
00095 *     ..
00096 *     .. Executable Statements ..
00097 
00098       NEGCNT = 0
00099 
00100 *     I) upper part: L D L^T - SIGMA I = L+ D+ L+^T
00101       T = -SIGMA
00102       DO 210 BJ = 1, R-1, BLKLEN
00103          NEG1 = 0
00104          BSAV = T
00105          DO 21 J = BJ, MIN(BJ+BLKLEN-1, R-1)
00106             DPLUS = D( J ) + T
00107             IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
00108             TMP = T / DPLUS
00109             T = TMP * LLD( J ) - SIGMA
00110  21      CONTINUE
00111          SAWNAN = DISNAN( T )
00112 *     Run a slower version of the above loop if a NaN is detected.
00113 *     A NaN should occur only with a zero pivot after an infinite
00114 *     pivot.  In that case, substituting 1 for T/DPLUS is the
00115 *     correct limit.
00116          IF( SAWNAN ) THEN
00117             NEG1 = 0
00118             T = BSAV
00119             DO 22 J = BJ, MIN(BJ+BLKLEN-1, R-1)
00120                DPLUS = D( J ) + T
00121                IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1
00122                TMP = T / DPLUS
00123                IF (DISNAN(TMP)) TMP = ONE
00124                T = TMP * LLD(J) - SIGMA
00125  22         CONTINUE
00126          END IF
00127          NEGCNT = NEGCNT + NEG1
00128  210  CONTINUE
00129 *
00130 *     II) lower part: L D L^T - SIGMA I = U- D- U-^T
00131       P = D( N ) - SIGMA
00132       DO 230 BJ = N-1, R, -BLKLEN
00133          NEG2 = 0
00134          BSAV = P
00135          DO 23 J = BJ, MAX(BJ-BLKLEN+1, R), -1
00136             DMINUS = LLD( J ) + P
00137             IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
00138             TMP = P / DMINUS
00139             P = TMP * D( J ) - SIGMA
00140  23      CONTINUE
00141          SAWNAN = DISNAN( P )
00142 *     As above, run a slower version that substitutes 1 for Inf/Inf.
00143 *
00144          IF( SAWNAN ) THEN
00145             NEG2 = 0
00146             P = BSAV
00147             DO 24 J = BJ, MAX(BJ-BLKLEN+1, R), -1
00148                DMINUS = LLD( J ) + P
00149                IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1
00150                TMP = P / DMINUS
00151                IF (DISNAN(TMP)) TMP = ONE
00152                P = TMP * D(J) - SIGMA
00153  24         CONTINUE
00154          END IF
00155          NEGCNT = NEGCNT + NEG2
00156  230  CONTINUE
00157 *
00158 *     III) Twist index
00159 *       T was shifted by SIGMA initially.
00160       GAMMA = (T + SIGMA) + P
00161       IF( GAMMA.LT.ZERO ) NEGCNT = NEGCNT+1
00162 
00163       DLANEG = NEGCNT
00164       END
 All Files Functions