slarrb2.f

Go to the documentation of this file.

      SUBROUTINE slarrb2( N, D, LLD, IFIRST, ILAST, RTOL1,
     $                   RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
     $                   PIVMIN, LGPVMN, LGSPDM, TWIST, INFO )
*
*  -- ScaLAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver
*     July 4, 2010
*
      IMPLICIT NONE
*
*     .. Scalar Arguments ..
      INTEGER            IFIRST, ILAST, INFO, N, OFFSET, TWIST
      REAL               LGPVMN, LGSPDM, PIVMIN, 
     $                   rtol1, rtol2
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      REAL               D( * ), LLD( * ), W( * ),
     $                   werr( * ), wgap( * ), work( * )
*     ..
*
*  Purpose
*  =======
*
*  Given the relatively robust representation(RRR) L D L^T, SLARRB2
*  does "limited" bisection to refine the eigenvalues of L D L^T,
*  W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
*  guesses for these eigenvalues are input in W, the corresponding estimate
*  of the error in these guesses and their gaps are input in WERR
*  and WGAP, respectively. During bisection, intervals
*  [left, right] are maintained by storing their mid-points and
*  semi-widths in the arrays W and WERR respectively.
*
*  NOTE: 
*  There are very few minor differences between SLARRB from LAPACK
*  and this current subroutine SLARRB2.
*  The most important reason for creating this nearly identical copy
*  is profiling: in the ScaLAPACK MRRR algorithm, eigenvalue computation 
*  using SLARRB2 is used for refinement in the construction of 
*  the representation tree, as opposed to the initial computation of the
*  eigenvalues for the root RRR which uses SLARRB. When profiling,
*  this allows an easy quantification of refinement work vs. computing
*  eigenvalues of the root.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The order of the matrix.
*
*  D       (input) REAL             array, dimension (N)
*          The N diagonal elements of the diagonal matrix D.
*
*  LLD     (input) REAL             array, dimension (N-1)
*          The (N-1) elements L(i)*L(i)*D(i).
*
*  IFIRST  (input) INTEGER
*          The index of the first eigenvalue to be computed.
*
*  ILAST   (input) INTEGER
*          The index of the last eigenvalue to be computed.
*
*  RTOL1   (input) REAL            
*  RTOL2   (input) REAL            
*          Tolerance for the convergence of the bisection intervals.
*          An interval [LEFT,RIGHT] has converged if
*          RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
*          where GAP is the (estimated) distance to the nearest
*          eigenvalue.
*
*  OFFSET  (input) INTEGER
*          Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET
*          through ILAST-OFFSET elements of these arrays are to be used.
*
*  W       (input/output) REAL             array, dimension (N)
*          On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
*          estimates of the eigenvalues of L D L^T indexed IFIRST through ILAST.
*          On output, these estimates are refined.
*
*  WGAP    (input/output) REAL             array, dimension (N-1)
*          On input, the (estimated) gaps between consecutive
*          eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between
*          eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST
*          then WGAP(IFIRST-OFFSET) must be set to ZERO.
*          On output, these gaps are refined.
*
*  WERR    (input/output) REAL             array, dimension (N)
*          On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
*          the errors in the estimates of the corresponding elements in W.
*          On output, these errors are refined.
*
*  WORK    (workspace) REAL             array, dimension (4*N)
*          Workspace.
*
*  IWORK   (workspace) INTEGER array, dimension (2*N)
*          Workspace.
*
*  PIVMIN  (input) REAL             
*          The minimum pivot in the sturm sequence.
*
*  LGPVMN  (input) REAL            
*          Logarithm of PIVMIN, precomputed.
*
*  LGSPDM  (input) REAL             
*          Logarithm of the spectral diameter, precomputed.
*
*  TWIST   (input) INTEGER
*          The twist index for the twisted factorization that is used
*          for the negcount. 
*          TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T
*          TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T
*          TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r)
*
*  INFO    (output) INTEGER
*          Error flag.
*
*     .. Parameters ..
      REAL               ZERO, TWO, HALF
      PARAMETER        ( ZERO = 0.0e0, two = 2.0e0,
     $                   half = 0.5e0 )
      INTEGER   MAXITR
*     ..
*     .. Local Scalars ..
      INTEGER            I, I1, II, INDLLD, IP, ITER, J, K, NEGCNT,
     $                   NEXT, NINT, OLNINT, PREV, R
      REAL               BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
     $                   RGAP, RIGHT, SAVGAP, TMP, WIDTH
      LOGICAL   PARANOID
*     ..
*     .. External Functions ..
      LOGICAL            SISNAN
      REAL               SLAMCH
      INTEGER            SLANEG2A
      EXTERNAL           sisnan, slamch, 
     $                   slaneg2a
*
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min
*     ..
*     .. Executable Statements ..
*
      info = 0
*     
*     Turn on paranoid check for rounding errors 
*     invalidating uncertainty intervals of eigenvalues
*
      paranoid = .true.
*
      maxitr = int( ( lgspdm - lgpvmn ) / log( two ) ) + 2
      mnwdth = two * pivmin
*
      r = twist
*
      indlld = 2*n     
      DO 5 j = 1, n-1 
         i=2*j
         work(indlld+i-1) = d(j)
         work(indlld+i) = lld(j)
  5   CONTINUE
      work(indlld+2*n-1) = d(n)
*
      IF((r.LT.1).OR.(r.GT.n)) r = n
*
*     Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
*     The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
*     Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
*     for an unconverged interval is set to the index of the next unconverged
*     interval, and is -1 or 0 for a converged interval. Thus a linked
*     list of unconverged intervals is set up.
*
      i1 = ifirst
*     The number of unconverged intervals 
      nint = 0
*     The last unconverged interval found
      prev = 0
     
      rgap = wgap( i1-offset )
      DO 75 i = i1, ilast
         k = 2*i
         ii = i - offset
         left = w( ii ) - werr( ii )
         right = w( ii ) + werr( ii )
         lgap = rgap
         rgap = wgap( ii )
         gap = min( lgap, rgap )
 
         IF((abs(left).LE.16*pivmin).OR.(abs(right).LE.16*pivmin))
     $      THEN
            info = -1
            RETURN
         ENDIF
 
         IF( paranoid ) THEN
*        Make sure that [LEFT,RIGHT] contains the desired eigenvalue
*        Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT 
*        
*        Do while( NEGCNT(LEFT).GT.I-1 )
*        
         back = werr( ii )
 20      CONTINUE
         negcnt = slaneg2a( n, work(indlld+1), left, pivmin, r )
         IF( negcnt.GT.i-1 ) THEN
            left = left - back
            back = two*back
            GO TO 20
         END IF
*
*        Do while( NEGCNT(RIGHT).LT.I )
*        Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT 
*        
         back = werr( ii )
 50      CONTINUE
         negcnt = slaneg2a( n, work(indlld+1),right, pivmin, r )
 
         IF( negcnt.LT.i ) THEN
             right = right + back
             back = two*back
             GO TO 50
         END IF
         ENDIF
 
         width = half*abs( left - right )
         tmp = max( abs( left ), abs( right ) )
         cvrgd = max(rtol1*gap,rtol2*tmp)
         IF( width.LE.cvrgd .OR. width.LE.mnwdth ) THEN
*           This interval has already converged and does not need refinement.
*           (Note that the gaps might change through refining the 
*            eigenvalues, however, they can only get bigger.)
*           Remove it from the list.
            iwork( k-1 ) = -1
*           Make sure that I1 always points to the first unconverged interval
            IF((i.EQ.i1).AND.(i.LT.ilast)) i1 = i + 1
            IF((prev.GE.i1).AND.(i.LE.ilast)) iwork( 2*prev-1 ) = i + 1
         ELSE
*           unconverged interval found
            prev = i
            nint = nint + 1
            iwork( k-1 ) = i + 1
            iwork( k ) = negcnt
         END IF
         work( k-1 ) = left
         work( k ) = right
 75   CONTINUE
 
*
*     Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
*     and while (ITER.LT.MAXITR)
*
      iter = 0 
 80   CONTINUE
      prev = i1 - 1
      i = i1
      olnint = nint
 
      DO 100 ip = 1, olnint
         k = 2*i
         ii = i - offset
         rgap = wgap( ii )
         lgap = rgap
         IF(ii.GT.1) lgap = wgap( ii-1 ) 
         gap = min( lgap, rgap )
         next = iwork( k-1 )
         left = work( k-1 )
         right = work( k )
         mid = half*( left + right ) 
*        semiwidth of interval
         width = right - mid
         tmp = max( abs( left ), abs( right ) )
         cvrgd = max(rtol1*gap,rtol2*tmp)
         IF( ( width.LE.cvrgd ) .OR. ( width.LE.mnwdth ).OR.
     $       ( iter.EQ.maxitr ) )THEN
*           reduce number of unconverged intervals
            nint = nint - 1
*           Mark interval as converged. 
            iwork( k-1 ) = 0
            IF( i1.EQ.i ) THEN
               i1 = next
            ELSE
*              Prev holds the last unconverged interval previously examined
               IF(prev.GE.i1) iwork( 2*prev-1 ) = next
            END IF
            i = next
            GO TO 100
         END IF
         prev = i
*
*        Perform one bisection step
*
         negcnt = slaneg2a( n, work(indlld+1), mid, pivmin, r )
         IF( negcnt.LE.i-1 ) THEN
            work( k-1 ) = mid
         ELSE
            work( k ) = mid
         END IF
         i = next
 100  CONTINUE
      iter = iter + 1
*     do another loop if there are still unconverged intervals
*     However, in the last iteration, all intervals are accepted
*     since this is the best we can do.
      IF( ( nint.GT.0 ).AND.(iter.LE.maxitr) ) GO TO 80
*
*
*     At this point, all the intervals have converged
*
*     save this gap to restore it after the loop
      savgap = wgap( ilast-offset )
*
      left = work( 2*ifirst-1 )
      DO 110 i = ifirst, ilast
         k = 2*i
         ii = i - offset
*        RIGHT is the right boundary of this current interval
         right = work( k ) 
*        All intervals marked by '0' have been refined.
         IF( iwork( k-1 ).EQ.0 ) THEN
            w( ii ) = half*( left+right )
            werr( ii ) = right - w( ii )
         END IF
*        Left is the boundary of the next interval
         left = work( k +1 ) 
         wgap( ii ) = max( zero, left - right )
 110  CONTINUE
*     restore the last gap which was overwritten by garbage
      wgap( ilast-offset ) = savgap
 
      RETURN
*
*     End of SLARRB2
*
      END
*
*
*
      FUNCTION slaneg2( N, D, LLD, SIGMA, PIVMIN, R )
*
      IMPLICIT NONE
*
      INTEGER slaneg2
*
*     .. Scalar Arguments ..
      INTEGER            n, r
      REAL               pivmin, sigma
*     ..
*     .. Array Arguments ..
      REAL               d( * ), lld( * )
*
      REAL               zero
      PARAMETER        ( zero = 0.0e0 )
 
      INTEGER blklen
      PARAMETER ( blklen = 2048 )
*     ..
*     .. Local Scalars ..
      INTEGER            bj, j, neg1, neg2, negcnt, to
      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 
*
      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 = sisnan( 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 = sisnan( 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
 
      slaneg2 = negcnt
      END
*
*
*
      FUNCTION slaneg2a( N, DLLD, SIGMA, PIVMIN, R )
*
      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
      END