LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

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