LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dlarrd()

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