LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slarrd()

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

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

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

Purpose:
 SLARRD computes the eigenvalues of a symmetric tridiagonal
 matrix T to suitable accuracy. This is an auxiliary code to be
 called from SSTEMR.
 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 REAL
          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 REAL
          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 REAL 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 REAL
          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 REAL array, dimension (N)
          The n diagonal elements of the tridiagonal matrix T.
[in]E
          E is REAL array, dimension (N-1)
          The (n-1) off-diagonal elements of the tridiagonal matrix T.
[in]E2
          E2 is REAL array, dimension (N-1)
          The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
[in]PIVMIN
          PIVMIN is REAL
          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 REAL array, dimension (N)
          On exit, the first M elements of W will contain the
          eigenvalue approximations. SLARRD 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 REAL array, dimension (N)
          The error bound on the corresponding eigenvalue approximation
          in W.
[out]WL
          WL is REAL
[out]WU
          WU is REAL
          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 SLAEBZ 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.  (SLARRD 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 REAL 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   REAL, 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.

Definition at line 325 of file slarrd.f.

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