LAPACK: SRC/dstebz.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
00002      $                   M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
00003      $                   INFO )
00004 *
00005 *  -- LAPACK routine (version 3.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     November 2006
00009 *     8-18-00:  Increase FUDGE factor for T3E (eca)
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          ORDER, RANGE
00013       INTEGER            IL, INFO, IU, M, N, NSPLIT
00014       DOUBLE PRECISION   ABSTOL, VL, VU
00015 *     ..
00016 *     .. Array Arguments ..
00017       INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * )
00018       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  DSTEBZ computes the eigenvalues of a symmetric tridiagonal
00025 *  matrix T.  The user may ask for all eigenvalues, all eigenvalues
00026 *  in the half-open interval (VL, VU], or the IL-th through IU-th
00027 *  eigenvalues.
00028 *
00029 *  To avoid overflow, the matrix must be scaled so that its
00030 *  largest element is no greater than overflow**(1/2) *
00031 *  underflow**(1/4) in absolute value, and for greatest
00032 *  accuracy, it should not be much smaller than that.
00033 *
00034 *  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
00035 *  Matrix", Report CS41, Computer Science Dept., Stanford
00036 *  University, July 21, 1966.
00037 *
00038 *  Arguments
00039 *  =========
00040 *
00041 *  RANGE   (input) CHARACTER*1
00042 *          = 'A': ("All")   all eigenvalues will be found.
00043 *          = 'V': ("Value") all eigenvalues in the half-open interval
00044 *                           (VL, VU] will be found.
00045 *          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
00046 *                           entire matrix) will be found.
00047 *
00048 *  ORDER   (input) CHARACTER*1
00049 *          = 'B': ("By Block") the eigenvalues will be grouped by
00050 *                              split-off block (see IBLOCK, ISPLIT) and
00051 *                              ordered from smallest to largest within
00052 *                              the block.
00053 *          = 'E': ("Entire matrix")
00054 *                              the eigenvalues for the entire matrix
00055 *                              will be ordered from smallest to
00056 *                              largest.
00057 *
00058 *  N       (input) INTEGER
00059 *          The order of the tridiagonal matrix T.  N >= 0.
00060 *
00061 *  VL      (input) DOUBLE PRECISION
00062 *  VU      (input) DOUBLE PRECISION
00063 *          If RANGE='V', the lower and upper bounds of the interval to
00064 *          be searched for eigenvalues.  Eigenvalues less than or equal
00065 *          to VL, or greater than VU, will not be returned.  VL < VU.
00066 *          Not referenced if RANGE = 'A' or 'I'.
00067 *
00068 *  IL      (input) INTEGER
00069 *  IU      (input) INTEGER
00070 *          If RANGE='I', the indices (in ascending order) of the
00071 *          smallest and largest eigenvalues to be returned.
00072 *          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
00073 *          Not referenced if RANGE = 'A' or 'V'.
00074 *
00075 *  ABSTOL  (input) DOUBLE PRECISION
00076 *          The absolute tolerance for the eigenvalues.  An eigenvalue
00077 *          (or cluster) is considered to be located if it has been
00078 *          determined to lie in an interval whose width is ABSTOL or
00079 *          less.  If ABSTOL is less than or equal to zero, then ULP*|T|
00080 *          will be used, where |T| means the 1-norm of T.
00081 *
00082 *          Eigenvalues will be computed most accurately when ABSTOL is
00083 *          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
00084 *
00085 *  D       (input) DOUBLE PRECISION array, dimension (N)
00086 *          The n diagonal elements of the tridiagonal matrix T.
00087 *
00088 *  E       (input) DOUBLE PRECISION array, dimension (N-1)
00089 *          The (n-1) off-diagonal elements of the tridiagonal matrix T.
00090 *
00091 *  M       (output) INTEGER
00092 *          The actual number of eigenvalues found. 0 <= M <= N.
00093 *          (See also the description of INFO=2,3.)
00094 *
00095 *  NSPLIT  (output) INTEGER
00096 *          The number of diagonal blocks in the matrix T.
00097 *          1 <= NSPLIT <= N.
00098 *
00099 *  W       (output) DOUBLE PRECISION array, dimension (N)
00100 *          On exit, the first M elements of W will contain the
00101 *          eigenvalues.  (DSTEBZ may use the remaining N-M elements as
00102 *          workspace.)
00103 *
00104 *  IBLOCK  (output) INTEGER array, dimension (N)
00105 *          At each row/column j where E(j) is zero or small, the
00106 *          matrix T is considered to split into a block diagonal
00107 *          matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which
00108 *          block (from 1 to the number of blocks) the eigenvalue W(i)
00109 *          belongs.  (DSTEBZ may use the remaining N-M elements as
00110 *          workspace.)
00111 *
00112 *  ISPLIT  (output) INTEGER array, dimension (N)
00113 *          The splitting points, at which T breaks up into submatrices.
00114 *          The first submatrix consists of rows/columns 1 to ISPLIT(1),
00115 *          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
00116 *          etc., and the NSPLIT-th consists of rows/columns
00117 *          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
00118 *          (Only the first NSPLIT elements will actually be used, but
00119 *          since the user cannot know a priori what value NSPLIT will
00120 *          have, N words must be reserved for ISPLIT.)
00121 *
00122 *  WORK    (workspace) DOUBLE PRECISION array, dimension (4*N)
00123 *
00124 *  IWORK   (workspace) INTEGER array, dimension (3*N)
00125 *
00126 *  INFO    (output) INTEGER
00127 *          = 0:  successful exit
00128 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00129 *          > 0:  some or all of the eigenvalues failed to converge or
00130 *                were not computed:
00131 *                =1 or 3: Bisection failed to converge for some
00132 *                        eigenvalues; these eigenvalues are flagged by a
00133 *                        negative block number.  The effect is that the
00134 *                        eigenvalues may not be as accurate as the
00135 *                        absolute and relative tolerances.  This is
00136 *                        generally caused by unexpectedly inaccurate
00137 *                        arithmetic.
00138 *                =2 or 3: RANGE='I' only: Not all of the eigenvalues
00139 *                        IL:IU were found.
00140 *                        Effect: M < IU+1-IL
00141 *                        Cause:  non-monotonic arithmetic, causing the
00142 *                                Sturm sequence to be non-monotonic.
00143 *                        Cure:   recalculate, using RANGE='A', and pick
00144 *                                out eigenvalues IL:IU.  In some cases,
00145 *                                increasing the PARAMETER "FUDGE" may
00146 *                                make things work.
00147 *                = 4:    RANGE='I', and the Gershgorin interval
00148 *                        initially used was too small.  No eigenvalues
00149 *                        were computed.
00150 *                        Probable cause: your machine has sloppy
00151 *                                        floating-point arithmetic.
00152 *                        Cure: Increase the PARAMETER "FUDGE",
00153 *                              recompile, and try again.
00154 *
00155 *  Internal Parameters
00156 *  ===================
00157 *
00158 *  RELFAC  DOUBLE PRECISION, default = 2.0e0
00159 *          The relative tolerance.  An interval (a,b] lies within
00160 *          "relative tolerance" if  b-a < RELFAC*ulp*max(|a|,|b|),
00161 *          where "ulp" is the machine precision (distance from 1 to
00162 *          the next larger floating point number.)
00163 *
00164 *  FUDGE   DOUBLE PRECISION, default = 2
00165 *          A "fudge factor" to widen the Gershgorin intervals.  Ideally,
00166 *          a value of 1 should work, but on machines with sloppy
00167 *          arithmetic, this needs to be larger.  The default for
00168 *          publicly released versions should be large enough to handle
00169 *          the worst machine around.  Note that this has no effect
00170 *          on accuracy of the solution.
00171 *
00172 *  =====================================================================
00173 *
00174 *     .. Parameters ..
00175       DOUBLE PRECISION   ZERO, ONE, TWO, HALF
00176       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
00177      $                   HALF = 1.0D0 / TWO )
00178       DOUBLE PRECISION   FUDGE, RELFAC
00179       PARAMETER          ( FUDGE = 2.1D0, RELFAC = 2.0D0 )
00180 *     ..
00181 *     .. Local Scalars ..
00182       LOGICAL            NCNVRG, TOOFEW
00183       INTEGER            IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
00184      $                   IM, IN, IOFF, IORDER, IOUT, IRANGE, ITMAX,
00185      $                   ITMP1, IW, IWOFF, J, JB, JDISC, JE, NB, NWL,
00186      $                   NWU
00187       DOUBLE PRECISION   ATOLI, BNORM, GL, GU, PIVMIN, RTOLI, SAFEMN,
00188      $                   TMP1, TMP2, TNORM, ULP, WKILL, WL, WLU, WU, WUL
00189 *     ..
00190 *     .. Local Arrays ..
00191       INTEGER            IDUMMA( 1 )
00192 *     ..
00193 *     .. External Functions ..
00194       LOGICAL            LSAME
00195       INTEGER            ILAENV
00196       DOUBLE PRECISION   DLAMCH
00197       EXTERNAL           LSAME, ILAENV, DLAMCH
00198 *     ..
00199 *     .. External Subroutines ..
00200       EXTERNAL           DLAEBZ, XERBLA
00201 *     ..
00202 *     .. Intrinsic Functions ..
00203       INTRINSIC          ABS, INT, LOG, MAX, MIN, SQRT
00204 *     ..
00205 *     .. Executable Statements ..
00206 *
00207       INFO = 0
00208 *
00209 *     Decode RANGE
00210 *
00211       IF( LSAME( RANGE, 'A' ) ) THEN
00212          IRANGE = 1
00213       ELSE IF( LSAME( RANGE, 'V' ) ) THEN
00214          IRANGE = 2
00215       ELSE IF( LSAME( RANGE, 'I' ) ) THEN
00216          IRANGE = 3
00217       ELSE
00218          IRANGE = 0
00219       END IF
00220 *
00221 *     Decode ORDER
00222 *
00223       IF( LSAME( ORDER, 'B' ) ) THEN
00224          IORDER = 2
00225       ELSE IF( LSAME( ORDER, 'E' ) ) THEN
00226          IORDER = 1
00227       ELSE
00228          IORDER = 0
00229       END IF
00230 *
00231 *     Check for Errors
00232 *
00233       IF( IRANGE.LE.0 ) THEN
00234          INFO = -1
00235       ELSE IF( IORDER.LE.0 ) THEN
00236          INFO = -2
00237       ELSE IF( N.LT.0 ) THEN
00238          INFO = -3
00239       ELSE IF( IRANGE.EQ.2 ) THEN
00240          IF( VL.GE.VU )
00241      $      INFO = -5
00242       ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) )
00243      $          THEN
00244          INFO = -6
00245       ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) )
00246      $          THEN
00247          INFO = -7
00248       END IF
00249 *
00250       IF( INFO.NE.0 ) THEN
00251          CALL XERBLA( 'DSTEBZ', -INFO )
00252          RETURN
00253       END IF
00254 *
00255 *     Initialize error flags
00256 *
00257       INFO = 0
00258       NCNVRG = .FALSE.
00259       TOOFEW = .FALSE.
00260 *
00261 *     Quick return if possible
00262 *
00263       M = 0
00264       IF( N.EQ.0 )
00265      $   RETURN
00266 *
00267 *     Simplifications:
00268 *
00269       IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N )
00270      $   IRANGE = 1
00271 *
00272 *     Get machine constants
00273 *     NB is the minimum vector length for vector bisection, or 0
00274 *     if only scalar is to be done.
00275 *
00276       SAFEMN = DLAMCH( 'S' )
00277       ULP = DLAMCH( 'P' )
00278       RTOLI = ULP*RELFAC
00279       NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 )
00280       IF( NB.LE.1 )
00281      $   NB = 0
00282 *
00283 *     Special Case when N=1
00284 *
00285       IF( N.EQ.1 ) THEN
00286          NSPLIT = 1
00287          ISPLIT( 1 ) = 1
00288          IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN
00289             M = 0
00290          ELSE
00291             W( 1 ) = D( 1 )
00292             IBLOCK( 1 ) = 1
00293             M = 1
00294          END IF
00295          RETURN
00296       END IF
00297 *
00298 *     Compute Splitting Points
00299 *
00300       NSPLIT = 1
00301       WORK( N ) = ZERO
00302       PIVMIN = ONE
00303 *
00304 *DIR$ NOVECTOR
00305       DO 10 J = 2, N
00306          TMP1 = E( J-1 )**2
00307          IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN
00308             ISPLIT( NSPLIT ) = J - 1
00309             NSPLIT = NSPLIT + 1
00310             WORK( J-1 ) = ZERO
00311          ELSE
00312             WORK( J-1 ) = TMP1
00313             PIVMIN = MAX( PIVMIN, TMP1 )
00314          END IF
00315    10 CONTINUE
00316       ISPLIT( NSPLIT ) = N
00317       PIVMIN = PIVMIN*SAFEMN
00318 *
00319 *     Compute Interval and ATOLI
00320 *
00321       IF( IRANGE.EQ.3 ) THEN
00322 *
00323 *        RANGE='I': Compute the interval containing eigenvalues
00324 *                   IL through IU.
00325 *
00326 *        Compute Gershgorin interval for entire (split) matrix
00327 *        and use it as the initial interval
00328 *
00329          GU = D( 1 )
00330          GL = D( 1 )
00331          TMP1 = ZERO
00332 *
00333          DO 20 J = 1, N - 1
00334             TMP2 = SQRT( WORK( J ) )
00335             GU = MAX( GU, D( J )+TMP1+TMP2 )
00336             GL = MIN( GL, D( J )-TMP1-TMP2 )
00337             TMP1 = TMP2
00338    20    CONTINUE
00339 *
00340          GU = MAX( GU, D( N )+TMP1 )
00341          GL = MIN( GL, D( N )-TMP1 )
00342          TNORM = MAX( ABS( GL ), ABS( GU ) )
00343          GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
00344          GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
00345 *
00346 *        Compute Iteration parameters
00347 *
00348          ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
00349      $           LOG( TWO ) ) + 2
00350          IF( ABSTOL.LE.ZERO ) THEN
00351             ATOLI = ULP*TNORM
00352          ELSE
00353             ATOLI = ABSTOL
00354          END IF
00355 *
00356          WORK( N+1 ) = GL
00357          WORK( N+2 ) = GL
00358          WORK( N+3 ) = GU
00359          WORK( N+4 ) = GU
00360          WORK( N+5 ) = GL
00361          WORK( N+6 ) = GU
00362          IWORK( 1 ) = -1
00363          IWORK( 2 ) = -1
00364          IWORK( 3 ) = N + 1
00365          IWORK( 4 ) = N + 1
00366          IWORK( 5 ) = IL - 1
00367          IWORK( 6 ) = IU
00368 *
00369          CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E,
00370      $                WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
00371      $                IWORK, W, IBLOCK, IINFO )
00372 *
00373          IF( IWORK( 6 ).EQ.IU ) THEN
00374             WL = WORK( N+1 )
00375             WLU = WORK( N+3 )
00376             NWL = IWORK( 1 )
00377             WU = WORK( N+4 )
00378             WUL = WORK( N+2 )
00379             NWU = IWORK( 4 )
00380          ELSE
00381             WL = WORK( N+2 )
00382             WLU = WORK( N+4 )
00383             NWL = IWORK( 2 )
00384             WU = WORK( N+3 )
00385             WUL = WORK( N+1 )
00386             NWU = IWORK( 3 )
00387          END IF
00388 *
00389          IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN
00390             INFO = 4
00391             RETURN
00392          END IF
00393       ELSE
00394 *
00395 *        RANGE='A' or 'V' -- Set ATOLI
00396 *
00397          TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ),
00398      $           ABS( D( N ) )+ABS( E( N-1 ) ) )
00399 *
00400          DO 30 J = 2, N - 1
00401             TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+
00402      $              ABS( E( J ) ) )
00403    30    CONTINUE
00404 *
00405          IF( ABSTOL.LE.ZERO ) THEN
00406             ATOLI = ULP*TNORM
00407          ELSE
00408             ATOLI = ABSTOL
00409          END IF
00410 *
00411          IF( IRANGE.EQ.2 ) THEN
00412             WL = VL
00413             WU = VU
00414          ELSE
00415             WL = ZERO
00416             WU = ZERO
00417          END IF
00418       END IF
00419 *
00420 *     Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU.
00421 *     NWL accumulates the number of eigenvalues .le. WL,
00422 *     NWU accumulates the number of eigenvalues .le. WU
00423 *
00424       M = 0
00425       IEND = 0
00426       INFO = 0
00427       NWL = 0
00428       NWU = 0
00429 *
00430       DO 70 JB = 1, NSPLIT
00431          IOFF = IEND
00432          IBEGIN = IOFF + 1
00433          IEND = ISPLIT( JB )
00434          IN = IEND - IOFF
00435 *
00436          IF( IN.EQ.1 ) THEN
00437 *
00438 *           Special Case -- IN=1
00439 *
00440             IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN )
00441      $         NWL = NWL + 1
00442             IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN )
00443      $         NWU = NWU + 1
00444             IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE.
00445      $          D( IBEGIN )-PIVMIN ) ) THEN
00446                M = M + 1
00447                W( M ) = D( IBEGIN )
00448                IBLOCK( M ) = JB
00449             END IF
00450          ELSE
00451 *
00452 *           General Case -- IN > 1
00453 *
00454 *           Compute Gershgorin Interval
00455 *           and use it as the initial interval
00456 *
00457             GU = D( IBEGIN )
00458             GL = D( IBEGIN )
00459             TMP1 = ZERO
00460 *
00461             DO 40 J = IBEGIN, IEND - 1
00462                TMP2 = ABS( E( J ) )
00463                GU = MAX( GU, D( J )+TMP1+TMP2 )
00464                GL = MIN( GL, D( J )-TMP1-TMP2 )
00465                TMP1 = TMP2
00466    40       CONTINUE
00467 *
00468             GU = MAX( GU, D( IEND )+TMP1 )
00469             GL = MIN( GL, D( IEND )-TMP1 )
00470             BNORM = MAX( ABS( GL ), ABS( GU ) )
00471             GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN
00472             GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN
00473 *
00474 *           Compute ATOLI for the current submatrix
00475 *
00476             IF( ABSTOL.LE.ZERO ) THEN
00477                ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) )
00478             ELSE
00479                ATOLI = ABSTOL
00480             END IF
00481 *
00482             IF( IRANGE.GT.1 ) THEN
00483                IF( GU.LT.WL ) THEN
00484                   NWL = NWL + IN
00485                   NWU = NWU + IN
00486                   GO TO 70
00487                END IF
00488                GL = MAX( GL, WL )
00489                GU = MIN( GU, WU )
00490                IF( GL.GE.GU )
00491      $            GO TO 70
00492             END IF
00493 *
00494 *           Set Up Initial Interval
00495 *
00496             WORK( N+1 ) = GL
00497             WORK( N+IN+1 ) = GU
00498             CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
00499      $                   D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
00500      $                   IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
00501      $                   IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
00502 *
00503             NWL = NWL + IWORK( 1 )
00504             NWU = NWU + IWORK( IN+1 )
00505             IWOFF = M - IWORK( 1 )
00506 *
00507 *           Compute Eigenvalues
00508 *
00509             ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) /
00510      $              LOG( TWO ) ) + 2
00511             CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
00512      $                   D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
00513      $                   IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
00514      $                   IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
00515 *
00516 *           Copy Eigenvalues Into W and IBLOCK
00517 *           Use -JB for block number for unconverged eigenvalues.
00518 *
00519             DO 60 J = 1, IOUT
00520                TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) )
00521 *
00522 *              Flag non-convergence.
00523 *
00524                IF( J.GT.IOUT-IINFO ) THEN
00525                   NCNVRG = .TRUE.
00526                   IB = -JB
00527                ELSE
00528                   IB = JB
00529                END IF
00530                DO 50 JE = IWORK( J ) + 1 + IWOFF,
00531      $                 IWORK( J+IN ) + IWOFF
00532                   W( JE ) = TMP1
00533                   IBLOCK( JE ) = IB
00534    50          CONTINUE
00535    60       CONTINUE
00536 *
00537             M = M + IM
00538          END IF
00539    70 CONTINUE
00540 *
00541 *     If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
00542 *     If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
00543 *
00544       IF( IRANGE.EQ.3 ) THEN
00545          IM = 0
00546          IDISCL = IL - 1 - NWL
00547          IDISCU = NWU - IU
00548 *
00549          IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
00550             DO 80 JE = 1, M
00551                IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN
00552                   IDISCL = IDISCL - 1
00553                ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN
00554                   IDISCU = IDISCU - 1
00555                ELSE
00556                   IM = IM + 1
00557                   W( IM ) = W( JE )
00558                   IBLOCK( IM ) = IBLOCK( JE )
00559                END IF
00560    80       CONTINUE
00561             M = IM
00562          END IF
00563          IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
00564 *
00565 *           Code to deal with effects of bad arithmetic:
00566 *           Some low eigenvalues to be discarded are not in (WL,WLU],
00567 *           or high eigenvalues to be discarded are not in (WUL,WU]
00568 *           so just kill off the smallest IDISCL/largest IDISCU
00569 *           eigenvalues, by simply finding the smallest/largest
00570 *           eigenvalue(s).
00571 *
00572 *           (If N(w) is monotone non-decreasing, this should never
00573 *               happen.)
00574 *
00575             IF( IDISCL.GT.0 ) THEN
00576                WKILL = WU
00577                DO 100 JDISC = 1, IDISCL
00578                   IW = 0
00579                   DO 90 JE = 1, M
00580                      IF( IBLOCK( JE ).NE.0 .AND.
00581      $                   ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN
00582                         IW = JE
00583                         WKILL = W( JE )
00584                      END IF
00585    90             CONTINUE
00586                   IBLOCK( IW ) = 0
00587   100          CONTINUE
00588             END IF
00589             IF( IDISCU.GT.0 ) THEN
00590 *
00591                WKILL = WL
00592                DO 120 JDISC = 1, IDISCU
00593                   IW = 0
00594                   DO 110 JE = 1, M
00595                      IF( IBLOCK( JE ).NE.0 .AND.
00596      $                   ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN
00597                         IW = JE
00598                         WKILL = W( JE )
00599                      END IF
00600   110             CONTINUE
00601                   IBLOCK( IW ) = 0
00602   120          CONTINUE
00603             END IF
00604             IM = 0
00605             DO 130 JE = 1, M
00606                IF( IBLOCK( JE ).NE.0 ) THEN
00607                   IM = IM + 1
00608                   W( IM ) = W( JE )
00609                   IBLOCK( IM ) = IBLOCK( JE )
00610                END IF
00611   130       CONTINUE
00612             M = IM
00613          END IF
00614          IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN
00615             TOOFEW = .TRUE.
00616          END IF
00617       END IF
00618 *
00619 *     If ORDER='B', do nothing -- the eigenvalues are already sorted
00620 *        by block.
00621 *     If ORDER='E', sort the eigenvalues from smallest to largest
00622 *
00623       IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN
00624          DO 150 JE = 1, M - 1
00625             IE = 0
00626             TMP1 = W( JE )
00627             DO 140 J = JE + 1, M
00628                IF( W( J ).LT.TMP1 ) THEN
00629                   IE = J
00630                   TMP1 = W( J )
00631                END IF
00632   140       CONTINUE
00633 *
00634             IF( IE.NE.0 ) THEN
00635                ITMP1 = IBLOCK( IE )
00636                W( IE ) = W( JE )
00637                IBLOCK( IE ) = IBLOCK( JE )
00638                W( JE ) = TMP1
00639                IBLOCK( JE ) = ITMP1
00640             END IF
00641   150    CONTINUE
00642       END IF
00643 *
00644       INFO = 0
00645       IF( NCNVRG )
00646      $   INFO = INFO + 1
00647       IF( TOOFEW )
00648      $   INFO = INFO + 2
00649       RETURN
00650 *
00651 *     End of DSTEBZ
00652 *
00653       END