slaneg2a()

integer function slaneg2a	(	integer	n,
		real, dimension( * )	dlld,
		real	sigma,
		real	pivmin,
		integer	r
	)

Definition at line 490 of file slarrb2.f.

*
      IMPLICIT NONE
*
      INTEGER SLANEG2A
*
*     .. Scalar Arguments ..
      INTEGER            N, R
      REAL               PIVMIN, SIGMA
*     ..
*     .. Array Arguments ..
      REAL               DLLD( * )
*
      REAL               ZERO
      parameter( zero = 0.0e0 )
 
      INTEGER BLKLEN
      parameter( blklen = 512 )
*
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          int
*     ..
*     .. Local Scalars ..
      INTEGER            BJ, I, J, NB, NEG1, NEG2, NEGCNT, NX
      REAL               DMINUS, DPLUS, GAMMA, P, S, T, TMP, XSAV
      LOGICAL SAWNAN
*     ..
*     .. External Functions ..
      LOGICAL SISNAN
      EXTERNAL sisnan
      
      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, 
*     first in chunks of size BLKLEN and then the remainder
*
      nb = int((r-1)/blklen)
      nx = nb*blklen
      s = zero      
      DO 210 bj = 1, nx, blklen
         neg1 = 0
         xsav = s
         DO 21 j = bj, bj+blklen-1 
            i = 2*j
            t = s - sigma
            dplus = dlld( i-1 ) + t
            IF( dplus.LT.zero ) neg1=neg1 + 1
            s = t*dlld( i ) / dplus 
 21      CONTINUE
         sawnan = sisnan( s )
*         
         IF( sawnan ) THEN
            neg1 = 0
            s = xsav
            DO 23 j = bj, bj+blklen-1 
               i = 2*j
               t = s - sigma
               dplus = dlld( i-1 ) + t
               IF(abs(dplus).LT.pivmin) 
     $            dplus = -pivmin
               tmp = dlld( i ) / dplus
               IF( dplus.LT.zero ) 
     $            neg1 = neg1 + 1
               s = t*tmp
               IF( tmp.EQ.zero ) s = dlld( i )
 23         CONTINUE
         END IF
         negcnt = negcnt + neg1
 210  CONTINUE
*
      neg1 = 0
      xsav = s
      DO 22 j = nx+1, r-1
         i = 2*j
         t = s - sigma
         dplus = dlld( i-1 ) + t
         IF( dplus.LT.zero ) neg1=neg1 + 1
         s = t*dlld( i ) / dplus 
 22   CONTINUE
      sawnan = sisnan( s )
*         
      IF( sawnan ) THEN
         neg1 = 0
         s = xsav
         DO 24 j = nx+1, r-1
            i = 2*j
            t = s - sigma
            dplus = dlld( i-1 ) + t
            IF(abs(dplus).LT.pivmin) 
     $         dplus = -pivmin
            tmp = dlld( i ) / dplus
            IF( dplus.LT.zero ) neg1=neg1+1
            s = t*tmp
            IF( tmp.EQ.zero ) s = dlld( i )
 24      CONTINUE
      ENDIF
      negcnt = negcnt + neg1
*
*     II) lower part: L D L^T - SIGMA I = U- D- U-^T
*     
      nb = int((n-r)/blklen)
      nx = n-nb*blklen
      p = dlld( 2*n-1 ) - sigma
      DO 230 bj = n-1, nx, -blklen
         neg2 = 0
         xsav = p
         DO 25 j = bj, bj-blklen+1, -1
            i = 2*j
            dminus = dlld( i ) + p
            IF( dminus.LT.zero ) neg2=neg2+1
            tmp = p / dminus
            p = tmp * dlld( i-1 ) - sigma
 25      CONTINUE
         sawnan = sisnan( p )
*
         IF( sawnan ) THEN
            neg2 = 0
            p = xsav
            DO 27 j = bj, bj-blklen+1, -1
               i = 2*j
               dminus = dlld( i ) + p
               IF(abs(dminus).LT.pivmin) 
     $            dminus = -pivmin
               tmp = dlld( i-1 ) / dminus
               IF( dminus.LT.zero ) 
     $            neg2 = neg2 + 1
               p = p*tmp - sigma
               IF( tmp.EQ.zero ) 
     $            p = dlld( i-1 ) - sigma
 27         CONTINUE
         END IF
         negcnt = negcnt + neg2
 230  CONTINUE
 
      neg2 = 0
      xsav = p
      DO 26 j = nx-1, r, -1
         i = 2*j
         dminus = dlld( i ) + p
         IF( dminus.LT.zero ) neg2=neg2+1
         tmp = p / dminus
         p = tmp * dlld( i-1 ) - sigma
 26   CONTINUE
      sawnan = sisnan( p )
*
      IF( sawnan ) THEN
         neg2 = 0
         p = xsav
         DO 28 j = nx-1, r, -1
            i = 2*j
            dminus = dlld( i ) + p
            IF(abs(dminus).LT.pivmin) 
     $         dminus = -pivmin
            tmp = dlld( i-1 ) / dminus
            IF( dminus.LT.zero ) 
     $         neg2 = neg2 + 1
            p = p*tmp - sigma
            IF( tmp.EQ.zero ) 
     $         p = dlld( i-1 ) - sigma
 28      CONTINUE
      END IF
      negcnt = negcnt + neg2
*     
*     III) Twist index
*
      gamma = s + p
      IF( gamma.LT.zero ) negcnt = negcnt+1
 
      slaneg2a = negcnt

Here is the call graph for this function:

Here is the caller graph for this function: