 LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 subroutine dlarrd ( character RANGE, character ORDER, integer N, double precision VL, double precision VU, integer IL, integer IU, double precision, dimension( * ) GERS, double precision RELTOL, double precision, dimension( * ) D, double precision, dimension( * ) E, double precision, dimension( * ) E2, double precision PIVMIN, integer NSPLIT, integer, dimension( * ) ISPLIT, integer M, double precision, dimension( * ) W, double precision, dimension( * ) WERR, double precision WL, double precision WU, integer, dimension( * ) IBLOCK, integer, dimension( * ) INDEXW, double precision, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO )

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

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

Definition at line 331 of file dlarrd.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: