LAPACK  3.8.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.

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
Date
June 2016

Definition at line 331 of file slarrd.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  REAL pivmin, reltol, vl, vu, wl, wu
341 * ..
342 * .. Array Arguments ..
343  INTEGER iblock( * ), indexw( * ),
344  \$ isplit( * ), iwork( * )
345  REAL d( * ), e( * ), e2( * ),
346  \$ gers( * ), w( * ), werr( * ), work( * )
347 * ..
348 *
349 * =====================================================================
350 *
351 * .. Parameters ..
352  REAL zero, one, two, half, fudge
353  parameter( zero = 0.0e0, one = 1.0e0,
354  \$ two = 2.0e0, 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  REAL 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  REAL slamch
376  EXTERNAL lsame, ilaenv, slamch
377 * ..
378 * .. External Subroutines ..
379  EXTERNAL slaebz
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 = slamch( 'P' )
443  uflow = slamch( '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, 'SSTEBZ', ' ', 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 SLAEBZ:
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 SLAEBZ.
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 slaebz( 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 SLAEBZ
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 slaebz( 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 slaebz( 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 SLAEBZ
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 SLARRD
868 *
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
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:321
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
Here is the call graph for this function:
Here is the caller graph for this function: