dlaneg2()

integer function dlaneg2	(	integer	n,
		double precision, dimension( * )	d,
		double precision, dimension( * )	lld,
		double precision	sigma,
		double precision	pivmin,
		integer	r
	)

Definition at line 336 of file dlarrb2.f.

*
      IMPLICIT NONE
*
      INTEGER DLANEG2
*
*     .. Scalar Arguments ..
      INTEGER            N, R
      DOUBLE PRECISION   PIVMIN, SIGMA
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   D( * ), LLD( * )
*
      DOUBLE PRECISION   ZERO
      parameter( zero = 0.0d0 )
 
      INTEGER BLKLEN
      parameter( blklen = 2048 )
*     ..
*     .. Local Scalars ..
      INTEGER            BJ, J, NEG1, NEG2, NEGCNT, TO
      DOUBLE PRECISION   DMINUS, DPLUS, GAMMA, P, S, T, TMP, XSAV
      LOGICAL SAWNAN
*     ..
*     .. External Functions ..
      LOGICAL DISNAN
      EXTERNAL disnan
      
      negcnt = 0
*      
*     I) upper part: L D L^T - SIGMA I = L+ D+ L+^T
*     run dstqds block-wise to avoid excessive work when NaNs occur 
*
      s = zero
      DO 210 bj = 1, r-1, blklen
         neg1 = 0
         xsav = s
         to = bj+blklen-1 
         IF ( to.LE.r-1 ) THEN
            DO 21 j = bj, to
               t = s - sigma
               dplus = d( j ) + t
               IF( dplus.LT.zero ) neg1=neg1 + 1
               s = t*lld( j ) / dplus 
 21         CONTINUE
         ELSE
            DO 22 j = bj, r-1
               t = s - sigma
               dplus = d( j ) + t
               IF( dplus.LT.zero ) neg1=neg1 + 1
               s = t*lld( j ) / dplus 
 22         CONTINUE
         ENDIF
         sawnan = disnan( s )
*         
         IF( sawnan ) THEN
            neg1 = 0
            s = xsav
            to = bj+blklen-1 
            IF ( to.LE.r-1 ) THEN
               DO 23 j = bj, to
                  t = s - sigma
                  dplus = d( j ) + t
                  IF(abs(dplus).LT.pivmin) 
     $               dplus = -pivmin
                  tmp = lld( j ) / dplus
                  IF( dplus.LT.zero ) 
     $               neg1 = neg1 + 1
                  s = t*tmp
                  IF( tmp.EQ.zero ) s = lld( j )
 23            CONTINUE
            ELSE
               DO 24 j = bj, r-1
                  t = s - sigma
                  dplus = d( j ) + t
                  IF(abs(dplus).LT.pivmin) 
     $               dplus = -pivmin
                  tmp = lld( j ) / dplus
                  IF( dplus.LT.zero ) neg1=neg1+1
                  s = t*tmp
                  IF( tmp.EQ.zero ) s = lld( j )
 24            CONTINUE
            ENDIF
         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
         xsav = p
         to = bj-blklen+1
         IF ( to.GE.r ) THEN
            DO 25 j = bj, to, -1
               dminus = lld( j ) + p
               IF( dminus.LT.zero ) neg2=neg2+1
               tmp = p / dminus
               p = tmp * d( j ) - sigma
 25         CONTINUE
         ELSE
            DO 26 j = bj, r, -1
               dminus = lld( j ) + p
               IF( dminus.LT.zero ) neg2=neg2+1
               tmp = p / dminus
               p = tmp * d( j ) - sigma
 26         CONTINUE
         ENDIF
         sawnan = disnan( p )
*
         IF( sawnan ) THEN
            neg2 = 0
            p = xsav
            to = bj-blklen+1
            IF ( to.GE.r ) THEN
               DO 27 j = bj, to, -1
                  dminus = lld( j ) + p
                  IF(abs(dminus).LT.pivmin) 
     $               dminus = -pivmin
                  tmp = d( j ) / dminus
                  IF( dminus.LT.zero ) 
     $               neg2 = neg2 + 1
                  p = p*tmp - sigma
                  IF( tmp.EQ.zero ) 
     $               p = d( j ) - sigma
 27            CONTINUE
            ELSE
               DO 28 j = bj, r, -1
                  dminus = lld( j ) + p
                  IF(abs(dminus).LT.pivmin) 
     $               dminus = -pivmin
                  tmp = d( j ) / dminus
                  IF( dminus.LT.zero ) 
     $               neg2 = neg2 + 1
                  p = p*tmp - sigma
                  IF( tmp.EQ.zero ) 
     $               p = d( j ) - sigma
 28            CONTINUE
            ENDIF
         END IF
         negcnt = negcnt + neg2
 230  CONTINUE
*     
*     III) Twist index
*
      gamma = s + p
      IF( gamma.LT.zero ) negcnt = negcnt+1
 
      dlaneg2 = negcnt

Here is the call graph for this function:

Here is the caller graph for this function: