LAPACK 3.3.0

dlarrd.f

Go to the documentation of this file.
00001       SUBROUTINE DLARRD( RANGE, ORDER, N, VL, VU, IL, IU, GERS,
00002      $                    RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
00003      $                    M, W, WERR, WL, WU, IBLOCK, INDEXW,
00004      $                    WORK, IWORK, INFO )
00005 *
00006 *  -- LAPACK auxiliary routine (version 3.3.0)                        --
00007 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00008 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00009 *     November 2010
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          ORDER, RANGE
00013       INTEGER            IL, INFO, IU, M, N, NSPLIT
00014       DOUBLE PRECISION    PIVMIN, RELTOL, VL, VU, WL, WU
00015 *     ..
00016 *     .. Array Arguments ..
00017       INTEGER            IBLOCK( * ), INDEXW( * ),
00018      $                   ISPLIT( * ), IWORK( * )
00019       DOUBLE PRECISION   D( * ), E( * ), E2( * ),
00020      $                   GERS( * ), W( * ), WERR( * ), WORK( * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  DLARRD computes the eigenvalues of a symmetric tridiagonal
00027 *  matrix T to suitable accuracy. This is an auxiliary code to be
00028 *  called from DSTEMR.
00029 *  The user may ask for all eigenvalues, all eigenvalues
00030 *  in the half-open interval (VL, VU], or the IL-th through IU-th
00031 *  eigenvalues.
00032 *
00033 *  To avoid overflow, the matrix must be scaled so that its
00034 *  largest element is no greater than overflow**(1/2) *
00035 *  underflow**(1/4) in absolute value, and for greatest
00036 *  accuracy, it should not be much smaller than that.
00037 *
00038 *  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
00039 *  Matrix", Report CS41, Computer Science Dept., Stanford
00040 *  University, July 21, 1966.
00041 *
00042 *  Arguments
00043 *  =========
00044 *
00045 *  RANGE   (input) CHARACTER*1
00046 *          = 'A': ("All")   all eigenvalues will be found.
00047 *          = 'V': ("Value") all eigenvalues in the half-open interval
00048 *                           (VL, VU] will be found.
00049 *          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
00050 *                           entire matrix) will be found.
00051 *
00052 *  ORDER   (input) CHARACTER*1
00053 *          = 'B': ("By Block") the eigenvalues will be grouped by
00054 *                              split-off block (see IBLOCK, ISPLIT) and
00055 *                              ordered from smallest to largest within
00056 *                              the block.
00057 *          = 'E': ("Entire matrix")
00058 *                              the eigenvalues for the entire matrix
00059 *                              will be ordered from smallest to
00060 *                              largest.
00061 *
00062 *  N       (input) INTEGER
00063 *          The order of the tridiagonal matrix T.  N >= 0.
00064 *
00065 *  VL      (input) DOUBLE PRECISION
00066 *  VU      (input) DOUBLE PRECISION
00067 *          If RANGE='V', the lower and upper bounds of the interval to
00068 *          be searched for eigenvalues.  Eigenvalues less than or equal
00069 *          to VL, or greater than VU, will not be returned.  VL < VU.
00070 *          Not referenced if RANGE = 'A' or 'I'.
00071 *
00072 *  IL      (input) INTEGER
00073 *  IU      (input) INTEGER
00074 *          If RANGE='I', the indices (in ascending order) of the
00075 *          smallest and largest eigenvalues to be returned.
00076 *          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
00077 *          Not referenced if RANGE = 'A' or 'V'.
00078 *
00079 *  GERS    (input) DOUBLE PRECISION array, dimension (2*N)
00080 *          The N Gerschgorin intervals (the i-th Gerschgorin interval
00081 *          is (GERS(2*i-1), GERS(2*i)).
00082 *
00083 *  RELTOL  (input) DOUBLE PRECISION
00084 *          The minimum relative width of an interval.  When an interval
00085 *          is narrower than RELTOL times the larger (in
00086 *          magnitude) endpoint, then it is considered to be
00087 *          sufficiently small, i.e., converged.  Note: this should
00088 *          always be at least radix*machine epsilon.
00089 *
00090 *  D       (input) DOUBLE PRECISION array, dimension (N)
00091 *          The n diagonal elements of the tridiagonal matrix T.
00092 *
00093 *  E       (input) DOUBLE PRECISION array, dimension (N-1)
00094 *          The (n-1) off-diagonal elements of the tridiagonal matrix T.
00095 *
00096 *  E2      (input) DOUBLE PRECISION array, dimension (N-1)
00097 *          The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
00098 *
00099 *  PIVMIN  (input) DOUBLE PRECISION
00100 *          The minimum pivot allowed in the Sturm sequence for T.
00101 *
00102 *  NSPLIT  (input) INTEGER
00103 *          The number of diagonal blocks in the matrix T.
00104 *          1 <= NSPLIT <= N.
00105 *
00106 *  ISPLIT  (input) INTEGER array, dimension (N)
00107 *          The splitting points, at which T breaks up into submatrices.
00108 *          The first submatrix consists of rows/columns 1 to ISPLIT(1),
00109 *          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
00110 *          etc., and the NSPLIT-th consists of rows/columns
00111 *          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
00112 *          (Only the first NSPLIT elements will actually be used, but
00113 *          since the user cannot know a priori what value NSPLIT will
00114 *          have, N words must be reserved for ISPLIT.)
00115 *
00116 *  M       (output) INTEGER
00117 *          The actual number of eigenvalues found. 0 <= M <= N.
00118 *          (See also the description of INFO=2,3.)
00119 *
00120 *  W       (output) DOUBLE PRECISION array, dimension (N)
00121 *          On exit, the first M elements of W will contain the
00122 *          eigenvalue approximations. DLARRD computes an interval
00123 *          I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue
00124 *          approximation is given as the interval midpoint
00125 *          W(j)= ( a_j + b_j)/2. The corresponding error is bounded by
00126 *          WERR(j) = abs( a_j - b_j)/2
00127 *
00128 *  WERR    (output) DOUBLE PRECISION array, dimension (N)
00129 *          The error bound on the corresponding eigenvalue approximation
00130 *          in W.
00131 *
00132 *  WL      (output) DOUBLE PRECISION
00133 *  WU      (output) DOUBLE PRECISION
00134 *          The interval (WL, WU] contains all the wanted eigenvalues.
00135 *          If RANGE='V', then WL=VL and WU=VU.
00136 *          If RANGE='A', then WL and WU are the global Gerschgorin bounds
00137 *                        on the spectrum.
00138 *          If RANGE='I', then WL and WU are computed by DLAEBZ from the
00139 *                        index range specified.
00140 *
00141 *  IBLOCK  (output) INTEGER array, dimension (N)
00142 *          At each row/column j where E(j) is zero or small, the
00143 *          matrix T is considered to split into a block diagonal
00144 *          matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which
00145 *          block (from 1 to the number of blocks) the eigenvalue W(i)
00146 *          belongs.  (DLARRD may use the remaining N-M elements as
00147 *          workspace.)
00148 *
00149 *  INDEXW  (output) INTEGER array, dimension (N)
00150 *          The indices of the eigenvalues within each block (submatrix);
00151 *          for example, INDEXW(i)= j and IBLOCK(i)=k imply that the
00152 *          i-th eigenvalue W(i) is the j-th eigenvalue in block k.
00153 *
00154 *  WORK    (workspace) DOUBLE PRECISION array, dimension (4*N)
00155 *
00156 *  IWORK   (workspace) INTEGER array, dimension (3*N)
00157 *
00158 *  INFO    (output) INTEGER
00159 *          = 0:  successful exit
00160 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00161 *          > 0:  some or all of the eigenvalues failed to converge or
00162 *                were not computed:
00163 *                =1 or 3: Bisection failed to converge for some
00164 *                        eigenvalues; these eigenvalues are flagged by a
00165 *                        negative block number.  The effect is that the
00166 *                        eigenvalues may not be as accurate as the
00167 *                        absolute and relative tolerances.  This is
00168 *                        generally caused by unexpectedly inaccurate
00169 *                        arithmetic.
00170 *                =2 or 3: RANGE='I' only: Not all of the eigenvalues
00171 *                        IL:IU were found.
00172 *                        Effect: M < IU+1-IL
00173 *                        Cause:  non-monotonic arithmetic, causing the
00174 *                                Sturm sequence to be non-monotonic.
00175 *                        Cure:   recalculate, using RANGE='A', and pick
00176 *                                out eigenvalues IL:IU.  In some cases,
00177 *                                increasing the PARAMETER "FUDGE" may
00178 *                                make things work.
00179 *                = 4:    RANGE='I', and the Gershgorin interval
00180 *                        initially used was too small.  No eigenvalues
00181 *                        were computed.
00182 *                        Probable cause: your machine has sloppy
00183 *                                        floating-point arithmetic.
00184 *                        Cure: Increase the PARAMETER "FUDGE",
00185 *                              recompile, and try again.
00186 *
00187 *  Internal Parameters
00188 *  ===================
00189 *
00190 *  FUDGE   DOUBLE PRECISION, default = 2
00191 *          A "fudge factor" to widen the Gershgorin intervals.  Ideally,
00192 *          a value of 1 should work, but on machines with sloppy
00193 *          arithmetic, this needs to be larger.  The default for
00194 *          publicly released versions should be large enough to handle
00195 *          the worst machine around.  Note that this has no effect
00196 *          on accuracy of the solution.
00197 *
00198 *  Based on contributions by
00199 *     W. Kahan, University of California, Berkeley, USA
00200 *     Beresford Parlett, University of California, Berkeley, USA
00201 *     Jim Demmel, University of California, Berkeley, USA
00202 *     Inderjit Dhillon, University of Texas, Austin, USA
00203 *     Osni Marques, LBNL/NERSC, USA
00204 *     Christof Voemel, University of California, Berkeley, USA
00205 *
00206 *  =====================================================================
00207 *
00208 *     .. Parameters ..
00209       DOUBLE PRECISION   ZERO, ONE, TWO, HALF, FUDGE
00210       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0,
00211      $                     TWO = 2.0D0, HALF = ONE/TWO,
00212      $                     FUDGE = TWO )
00213       INTEGER   ALLRNG, VALRNG, INDRNG
00214       PARAMETER ( ALLRNG = 1, VALRNG = 2, INDRNG = 3 )
00215 *     ..
00216 *     .. Local Scalars ..
00217       LOGICAL            NCNVRG, TOOFEW
00218       INTEGER            I, IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
00219      $                   IM, IN, IOFF, IOUT, IRANGE, ITMAX, ITMP1,
00220      $                   ITMP2, IW, IWOFF, J, JBLK, JDISC, JE, JEE, NB,
00221      $                   NWL, NWU
00222       DOUBLE PRECISION   ATOLI, EPS, GL, GU, RTOLI, TMP1, TMP2,
00223      $                   TNORM, UFLOW, WKILL, WLU, WUL
00224 
00225 *     ..
00226 *     .. Local Arrays ..
00227       INTEGER            IDUMMA( 1 )
00228 *     ..
00229 *     .. External Functions ..
00230       LOGICAL            LSAME
00231       INTEGER            ILAENV
00232       DOUBLE PRECISION   DLAMCH
00233       EXTERNAL           LSAME, ILAENV, DLAMCH
00234 *     ..
00235 *     .. External Subroutines ..
00236       EXTERNAL           DLAEBZ
00237 *     ..
00238 *     .. Intrinsic Functions ..
00239       INTRINSIC          ABS, INT, LOG, MAX, MIN
00240 *     ..
00241 *     .. Executable Statements ..
00242 *
00243       INFO = 0
00244 *
00245 *     Decode RANGE
00246 *
00247       IF( LSAME( RANGE, 'A' ) ) THEN
00248          IRANGE = ALLRNG
00249       ELSE IF( LSAME( RANGE, 'V' ) ) THEN
00250          IRANGE = VALRNG
00251       ELSE IF( LSAME( RANGE, 'I' ) ) THEN
00252          IRANGE = INDRNG
00253       ELSE
00254          IRANGE = 0
00255       END IF
00256 *
00257 *     Check for Errors
00258 *
00259       IF( IRANGE.LE.0 ) THEN
00260          INFO = -1
00261       ELSE IF( .NOT.(LSAME(ORDER,'B').OR.LSAME(ORDER,'E')) ) THEN
00262          INFO = -2
00263       ELSE IF( N.LT.0 ) THEN
00264          INFO = -3
00265       ELSE IF( IRANGE.EQ.VALRNG ) THEN
00266          IF( VL.GE.VU )
00267      $      INFO = -5
00268       ELSE IF( IRANGE.EQ.INDRNG .AND.
00269      $        ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) THEN
00270          INFO = -6
00271       ELSE IF( IRANGE.EQ.INDRNG .AND.
00272      $        ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN
00273          INFO = -7
00274       END IF
00275 *
00276       IF( INFO.NE.0 ) THEN
00277          RETURN
00278       END IF
00279 
00280 *     Initialize error flags
00281       INFO = 0
00282       NCNVRG = .FALSE.
00283       TOOFEW = .FALSE.
00284 
00285 *     Quick return if possible
00286       M = 0
00287       IF( N.EQ.0 ) RETURN
00288 
00289 *     Simplification:
00290       IF( IRANGE.EQ.INDRNG .AND. IL.EQ.1 .AND. IU.EQ.N ) IRANGE = 1
00291 
00292 *     Get machine constants
00293       EPS = DLAMCH( 'P' )
00294       UFLOW = DLAMCH( 'U' )
00295 
00296 
00297 *     Special Case when N=1
00298 *     Treat case of 1x1 matrix for quick return
00299       IF( N.EQ.1 ) THEN
00300          IF( (IRANGE.EQ.ALLRNG).OR.
00301      $       ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
00302      $       ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
00303             M = 1
00304             W(1) = D(1)
00305 *           The computation error of the eigenvalue is zero
00306             WERR(1) = ZERO
00307             IBLOCK( 1 ) = 1
00308             INDEXW( 1 ) = 1
00309          ENDIF
00310          RETURN
00311       END IF
00312 
00313 *     NB is the minimum vector length for vector bisection, or 0
00314 *     if only scalar is to be done.
00315       NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 )
00316       IF( NB.LE.1 ) NB = 0
00317 
00318 *     Find global spectral radius
00319       GL = D(1)
00320       GU = D(1)
00321       DO 5 I = 1,N
00322          GL =  MIN( GL, GERS( 2*I - 1))
00323          GU = MAX( GU, GERS(2*I) )
00324  5    CONTINUE
00325 *     Compute global Gerschgorin bounds and spectral diameter
00326       TNORM = MAX( ABS( GL ), ABS( GU ) )
00327       GL = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
00328       GU = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
00329 *     [JAN/28/2009] remove the line below since SPDIAM variable not use
00330 *     SPDIAM = GU - GL
00331 *     Input arguments for DLAEBZ:
00332 *     The relative tolerance.  An interval (a,b] lies within
00333 *     "relative tolerance" if  b-a < RELTOL*max(|a|,|b|),
00334       RTOLI = RELTOL
00335 *     Set the absolute tolerance for interval convergence to zero to force
00336 *     interval convergence based on relative size of the interval.
00337 *     This is dangerous because intervals might not converge when RELTOL is
00338 *     small. But at least a very small number should be selected so that for
00339 *     strongly graded matrices, the code can get relatively accurate
00340 *     eigenvalues.
00341       ATOLI = FUDGE*TWO*UFLOW + FUDGE*TWO*PIVMIN
00342 
00343       IF( IRANGE.EQ.INDRNG ) THEN
00344 
00345 *        RANGE='I': Compute an interval containing eigenvalues
00346 *        IL through IU. The initial interval [GL,GU] from the global
00347 *        Gerschgorin bounds GL and GU is refined by DLAEBZ.
00348          ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
00349      $           LOG( TWO ) ) + 2
00350          WORK( N+1 ) = GL
00351          WORK( N+2 ) = GL
00352          WORK( N+3 ) = GU
00353          WORK( N+4 ) = GU
00354          WORK( N+5 ) = GL
00355          WORK( N+6 ) = GU
00356          IWORK( 1 ) = -1
00357          IWORK( 2 ) = -1
00358          IWORK( 3 ) = N + 1
00359          IWORK( 4 ) = N + 1
00360          IWORK( 5 ) = IL - 1
00361          IWORK( 6 ) = IU
00362 *
00363          CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN,
00364      $         D, E, E2, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
00365      $                IWORK, W, IBLOCK, IINFO )
00366          IF( IINFO .NE. 0 ) THEN
00367             INFO = IINFO
00368             RETURN
00369          END IF
00370 *        On exit, output intervals may not be ordered by ascending negcount
00371          IF( IWORK( 6 ).EQ.IU ) THEN
00372             WL = WORK( N+1 )
00373             WLU = WORK( N+3 )
00374             NWL = IWORK( 1 )
00375             WU = WORK( N+4 )
00376             WUL = WORK( N+2 )
00377             NWU = IWORK( 4 )
00378          ELSE
00379             WL = WORK( N+2 )
00380             WLU = WORK( N+4 )
00381             NWL = IWORK( 2 )
00382             WU = WORK( N+3 )
00383             WUL = WORK( N+1 )
00384             NWU = IWORK( 3 )
00385          END IF
00386 *        On exit, the interval [WL, WLU] contains a value with negcount NWL,
00387 *        and [WUL, WU] contains a value with negcount NWU.
00388          IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN
00389             INFO = 4
00390             RETURN
00391          END IF
00392 
00393       ELSEIF( IRANGE.EQ.VALRNG ) THEN
00394          WL = VL
00395          WU = VU
00396 
00397       ELSEIF( IRANGE.EQ.ALLRNG ) THEN
00398          WL = GL
00399          WU = GU
00400       ENDIF
00401 
00402 
00403 
00404 *     Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU.
00405 *     NWL accumulates the number of eigenvalues .le. WL,
00406 *     NWU accumulates the number of eigenvalues .le. WU
00407       M = 0
00408       IEND = 0
00409       INFO = 0
00410       NWL = 0
00411       NWU = 0
00412 *
00413       DO 70 JBLK = 1, NSPLIT
00414          IOFF = IEND
00415          IBEGIN = IOFF + 1
00416          IEND = ISPLIT( JBLK )
00417          IN = IEND - IOFF
00418 *
00419          IF( IN.EQ.1 ) THEN
00420 *           1x1 block
00421             IF( WL.GE.D( IBEGIN )-PIVMIN )
00422      $         NWL = NWL + 1
00423             IF( WU.GE.D( IBEGIN )-PIVMIN )
00424      $         NWU = NWU + 1
00425             IF( IRANGE.EQ.ALLRNG .OR.
00426      $           ( WL.LT.D( IBEGIN )-PIVMIN
00427      $             .AND. WU.GE. D( IBEGIN )-PIVMIN ) ) THEN
00428                M = M + 1
00429                W( M ) = D( IBEGIN )
00430                WERR(M) = ZERO
00431 *              The gap for a single block doesn't matter for the later
00432 *              algorithm and is assigned an arbitrary large value
00433                IBLOCK( M ) = JBLK
00434                INDEXW( M ) = 1
00435             END IF
00436 
00437 *        Disabled 2x2 case because of a failure on the following matrix
00438 *        RANGE = 'I', IL = IU = 4
00439 *          Original Tridiagonal, d = [
00440 *           -0.150102010615740E+00
00441 *           -0.849897989384260E+00
00442 *           -0.128208148052635E-15
00443 *            0.128257718286320E-15
00444 *          ];
00445 *          e = [
00446 *           -0.357171383266986E+00
00447 *           -0.180411241501588E-15
00448 *           -0.175152352710251E-15
00449 *          ];
00450 *
00451 *         ELSE IF( IN.EQ.2 ) THEN
00452 **           2x2 block
00453 *            DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 )
00454 *            TMP1 = HALF*(D(IBEGIN)+D(IEND))
00455 *            L1 = TMP1 - DISC
00456 *            IF( WL.GE. L1-PIVMIN )
00457 *     $         NWL = NWL + 1
00458 *            IF( WU.GE. L1-PIVMIN )
00459 *     $         NWU = NWU + 1
00460 *            IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE.
00461 *     $          L1-PIVMIN ) ) THEN
00462 *               M = M + 1
00463 *               W( M ) = L1
00464 **              The uncertainty of eigenvalues of a 2x2 matrix is very small
00465 *               WERR( M ) = EPS * ABS( W( M ) ) * TWO
00466 *               IBLOCK( M ) = JBLK
00467 *               INDEXW( M ) = 1
00468 *            ENDIF
00469 *            L2 = TMP1 + DISC
00470 *            IF( WL.GE. L2-PIVMIN )
00471 *     $         NWL = NWL + 1
00472 *            IF( WU.GE. L2-PIVMIN )
00473 *     $         NWU = NWU + 1
00474 *            IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE.
00475 *     $          L2-PIVMIN ) ) THEN
00476 *               M = M + 1
00477 *               W( M ) = L2
00478 **              The uncertainty of eigenvalues of a 2x2 matrix is very small
00479 *               WERR( M ) = EPS * ABS( W( M ) ) * TWO
00480 *               IBLOCK( M ) = JBLK
00481 *               INDEXW( M ) = 2
00482 *            ENDIF
00483          ELSE
00484 *           General Case - block of size IN >= 2
00485 *           Compute local Gerschgorin interval and use it as the initial
00486 *           interval for DLAEBZ
00487             GU = D( IBEGIN )
00488             GL = D( IBEGIN )
00489             TMP1 = ZERO
00490 
00491             DO 40 J = IBEGIN, IEND
00492                GL =  MIN( GL, GERS( 2*J - 1))
00493                GU = MAX( GU, GERS(2*J) )
00494    40       CONTINUE
00495 *           [JAN/28/2009]
00496 *           change SPDIAM by TNORM in lines 2 and 3 thereafter
00497 *           line 1: remove computation of SPDIAM (not useful anymore)
00498 *           SPDIAM = GU - GL
00499 *           GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN
00500 *           GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN
00501             GL = GL - FUDGE*TNORM*EPS*IN - FUDGE*PIVMIN
00502             GU = GU + FUDGE*TNORM*EPS*IN + FUDGE*PIVMIN
00503 *
00504             IF( IRANGE.GT.1 ) THEN
00505                IF( GU.LT.WL ) THEN
00506 *                 the local block contains none of the wanted eigenvalues
00507                   NWL = NWL + IN
00508                   NWU = NWU + IN
00509                   GO TO 70
00510                END IF
00511 *              refine search interval if possible, only range (WL,WU] matters
00512                GL = MAX( GL, WL )
00513                GU = MIN( GU, WU )
00514                IF( GL.GE.GU )
00515      $            GO TO 70
00516             END IF
00517 
00518 *           Find negcount of initial interval boundaries GL and GU
00519             WORK( N+1 ) = GL
00520             WORK( N+IN+1 ) = GU
00521             CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
00522      $                   D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ),
00523      $                   IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
00524      $                   IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
00525             IF( IINFO .NE. 0 ) THEN
00526                INFO = IINFO
00527                RETURN
00528             END IF
00529 *
00530             NWL = NWL + IWORK( 1 )
00531             NWU = NWU + IWORK( IN+1 )
00532             IWOFF = M - IWORK( 1 )
00533 
00534 *           Compute Eigenvalues
00535             ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) /
00536      $              LOG( TWO ) ) + 2
00537             CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
00538      $                   D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ),
00539      $                   IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
00540      $                   IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
00541             IF( IINFO .NE. 0 ) THEN
00542                INFO = IINFO
00543                RETURN
00544             END IF
00545 *
00546 *           Copy eigenvalues into W and IBLOCK
00547 *           Use -JBLK for block number for unconverged eigenvalues.
00548 *           Loop over the number of output intervals from DLAEBZ
00549             DO 60 J = 1, IOUT
00550 *              eigenvalue approximation is middle point of interval
00551                TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) )
00552 *              semi length of error interval
00553                TMP2 = HALF*ABS( WORK( J+N )-WORK( J+IN+N ) )
00554                IF( J.GT.IOUT-IINFO ) THEN
00555 *                 Flag non-convergence.
00556                   NCNVRG = .TRUE.
00557                   IB = -JBLK
00558                ELSE
00559                   IB = JBLK
00560                END IF
00561                DO 50 JE = IWORK( J ) + 1 + IWOFF,
00562      $                 IWORK( J+IN ) + IWOFF
00563                   W( JE ) = TMP1
00564                   WERR( JE ) = TMP2
00565                   INDEXW( JE ) = JE - IWOFF
00566                   IBLOCK( JE ) = IB
00567    50          CONTINUE
00568    60       CONTINUE
00569 *
00570             M = M + IM
00571          END IF
00572    70 CONTINUE
00573 
00574 *     If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
00575 *     If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
00576       IF( IRANGE.EQ.INDRNG ) THEN
00577          IDISCL = IL - 1 - NWL
00578          IDISCU = NWU - IU
00579 *
00580          IF( IDISCL.GT.0 ) THEN
00581             IM = 0
00582             DO 80 JE = 1, M
00583 *              Remove some of the smallest eigenvalues from the left so that
00584 *              at the end IDISCL =0. Move all eigenvalues up to the left.
00585                IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN
00586                   IDISCL = IDISCL - 1
00587                ELSE
00588                   IM = IM + 1
00589                   W( IM ) = W( JE )
00590                   WERR( IM ) = WERR( JE )
00591                   INDEXW( IM ) = INDEXW( JE )
00592                   IBLOCK( IM ) = IBLOCK( JE )
00593                END IF
00594  80         CONTINUE
00595             M = IM
00596          END IF
00597          IF( IDISCU.GT.0 ) THEN
00598 *           Remove some of the largest eigenvalues from the right so that
00599 *           at the end IDISCU =0. Move all eigenvalues up to the left.
00600             IM=M+1
00601             DO 81 JE = M, 1, -1
00602                IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN
00603                   IDISCU = IDISCU - 1
00604                ELSE
00605                   IM = IM - 1
00606                   W( IM ) = W( JE )
00607                   WERR( IM ) = WERR( JE )
00608                   INDEXW( IM ) = INDEXW( JE )
00609                   IBLOCK( IM ) = IBLOCK( JE )
00610                END IF
00611  81         CONTINUE
00612             JEE = 0
00613             DO 82 JE = IM, M
00614                JEE = JEE + 1
00615                W( JEE ) = W( JE )
00616                WERR( JEE ) = WERR( JE )
00617                INDEXW( JEE ) = INDEXW( JE )
00618                IBLOCK( JEE ) = IBLOCK( JE )
00619  82         CONTINUE
00620             M = M-IM+1
00621          END IF
00622 
00623          IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
00624 *           Code to deal with effects of bad arithmetic. (If N(w) is
00625 *           monotone non-decreasing, this should never happen.)
00626 *           Some low eigenvalues to be discarded are not in (WL,WLU],
00627 *           or high eigenvalues to be discarded are not in (WUL,WU]
00628 *           so just kill off the smallest IDISCL/largest IDISCU
00629 *           eigenvalues, by marking the corresponding IBLOCK = 0
00630             IF( IDISCL.GT.0 ) THEN
00631                WKILL = WU
00632                DO 100 JDISC = 1, IDISCL
00633                   IW = 0
00634                   DO 90 JE = 1, M
00635                      IF( IBLOCK( JE ).NE.0 .AND.
00636      $                    ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN
00637                         IW = JE
00638                         WKILL = W( JE )
00639                      END IF
00640  90               CONTINUE
00641                   IBLOCK( IW ) = 0
00642  100           CONTINUE
00643             END IF
00644             IF( IDISCU.GT.0 ) THEN
00645                WKILL = WL
00646                DO 120 JDISC = 1, IDISCU
00647                   IW = 0
00648                   DO 110 JE = 1, M
00649                      IF( IBLOCK( JE ).NE.0 .AND.
00650      $                    ( W( JE ).GE.WKILL .OR. IW.EQ.0 ) ) THEN
00651                         IW = JE
00652                         WKILL = W( JE )
00653                      END IF
00654  110              CONTINUE
00655                   IBLOCK( IW ) = 0
00656  120           CONTINUE
00657             END IF
00658 *           Now erase all eigenvalues with IBLOCK set to zero
00659             IM = 0
00660             DO 130 JE = 1, M
00661                IF( IBLOCK( JE ).NE.0 ) THEN
00662                   IM = IM + 1
00663                   W( IM ) = W( JE )
00664                   WERR( IM ) = WERR( JE )
00665                   INDEXW( IM ) = INDEXW( JE )
00666                   IBLOCK( IM ) = IBLOCK( JE )
00667                END IF
00668  130        CONTINUE
00669             M = IM
00670          END IF
00671          IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN
00672             TOOFEW = .TRUE.
00673          END IF
00674       END IF
00675 *
00676       IF(( IRANGE.EQ.ALLRNG .AND. M.NE.N ).OR.
00677      $   ( IRANGE.EQ.INDRNG .AND. M.NE.IU-IL+1 ) ) THEN
00678          TOOFEW = .TRUE.
00679       END IF
00680 
00681 *     If ORDER='B', do nothing the eigenvalues are already sorted by
00682 *        block.
00683 *     If ORDER='E', sort the eigenvalues from smallest to largest
00684 
00685       IF( LSAME(ORDER,'E') .AND. NSPLIT.GT.1 ) THEN
00686          DO 150 JE = 1, M - 1
00687             IE = 0
00688             TMP1 = W( JE )
00689             DO 140 J = JE + 1, M
00690                IF( W( J ).LT.TMP1 ) THEN
00691                   IE = J
00692                   TMP1 = W( J )
00693                END IF
00694   140       CONTINUE
00695             IF( IE.NE.0 ) THEN
00696                TMP2 = WERR( IE )
00697                ITMP1 = IBLOCK( IE )
00698                ITMP2 = INDEXW( IE )
00699                W( IE ) = W( JE )
00700                WERR( IE ) = WERR( JE )
00701                IBLOCK( IE ) = IBLOCK( JE )
00702                INDEXW( IE ) = INDEXW( JE )
00703                W( JE ) = TMP1
00704                WERR( JE ) = TMP2
00705                IBLOCK( JE ) = ITMP1
00706                INDEXW( JE ) = ITMP2
00707             END IF
00708   150    CONTINUE
00709       END IF
00710 *
00711       INFO = 0
00712       IF( NCNVRG )
00713      $   INFO = INFO + 1
00714       IF( TOOFEW )
00715      $   INFO = INFO + 2
00716       RETURN
00717 *
00718 *     End of DLARRD
00719 *
00720       END
 All Files Functions