LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
dlarrd.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine dlarrd (RANGE, ORDER, N, VL, VU, IL, IU, GERS, RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, M, W, WERR, WL, WU, IBLOCK, INDEXW, WORK, IWORK, INFO)
 DLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy. More...
 

Function/Subroutine Documentation

subroutine dlarrd ( character  RANGE,
character  ORDER,
integer  N,
double precision  VL,
double precision  VU,
integer  IL,
integer  IU,
double precision, dimension( * )  GERS,
double precision  RELTOL,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( * )  E2,
double precision  PIVMIN,
integer  NSPLIT,
integer, dimension( * )  ISPLIT,
integer  M,
double precision, dimension( * )  W,
double precision, dimension( * )  WERR,
double precision  WL,
double precision  WU,
integer, dimension( * )  IBLOCK,
integer, dimension( * )  INDEXW,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy.

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

Purpose:
 DLARRD computes the eigenvalues of a symmetric tridiagonal
 matrix T to suitable accuracy. This is an auxiliary code to be
 called from DSTEMR.
 The user may ask for all eigenvalues, all eigenvalues
 in the half-open interval (VL, VU], or the IL-th through IU-th
 eigenvalues.

 To avoid overflow, the matrix must be scaled so that its
 largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
 accuracy, it should not be much smaller than that.

 See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
 Matrix", Report CS41, Computer Science Dept., Stanford
 University, July 21, 1966.
Parameters
[in]RANGE
          RANGE is CHARACTER*1
          = 'A': ("All")   all eigenvalues will be found.
          = 'V': ("Value") all eigenvalues in the half-open interval
                           (VL, VU] will be found.
          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
                           entire matrix) will be found.
[in]ORDER
          ORDER is CHARACTER*1
          = 'B': ("By Block") the eigenvalues will be grouped by
                              split-off block (see IBLOCK, ISPLIT) and
                              ordered from smallest to largest within
                              the block.
          = 'E': ("Entire matrix")
                              the eigenvalues for the entire matrix
                              will be ordered from smallest to
                              largest.
[in]N
          N is INTEGER
          The order of the tridiagonal matrix T.  N >= 0.
[in]VL
          VL is DOUBLE PRECISION
[in]VU
          VU is DOUBLE PRECISION
          If RANGE='V', the lower and upper bounds of the interval to
          be searched for eigenvalues.  Eigenvalues less than or equal
          to VL, or greater than VU, will not be returned.  VL < VU.
          Not referenced if RANGE = 'A' or 'I'.
[in]IL
          IL is INTEGER
[in]IU
          IU is INTEGER
          If RANGE='I', the indices (in ascending order) of the
          smallest and largest eigenvalues to be returned.
          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
          Not referenced if RANGE = 'A' or 'V'.
[in]GERS
          GERS is DOUBLE PRECISION array, dimension (2*N)
          The N Gerschgorin intervals (the i-th Gerschgorin interval
          is (GERS(2*i-1), GERS(2*i)).
[in]RELTOL
          RELTOL is DOUBLE PRECISION
          The minimum relative width of an interval.  When an interval
          is narrower than RELTOL times the larger (in
          magnitude) endpoint, then it is considered to be
          sufficiently small, i.e., converged.  Note: this should
          always be at least radix*machine epsilon.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          The n diagonal elements of the tridiagonal matrix T.
[in]E
          E is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) off-diagonal elements of the tridiagonal matrix T.
[in]E2
          E2 is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
[in]PIVMIN
          PIVMIN is DOUBLE PRECISION
          The minimum pivot allowed in the Sturm sequence for T.
[in]NSPLIT
          NSPLIT is INTEGER
          The number of diagonal blocks in the matrix T.
          1 <= NSPLIT <= N.
[in]ISPLIT
          ISPLIT is INTEGER array, dimension (N)
          The splitting points, at which T breaks up into submatrices.
          The first submatrix consists of rows/columns 1 to ISPLIT(1),
          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
          etc., and the NSPLIT-th consists of rows/columns
          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
          (Only the first NSPLIT elements will actually be used, but
          since the user cannot know a priori what value NSPLIT will
          have, N words must be reserved for ISPLIT.)
[out]M
          M is INTEGER
          The actual number of eigenvalues found. 0 <= M <= N.
          (See also the description of INFO=2,3.)
[out]W
          W is DOUBLE PRECISION array, dimension (N)
          On exit, the first M elements of W will contain the
          eigenvalue approximations. DLARRD computes an interval
          I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue
          approximation is given as the interval midpoint
          W(j)= ( a_j + b_j)/2. The corresponding error is bounded by
          WERR(j) = abs( a_j - b_j)/2
[out]WERR
          WERR is DOUBLE PRECISION array, dimension (N)
          The error bound on the corresponding eigenvalue approximation
          in W.
[out]WL
          WL is DOUBLE PRECISION
[out]WU
          WU is DOUBLE PRECISION
          The interval (WL, WU] contains all the wanted eigenvalues.
          If RANGE='V', then WL=VL and WU=VU.
          If RANGE='A', then WL and WU are the global Gerschgorin bounds
                        on the spectrum.
          If RANGE='I', then WL and WU are computed by DLAEBZ from the
                        index range specified.
[out]IBLOCK
          IBLOCK is INTEGER array, dimension (N)
          At each row/column j where E(j) is zero or small, the
          matrix T is considered to split into a block diagonal
          matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which
          block (from 1 to the number of blocks) the eigenvalue W(i)
          belongs.  (DLARRD may use the remaining N-M elements as
          workspace.)
[out]INDEXW
          INDEXW is INTEGER array, dimension (N)
          The indices of the eigenvalues within each block (submatrix);
          for example, INDEXW(i)= j and IBLOCK(i)=k imply that the
          i-th eigenvalue W(i) is the j-th eigenvalue in block k.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (4*N)
[out]IWORK
          IWORK is INTEGER array, dimension (3*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  some or all of the eigenvalues failed to converge or
                were not computed:
                =1 or 3: Bisection failed to converge for some
                        eigenvalues; these eigenvalues are flagged by a
                        negative block number.  The effect is that the
                        eigenvalues may not be as accurate as the
                        absolute and relative tolerances.  This is
                        generally caused by unexpectedly inaccurate
                        arithmetic.
                =2 or 3: RANGE='I' only: Not all of the eigenvalues
                        IL:IU were found.
                        Effect: M < IU+1-IL
                        Cause:  non-monotonic arithmetic, causing the
                                Sturm sequence to be non-monotonic.
                        Cure:   recalculate, using RANGE='A', and pick
                                out eigenvalues IL:IU.  In some cases,
                                increasing the PARAMETER "FUDGE" may
                                make things work.
                = 4:    RANGE='I', and the Gershgorin interval
                        initially used was too small.  No eigenvalues
                        were computed.
                        Probable cause: your machine has sloppy
                                        floating-point arithmetic.
                        Cure: Increase the PARAMETER "FUDGE",
                              recompile, and try again.
Internal Parameters:
  FUDGE   DOUBLE PRECISION, default = 2
          A "fudge factor" to widen the Gershgorin intervals.  Ideally,
          a value of 1 should work, but on machines with sloppy
          arithmetic, this needs to be larger.  The default for
          publicly released versions should be large enough to handle
          the worst machine around.  Note that this has no effect
          on accuracy of the solution.
Contributors:
W. Kahan, University of California, Berkeley, USA
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
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 323 of file dlarrd.f.

323 *
324 * -- LAPACK auxiliary routine (version 3.4.2) --
325 * -- LAPACK is a software package provided by Univ. of Tennessee, --
326 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
327 * September 2012
328 *
329 * .. Scalar Arguments ..
330  CHARACTER order, range
331  INTEGER il, info, iu, m, n, nsplit
332  DOUBLE PRECISION pivmin, reltol, vl, vu, wl, wu
333 * ..
334 * .. Array Arguments ..
335  INTEGER iblock( * ), indexw( * ),
336  $ isplit( * ), iwork( * )
337  DOUBLE PRECISION d( * ), e( * ), e2( * ),
338  $ gers( * ), w( * ), werr( * ), work( * )
339 * ..
340 *
341 * =====================================================================
342 *
343 * .. Parameters ..
344  DOUBLE PRECISION zero, one, two, half, fudge
345  parameter( zero = 0.0d0, one = 1.0d0,
346  $ two = 2.0d0, half = one/two,
347  $ fudge = two )
348  INTEGER allrng, valrng, indrng
349  parameter( allrng = 1, valrng = 2, indrng = 3 )
350 * ..
351 * .. Local Scalars ..
352  LOGICAL ncnvrg, toofew
353  INTEGER i, ib, ibegin, idiscl, idiscu, ie, iend, iinfo,
354  $ im, in, ioff, iout, irange, itmax, itmp1,
355  $ itmp2, iw, iwoff, j, jblk, jdisc, je, jee, nb,
356  $ nwl, nwu
357  DOUBLE PRECISION atoli, eps, gl, gu, rtoli, tmp1, tmp2,
358  $ tnorm, uflow, wkill, wlu, wul
359 
360 * ..
361 * .. Local Arrays ..
362  INTEGER idumma( 1 )
363 * ..
364 * .. External Functions ..
365  LOGICAL lsame
366  INTEGER ilaenv
367  DOUBLE PRECISION dlamch
368  EXTERNAL lsame, ilaenv, dlamch
369 * ..
370 * .. External Subroutines ..
371  EXTERNAL dlaebz
372 * ..
373 * .. Intrinsic Functions ..
374  INTRINSIC abs, int, log, max, min
375 * ..
376 * .. Executable Statements ..
377 *
378  info = 0
379 *
380 * Decode RANGE
381 *
382  IF( lsame( range, 'A' ) ) THEN
383  irange = allrng
384  ELSE IF( lsame( range, 'V' ) ) THEN
385  irange = valrng
386  ELSE IF( lsame( range, 'I' ) ) THEN
387  irange = indrng
388  ELSE
389  irange = 0
390  END IF
391 *
392 * Check for Errors
393 *
394  IF( irange.LE.0 ) THEN
395  info = -1
396  ELSE IF( .NOT.(lsame(order,'B').OR.lsame(order,'E')) ) THEN
397  info = -2
398  ELSE IF( n.LT.0 ) THEN
399  info = -3
400  ELSE IF( irange.EQ.valrng ) THEN
401  IF( vl.GE.vu )
402  $ info = -5
403  ELSE IF( irange.EQ.indrng .AND.
404  $ ( il.LT.1 .OR. il.GT.max( 1, n ) ) ) THEN
405  info = -6
406  ELSE IF( irange.EQ.indrng .AND.
407  $ ( iu.LT.min( n, il ) .OR. iu.GT.n ) ) THEN
408  info = -7
409  END IF
410 *
411  IF( info.NE.0 ) THEN
412  RETURN
413  END IF
414 
415 * Initialize error flags
416  info = 0
417  ncnvrg = .false.
418  toofew = .false.
419 
420 * Quick return if possible
421  m = 0
422  IF( n.EQ.0 ) RETURN
423 
424 * Simplification:
425  IF( irange.EQ.indrng .AND. il.EQ.1 .AND. iu.EQ.n ) irange = 1
426 
427 * Get machine constants
428  eps = dlamch( 'P' )
429  uflow = dlamch( 'U' )
430 
431 
432 * Special Case when N=1
433 * Treat case of 1x1 matrix for quick return
434  IF( n.EQ.1 ) THEN
435  IF( (irange.EQ.allrng).OR.
436  $ ((irange.EQ.valrng).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
437  $ ((irange.EQ.indrng).AND.(il.EQ.1).AND.(iu.EQ.1)) ) THEN
438  m = 1
439  w(1) = d(1)
440 * The computation error of the eigenvalue is zero
441  werr(1) = zero
442  iblock( 1 ) = 1
443  indexw( 1 ) = 1
444  ENDIF
445  RETURN
446  END IF
447 
448 * NB is the minimum vector length for vector bisection, or 0
449 * if only scalar is to be done.
450  nb = ilaenv( 1, 'DSTEBZ', ' ', n, -1, -1, -1 )
451  IF( nb.LE.1 ) nb = 0
452 
453 * Find global spectral radius
454  gl = d(1)
455  gu = d(1)
456  DO 5 i = 1,n
457  gl = min( gl, gers( 2*i - 1))
458  gu = max( gu, gers(2*i) )
459  5 CONTINUE
460 * Compute global Gerschgorin bounds and spectral diameter
461  tnorm = max( abs( gl ), abs( gu ) )
462  gl = gl - fudge*tnorm*eps*n - fudge*two*pivmin
463  gu = gu + fudge*tnorm*eps*n + fudge*two*pivmin
464 * [JAN/28/2009] remove the line below since SPDIAM variable not use
465 * SPDIAM = GU - GL
466 * Input arguments for DLAEBZ:
467 * The relative tolerance. An interval (a,b] lies within
468 * "relative tolerance" if b-a < RELTOL*max(|a|,|b|),
469  rtoli = reltol
470 * Set the absolute tolerance for interval convergence to zero to force
471 * interval convergence based on relative size of the interval.
472 * This is dangerous because intervals might not converge when RELTOL is
473 * small. But at least a very small number should be selected so that for
474 * strongly graded matrices, the code can get relatively accurate
475 * eigenvalues.
476  atoli = fudge*two*uflow + fudge*two*pivmin
477 
478  IF( irange.EQ.indrng ) THEN
479 
480 * RANGE='I': Compute an interval containing eigenvalues
481 * IL through IU. The initial interval [GL,GU] from the global
482 * Gerschgorin bounds GL and GU is refined by DLAEBZ.
483  itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /
484  $ log( two ) ) + 2
485  work( n+1 ) = gl
486  work( n+2 ) = gl
487  work( n+3 ) = gu
488  work( n+4 ) = gu
489  work( n+5 ) = gl
490  work( n+6 ) = gu
491  iwork( 1 ) = -1
492  iwork( 2 ) = -1
493  iwork( 3 ) = n + 1
494  iwork( 4 ) = n + 1
495  iwork( 5 ) = il - 1
496  iwork( 6 ) = iu
497 *
498  CALL dlaebz( 3, itmax, n, 2, 2, nb, atoli, rtoli, pivmin,
499  $ d, e, e2, iwork( 5 ), work( n+1 ), work( n+5 ), iout,
500  $ iwork, w, iblock, iinfo )
501  IF( iinfo .NE. 0 ) THEN
502  info = iinfo
503  RETURN
504  END IF
505 * On exit, output intervals may not be ordered by ascending negcount
506  IF( iwork( 6 ).EQ.iu ) THEN
507  wl = work( n+1 )
508  wlu = work( n+3 )
509  nwl = iwork( 1 )
510  wu = work( n+4 )
511  wul = work( n+2 )
512  nwu = iwork( 4 )
513  ELSE
514  wl = work( n+2 )
515  wlu = work( n+4 )
516  nwl = iwork( 2 )
517  wu = work( n+3 )
518  wul = work( n+1 )
519  nwu = iwork( 3 )
520  END IF
521 * On exit, the interval [WL, WLU] contains a value with negcount NWL,
522 * and [WUL, WU] contains a value with negcount NWU.
523  IF( nwl.LT.0 .OR. nwl.GE.n .OR. nwu.LT.1 .OR. nwu.GT.n ) THEN
524  info = 4
525  RETURN
526  END IF
527 
528  ELSEIF( irange.EQ.valrng ) THEN
529  wl = vl
530  wu = vu
531 
532  ELSEIF( irange.EQ.allrng ) THEN
533  wl = gl
534  wu = gu
535  ENDIF
536 
537 
538 
539 * Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU.
540 * NWL accumulates the number of eigenvalues .le. WL,
541 * NWU accumulates the number of eigenvalues .le. WU
542  m = 0
543  iend = 0
544  info = 0
545  nwl = 0
546  nwu = 0
547 *
548  DO 70 jblk = 1, nsplit
549  ioff = iend
550  ibegin = ioff + 1
551  iend = isplit( jblk )
552  in = iend - ioff
553 *
554  IF( in.EQ.1 ) THEN
555 * 1x1 block
556  IF( wl.GE.d( ibegin )-pivmin )
557  $ nwl = nwl + 1
558  IF( wu.GE.d( ibegin )-pivmin )
559  $ nwu = nwu + 1
560  IF( irange.EQ.allrng .OR.
561  $ ( wl.LT.d( ibegin )-pivmin
562  $ .AND. wu.GE. d( ibegin )-pivmin ) ) THEN
563  m = m + 1
564  w( m ) = d( ibegin )
565  werr(m) = zero
566 * The gap for a single block doesn't matter for the later
567 * algorithm and is assigned an arbitrary large value
568  iblock( m ) = jblk
569  indexw( m ) = 1
570  END IF
571 
572 * Disabled 2x2 case because of a failure on the following matrix
573 * RANGE = 'I', IL = IU = 4
574 * Original Tridiagonal, d = [
575 * -0.150102010615740E+00
576 * -0.849897989384260E+00
577 * -0.128208148052635E-15
578 * 0.128257718286320E-15
579 * ];
580 * e = [
581 * -0.357171383266986E+00
582 * -0.180411241501588E-15
583 * -0.175152352710251E-15
584 * ];
585 *
586 * ELSE IF( IN.EQ.2 ) THEN
587 ** 2x2 block
588 * DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 )
589 * TMP1 = HALF*(D(IBEGIN)+D(IEND))
590 * L1 = TMP1 - DISC
591 * IF( WL.GE. L1-PIVMIN )
592 * $ NWL = NWL + 1
593 * IF( WU.GE. L1-PIVMIN )
594 * $ NWU = NWU + 1
595 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE.
596 * $ L1-PIVMIN ) ) THEN
597 * M = M + 1
598 * W( M ) = L1
599 ** The uncertainty of eigenvalues of a 2x2 matrix is very small
600 * WERR( M ) = EPS * ABS( W( M ) ) * TWO
601 * IBLOCK( M ) = JBLK
602 * INDEXW( M ) = 1
603 * ENDIF
604 * L2 = TMP1 + DISC
605 * IF( WL.GE. L2-PIVMIN )
606 * $ NWL = NWL + 1
607 * IF( WU.GE. L2-PIVMIN )
608 * $ NWU = NWU + 1
609 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE.
610 * $ L2-PIVMIN ) ) THEN
611 * M = M + 1
612 * W( M ) = L2
613 ** The uncertainty of eigenvalues of a 2x2 matrix is very small
614 * WERR( M ) = EPS * ABS( W( M ) ) * TWO
615 * IBLOCK( M ) = JBLK
616 * INDEXW( M ) = 2
617 * ENDIF
618  ELSE
619 * General Case - block of size IN >= 2
620 * Compute local Gerschgorin interval and use it as the initial
621 * interval for DLAEBZ
622  gu = d( ibegin )
623  gl = d( ibegin )
624  tmp1 = zero
625 
626  DO 40 j = ibegin, iend
627  gl = min( gl, gers( 2*j - 1))
628  gu = max( gu, gers(2*j) )
629  40 CONTINUE
630 * [JAN/28/2009]
631 * change SPDIAM by TNORM in lines 2 and 3 thereafter
632 * line 1: remove computation of SPDIAM (not useful anymore)
633 * SPDIAM = GU - GL
634 * GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN
635 * GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN
636  gl = gl - fudge*tnorm*eps*in - fudge*pivmin
637  gu = gu + fudge*tnorm*eps*in + fudge*pivmin
638 *
639  IF( irange.GT.1 ) THEN
640  IF( gu.LT.wl ) THEN
641 * the local block contains none of the wanted eigenvalues
642  nwl = nwl + in
643  nwu = nwu + in
644  GO TO 70
645  END IF
646 * refine search interval if possible, only range (WL,WU] matters
647  gl = max( gl, wl )
648  gu = min( gu, wu )
649  IF( gl.GE.gu )
650  $ GO TO 70
651  END IF
652 
653 * Find negcount of initial interval boundaries GL and GU
654  work( n+1 ) = gl
655  work( n+in+1 ) = gu
656  CALL dlaebz( 1, 0, in, in, 1, nb, atoli, rtoli, pivmin,
657  $ d( ibegin ), e( ibegin ), e2( ibegin ),
658  $ idumma, work( n+1 ), work( n+2*in+1 ), im,
659  $ iwork, w( m+1 ), iblock( m+1 ), iinfo )
660  IF( iinfo .NE. 0 ) THEN
661  info = iinfo
662  RETURN
663  END IF
664 *
665  nwl = nwl + iwork( 1 )
666  nwu = nwu + iwork( in+1 )
667  iwoff = m - iwork( 1 )
668 
669 * Compute Eigenvalues
670  itmax = int( ( log( gu-gl+pivmin )-log( pivmin ) ) /
671  $ log( two ) ) + 2
672  CALL dlaebz( 2, itmax, in, in, 1, nb, atoli, rtoli, pivmin,
673  $ d( ibegin ), e( ibegin ), e2( ibegin ),
674  $ idumma, work( n+1 ), work( n+2*in+1 ), iout,
675  $ iwork, w( m+1 ), iblock( m+1 ), iinfo )
676  IF( iinfo .NE. 0 ) THEN
677  info = iinfo
678  RETURN
679  END IF
680 *
681 * Copy eigenvalues into W and IBLOCK
682 * Use -JBLK for block number for unconverged eigenvalues.
683 * Loop over the number of output intervals from DLAEBZ
684  DO 60 j = 1, iout
685 * eigenvalue approximation is middle point of interval
686  tmp1 = half*( work( j+n )+work( j+in+n ) )
687 * semi length of error interval
688  tmp2 = half*abs( work( j+n )-work( j+in+n ) )
689  IF( j.GT.iout-iinfo ) THEN
690 * Flag non-convergence.
691  ncnvrg = .true.
692  ib = -jblk
693  ELSE
694  ib = jblk
695  END IF
696  DO 50 je = iwork( j ) + 1 + iwoff,
697  $ iwork( j+in ) + iwoff
698  w( je ) = tmp1
699  werr( je ) = tmp2
700  indexw( je ) = je - iwoff
701  iblock( je ) = ib
702  50 CONTINUE
703  60 CONTINUE
704 *
705  m = m + im
706  END IF
707  70 CONTINUE
708 
709 * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
710 * If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
711  IF( irange.EQ.indrng ) THEN
712  idiscl = il - 1 - nwl
713  idiscu = nwu - iu
714 *
715  IF( idiscl.GT.0 ) THEN
716  im = 0
717  DO 80 je = 1, m
718 * Remove some of the smallest eigenvalues from the left so that
719 * at the end IDISCL =0. Move all eigenvalues up to the left.
720  IF( w( je ).LE.wlu .AND. idiscl.GT.0 ) THEN
721  idiscl = idiscl - 1
722  ELSE
723  im = im + 1
724  w( im ) = w( je )
725  werr( im ) = werr( je )
726  indexw( im ) = indexw( je )
727  iblock( im ) = iblock( je )
728  END IF
729  80 CONTINUE
730  m = im
731  END IF
732  IF( idiscu.GT.0 ) THEN
733 * Remove some of the largest eigenvalues from the right so that
734 * at the end IDISCU =0. Move all eigenvalues up to the left.
735  im=m+1
736  DO 81 je = m, 1, -1
737  IF( w( je ).GE.wul .AND. idiscu.GT.0 ) THEN
738  idiscu = idiscu - 1
739  ELSE
740  im = im - 1
741  w( im ) = w( je )
742  werr( im ) = werr( je )
743  indexw( im ) = indexw( je )
744  iblock( im ) = iblock( je )
745  END IF
746  81 CONTINUE
747  jee = 0
748  DO 82 je = im, m
749  jee = jee + 1
750  w( jee ) = w( je )
751  werr( jee ) = werr( je )
752  indexw( jee ) = indexw( je )
753  iblock( jee ) = iblock( je )
754  82 CONTINUE
755  m = m-im+1
756  END IF
757 
758  IF( idiscl.GT.0 .OR. idiscu.GT.0 ) THEN
759 * Code to deal with effects of bad arithmetic. (If N(w) is
760 * monotone non-decreasing, this should never happen.)
761 * Some low eigenvalues to be discarded are not in (WL,WLU],
762 * or high eigenvalues to be discarded are not in (WUL,WU]
763 * so just kill off the smallest IDISCL/largest IDISCU
764 * eigenvalues, by marking the corresponding IBLOCK = 0
765  IF( idiscl.GT.0 ) THEN
766  wkill = wu
767  DO 100 jdisc = 1, idiscl
768  iw = 0
769  DO 90 je = 1, m
770  IF( iblock( je ).NE.0 .AND.
771  $ ( w( je ).LT.wkill .OR. iw.EQ.0 ) ) THEN
772  iw = je
773  wkill = w( je )
774  END IF
775  90 CONTINUE
776  iblock( iw ) = 0
777  100 CONTINUE
778  END IF
779  IF( idiscu.GT.0 ) THEN
780  wkill = wl
781  DO 120 jdisc = 1, idiscu
782  iw = 0
783  DO 110 je = 1, m
784  IF( iblock( je ).NE.0 .AND.
785  $ ( w( je ).GE.wkill .OR. iw.EQ.0 ) ) THEN
786  iw = je
787  wkill = w( je )
788  END IF
789  110 CONTINUE
790  iblock( iw ) = 0
791  120 CONTINUE
792  END IF
793 * Now erase all eigenvalues with IBLOCK set to zero
794  im = 0
795  DO 130 je = 1, m
796  IF( iblock( je ).NE.0 ) THEN
797  im = im + 1
798  w( im ) = w( je )
799  werr( im ) = werr( je )
800  indexw( im ) = indexw( je )
801  iblock( im ) = iblock( je )
802  END IF
803  130 CONTINUE
804  m = im
805  END IF
806  IF( idiscl.LT.0 .OR. idiscu.LT.0 ) THEN
807  toofew = .true.
808  END IF
809  END IF
810 *
811  IF(( irange.EQ.allrng .AND. m.NE.n ).OR.
812  $ ( irange.EQ.indrng .AND. m.NE.iu-il+1 ) ) THEN
813  toofew = .true.
814  END IF
815 
816 * If ORDER='B', do nothing the eigenvalues are already sorted by
817 * block.
818 * If ORDER='E', sort the eigenvalues from smallest to largest
819 
820  IF( lsame(order,'E') .AND. nsplit.GT.1 ) THEN
821  DO 150 je = 1, m - 1
822  ie = 0
823  tmp1 = w( je )
824  DO 140 j = je + 1, m
825  IF( w( j ).LT.tmp1 ) THEN
826  ie = j
827  tmp1 = w( j )
828  END IF
829  140 CONTINUE
830  IF( ie.NE.0 ) THEN
831  tmp2 = werr( ie )
832  itmp1 = iblock( ie )
833  itmp2 = indexw( ie )
834  w( ie ) = w( je )
835  werr( ie ) = werr( je )
836  iblock( ie ) = iblock( je )
837  indexw( ie ) = indexw( je )
838  w( je ) = tmp1
839  werr( je ) = tmp2
840  iblock( je ) = itmp1
841  indexw( je ) = itmp2
842  END IF
843  150 CONTINUE
844  END IF
845 *
846  info = 0
847  IF( ncnvrg )
848  $ info = info + 1
849  IF( toofew )
850  $ info = info + 2
851  RETURN
852 *
853 * End of DLARRD
854 *
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dlaebz(IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, NAB, WORK, IWORK, INFO)
DLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than ...
Definition: dlaebz.f:321
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function: