 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ dlarrd()

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

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

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

Definition at line 325 of file dlarrd.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  DOUBLE PRECISION PIVMIN, RELTOL, VL, VU, WL, WU
338 * ..
339 * .. Array Arguments ..
340  INTEGER IBLOCK( * ), INDEXW( * ),
341  \$ ISPLIT( * ), IWORK( * )
342  DOUBLE PRECISION D( * ), E( * ), E2( * ),
343  \$ GERS( * ), W( * ), WERR( * ), WORK( * )
344 * ..
345 *
346 * =====================================================================
347 *
348 * .. Parameters ..
349  DOUBLE PRECISION ZERO, ONE, TWO, HALF, FUDGE
350  parameter( zero = 0.0d0, one = 1.0d0,
351  \$ two = 2.0d0, 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  DOUBLE PRECISION 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  DOUBLE PRECISION DLAMCH
373  EXTERNAL lsame, ilaenv, dlamch
374 * ..
375 * .. External Subroutines ..
376  EXTERNAL dlaebz
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 = dlamch( 'P' )
440  uflow = dlamch( '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, 'DSTEBZ', ' ', 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 DLAEBZ:
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 DLAEBZ.
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 dlaebz( 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 DLAEBZ
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 dlaebz( 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 dlaebz( 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 DLAEBZ
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 DLARRD
865 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
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: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
Here is the call graph for this function:
Here is the caller graph for this function: