slaneg()

integer function slaneg	(	integer	n,
		real, dimension( * )	d,
		real, dimension( * )	lld,
		real	sigma,
		real	pivmin,
		integer	r )

SLANEG computes the Sturm count.

Download SLANEG + dependencies [TGZ] [ZIP] [TXT]

Purpose:

!>
!> SLANEG 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 SLARRB.
!>
!> 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,  SIAM Journal on
!>   Scientific Computing, v28, n5, 2006.  DOI 10.1137/050641624
!>   (Tech report version in LAWN 172 with the same title.)
!> 

Parameters

[in]	N	!> N is INTEGER !> The order of the matrix. !>
[in]	D	!> D is REAL array, dimension (N) !> The N diagonal elements of the diagonal matrix D. !>
[in]	LLD	!> LLD is REAL array, dimension (N-1) !> The (N-1) elements L(i)*L(i)*D(i). !>
[in]	SIGMA	!> SIGMA is REAL !> Shift amount in T - sigma I = L D L^T. !>
[in]	PIVMIN	!> PIVMIN is REAL !> The minimum pivot in the Sturm sequence. May be used !> when zero pivots are encountered on non-IEEE-754 !> architectures. !>
[in]	R	!> R is INTEGER !> The twist index for the twisted factorization that is used !> for the negcount. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA
Jason Riedy, University of California, Berkeley, USA

Definition at line 115 of file slaneg.f.

*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            N, R
      REAL               PIVMIN, SIGMA
*     ..
*     .. Array Arguments ..
      REAL               D( * ), LLD( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter( zero = 0.0e0, one = 1.0e0 )
*     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
      REAL               BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP
      LOGICAL SAWNAN
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC min, max
*     ..
*     .. External Functions ..
      LOGICAL SISNAN
      EXTERNAL sisnan
*     ..
*     .. 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 = sisnan( 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 (sisnan(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 = sisnan( 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 (sisnan(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
 
      slaneg = negcnt

Here is the call graph for this function:

Here is the caller graph for this function: