LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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
          If RANGE='V', the lower bound 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]VU
          VU is DOUBLE PRECISION
          If RANGE='V', the upper bound 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
          If RANGE='I', the index of the
          smallest eigenvalue 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]IU
          IU is INTEGER
          If RANGE='I', the index of the
          largest eigenvalue 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
June 2016

Definition at line 331 of file dlarrd.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: