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