LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dlarrb()

subroutine dlarrb ( integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  LLD,
integer  IFIRST,
integer  ILAST,
double precision  RTOL1,
double precision  RTOL2,
integer  OFFSET,
double precision, dimension( * )  W,
double precision, dimension( * )  WGAP,
double precision, dimension( * )  WERR,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
double precision  PIVMIN,
double precision  SPDIAM,
integer  TWIST,
integer  INFO 
)

DLARRB provides limited bisection to locate eigenvalues for more accuracy.

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

Purpose:
 Given the relatively robust representation(RRR) L D L^T, DLARRB
 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.
Parameters
[in]N
          N is INTEGER
          The order of the matrix.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          The N diagonal elements of the diagonal matrix D.
[in]LLD
          LLD is DOUBLE PRECISION array, dimension (N-1)
          The (N-1) elements L(i)*L(i)*D(i).
[in]IFIRST
          IFIRST is INTEGER
          The index of the first eigenvalue to be computed.
[in]ILAST
          ILAST is INTEGER
          The index of the last eigenvalue to be computed.
[in]RTOL1
          RTOL1 is DOUBLE PRECISION
[in]RTOL2
          RTOL2 is DOUBLE PRECISION
          Tolerance for the convergence of the bisection intervals.
          An interval [LEFT,RIGHT] has converged if
          RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
          where GAP is the (estimated) distance to the nearest
          eigenvalue.
[in]OFFSET
          OFFSET is 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.
[in,out]W
          W is DOUBLE PRECISION 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.
[in,out]WGAP
          WGAP is DOUBLE PRECISION 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 = ILAST
          then WGAP(IFIRST-OFFSET) must be set to ZERO.
          On output, these gaps are refined.
[in,out]WERR
          WERR is DOUBLE PRECISION 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.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (2*N)
          Workspace.
[out]IWORK
          IWORK is INTEGER array, dimension (2*N)
          Workspace.
[in]PIVMIN
          PIVMIN is DOUBLE PRECISION
          The minimum pivot in the Sturm sequence.
[in]SPDIAM
          SPDIAM is DOUBLE PRECISION
          The spectral diameter of the matrix.
[in]TWIST
          TWIST is 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)
[out]INFO
          INFO is INTEGER
          Error flag.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

Definition at line 193 of file dlarrb.f.

196 *
197 * -- LAPACK auxiliary routine --
198 * -- LAPACK is a software package provided by Univ. of Tennessee, --
199 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
200 *
201 * .. Scalar Arguments ..
202  INTEGER IFIRST, ILAST, INFO, N, OFFSET, TWIST
203  DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPDIAM
204 * ..
205 * .. Array Arguments ..
206  INTEGER IWORK( * )
207  DOUBLE PRECISION D( * ), LLD( * ), W( * ),
208  $ WERR( * ), WGAP( * ), WORK( * )
209 * ..
210 *
211 * =====================================================================
212 *
213 * .. Parameters ..
214  DOUBLE PRECISION ZERO, TWO, HALF
215  parameter( zero = 0.0d0, two = 2.0d0,
216  $ half = 0.5d0 )
217  INTEGER MAXITR
218 * ..
219 * .. Local Scalars ..
220  INTEGER I, I1, II, IP, ITER, K, NEGCNT, NEXT, NINT,
221  $ OLNINT, PREV, R
222  DOUBLE PRECISION BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
223  $ RGAP, RIGHT, TMP, WIDTH
224 * ..
225 * .. External Functions ..
226  INTEGER DLANEG
227  EXTERNAL dlaneg
228 *
229 * ..
230 * .. Intrinsic Functions ..
231  INTRINSIC abs, max, min
232 * ..
233 * .. Executable Statements ..
234 *
235  info = 0
236 *
237 * Quick return if possible
238 *
239  IF( n.LE.0 ) THEN
240  RETURN
241  END IF
242 *
243  maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /
244  $ log( two ) ) + 2
245  mnwdth = two * pivmin
246 *
247  r = twist
248  IF((r.LT.1).OR.(r.GT.n)) r = n
249 *
250 * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
251 * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
252 * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
253 * for an unconverged interval is set to the index of the next unconverged
254 * interval, and is -1 or 0 for a converged interval. Thus a linked
255 * list of unconverged intervals is set up.
256 *
257  i1 = ifirst
258 * The number of unconverged intervals
259  nint = 0
260 * The last unconverged interval found
261  prev = 0
262 
263  rgap = wgap( i1-offset )
264  DO 75 i = i1, ilast
265  k = 2*i
266  ii = i - offset
267  left = w( ii ) - werr( ii )
268  right = w( ii ) + werr( ii )
269  lgap = rgap
270  rgap = wgap( ii )
271  gap = min( lgap, rgap )
272 
273 * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
274 * Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT
275 *
276 * Do while( NEGCNT(LEFT).GT.I-1 )
277 *
278  back = werr( ii )
279  20 CONTINUE
280  negcnt = dlaneg( n, d, lld, left, pivmin, r )
281  IF( negcnt.GT.i-1 ) THEN
282  left = left - back
283  back = two*back
284  GO TO 20
285  END IF
286 *
287 * Do while( NEGCNT(RIGHT).LT.I )
288 * Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
289 *
290  back = werr( ii )
291  50 CONTINUE
292 
293  negcnt = dlaneg( n, d, lld, right, pivmin, r )
294  IF( negcnt.LT.i ) THEN
295  right = right + back
296  back = two*back
297  GO TO 50
298  END IF
299  width = half*abs( left - right )
300  tmp = max( abs( left ), abs( right ) )
301  cvrgd = max(rtol1*gap,rtol2*tmp)
302  IF( width.LE.cvrgd .OR. width.LE.mnwdth ) THEN
303 * This interval has already converged and does not need refinement.
304 * (Note that the gaps might change through refining the
305 * eigenvalues, however, they can only get bigger.)
306 * Remove it from the list.
307  iwork( k-1 ) = -1
308 * Make sure that I1 always points to the first unconverged interval
309  IF((i.EQ.i1).AND.(i.LT.ilast)) i1 = i + 1
310  IF((prev.GE.i1).AND.(i.LE.ilast)) iwork( 2*prev-1 ) = i + 1
311  ELSE
312 * unconverged interval found
313  prev = i
314  nint = nint + 1
315  iwork( k-1 ) = i + 1
316  iwork( k ) = negcnt
317  END IF
318  work( k-1 ) = left
319  work( k ) = right
320  75 CONTINUE
321 
322 *
323 * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
324 * and while (ITER.LT.MAXITR)
325 *
326  iter = 0
327  80 CONTINUE
328  prev = i1 - 1
329  i = i1
330  olnint = nint
331 
332  DO 100 ip = 1, olnint
333  k = 2*i
334  ii = i - offset
335  rgap = wgap( ii )
336  lgap = rgap
337  IF(ii.GT.1) lgap = wgap( ii-1 )
338  gap = min( lgap, rgap )
339  next = iwork( k-1 )
340  left = work( k-1 )
341  right = work( k )
342  mid = half*( left + right )
343 
344 * semiwidth of interval
345  width = right - mid
346  tmp = max( abs( left ), abs( right ) )
347  cvrgd = max(rtol1*gap,rtol2*tmp)
348  IF( ( width.LE.cvrgd ) .OR. ( width.LE.mnwdth ).OR.
349  $ ( iter.EQ.maxitr ) )THEN
350 * reduce number of unconverged intervals
351  nint = nint - 1
352 * Mark interval as converged.
353  iwork( k-1 ) = 0
354  IF( i1.EQ.i ) THEN
355  i1 = next
356  ELSE
357 * Prev holds the last unconverged interval previously examined
358  IF(prev.GE.i1) iwork( 2*prev-1 ) = next
359  END IF
360  i = next
361  GO TO 100
362  END IF
363  prev = i
364 *
365 * Perform one bisection step
366 *
367  negcnt = dlaneg( n, d, lld, mid, pivmin, r )
368  IF( negcnt.LE.i-1 ) THEN
369  work( k-1 ) = mid
370  ELSE
371  work( k ) = mid
372  END IF
373  i = next
374  100 CONTINUE
375  iter = iter + 1
376 * do another loop if there are still unconverged intervals
377 * However, in the last iteration, all intervals are accepted
378 * since this is the best we can do.
379  IF( ( nint.GT.0 ).AND.(iter.LE.maxitr) ) GO TO 80
380 *
381 *
382 * At this point, all the intervals have converged
383  DO 110 i = ifirst, ilast
384  k = 2*i
385  ii = i - offset
386 * All intervals marked by '0' have been refined.
387  IF( iwork( k-1 ).EQ.0 ) THEN
388  w( ii ) = half*( work( k-1 )+work( k ) )
389  werr( ii ) = work( k ) - w( ii )
390  END IF
391  110 CONTINUE
392 *
393  DO 111 i = ifirst+1, ilast
394  k = 2*i
395  ii = i - offset
396  wgap( ii-1 ) = max( zero,
397  $ w(ii) - werr(ii) - w( ii-1 ) - werr( ii-1 ))
398  111 CONTINUE
399 
400  RETURN
401 *
402 * End of DLARRB
403 *
integer function dlaneg(N, D, LLD, SIGMA, PIVMIN, R)
DLANEG computes the Sturm count.
Definition: dlaneg.f:118
Here is the caller graph for this function: