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