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