LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slarre.f
Go to the documentation of this file.
00001 *> \brief \b SLARRE
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLARRE + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarre.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarre.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarre.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2,
00022 *                           RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
00023 *                           W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
00024 *                           WORK, IWORK, INFO )
00025 * 
00026 *       .. Scalar Arguments ..
00027 *       CHARACTER          RANGE
00028 *       INTEGER            IL, INFO, IU, M, N, NSPLIT
00029 *       REAL               PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
00030 *       ..
00031 *       .. Array Arguments ..
00032 *       INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * ),
00033 *      $                   INDEXW( * )
00034 *       REAL               D( * ), E( * ), E2( * ), GERS( * ),
00035 *      $                   W( * ),WERR( * ), WGAP( * ), WORK( * )
00036 *       ..
00037 *  
00038 *
00039 *> \par Purpose:
00040 *  =============
00041 *>
00042 *> \verbatim
00043 *>
00044 *> To find the desired eigenvalues of a given real symmetric
00045 *> tridiagonal matrix T, SLARRE sets any "small" off-diagonal
00046 *> elements to zero, and for each unreduced block T_i, it finds
00047 *> (a) a suitable shift at one end of the block's spectrum,
00048 *> (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
00049 *> (c) eigenvalues of each L_i D_i L_i^T.
00050 *> The representations and eigenvalues found are then used by
00051 *> SSTEMR to compute the eigenvectors of T.
00052 *> The accuracy varies depending on whether bisection is used to
00053 *> find a few eigenvalues or the dqds algorithm (subroutine SLASQ2) to
00054 *> conpute all and then discard any unwanted one.
00055 *> As an added benefit, SLARRE also outputs the n
00056 *> Gerschgorin intervals for the matrices L_i D_i L_i^T.
00057 *> \endverbatim
00058 *
00059 *  Arguments:
00060 *  ==========
00061 *
00062 *> \param[in] RANGE
00063 *> \verbatim
00064 *>          RANGE is CHARACTER*1
00065 *>          = 'A': ("All")   all eigenvalues will be found.
00066 *>          = 'V': ("Value") all eigenvalues in the half-open interval
00067 *>                           (VL, VU] will be found.
00068 *>          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
00069 *>                           entire matrix) will be found.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] N
00073 *> \verbatim
00074 *>          N is INTEGER
00075 *>          The order of the matrix. N > 0.
00076 *> \endverbatim
00077 *>
00078 *> \param[in,out] VL
00079 *> \verbatim
00080 *>          VL is REAL
00081 *> \endverbatim
00082 *>
00083 *> \param[in,out] VU
00084 *> \verbatim
00085 *>          VU is REAL
00086 *>          If RANGE='V', the lower and upper bounds for the eigenvalues.
00087 *>          Eigenvalues less than or equal to VL, or greater than VU,
00088 *>          will not be returned.  VL < VU.
00089 *>          If RANGE='I' or ='A', SLARRE computes bounds on the desired
00090 *>          part of the spectrum.
00091 *> \endverbatim
00092 *>
00093 *> \param[in] IL
00094 *> \verbatim
00095 *>          IL is INTEGER
00096 *> \endverbatim
00097 *>
00098 *> \param[in] IU
00099 *> \verbatim
00100 *>          IU is INTEGER
00101 *>          If RANGE='I', the indices (in ascending order) of the
00102 *>          smallest and largest eigenvalues to be returned.
00103 *>          1 <= IL <= IU <= N.
00104 *> \endverbatim
00105 *>
00106 *> \param[in,out] D
00107 *> \verbatim
00108 *>          D is REAL array, dimension (N)
00109 *>          On entry, the N diagonal elements of the tridiagonal
00110 *>          matrix T.
00111 *>          On exit, the N diagonal elements of the diagonal
00112 *>          matrices D_i.
00113 *> \endverbatim
00114 *>
00115 *> \param[in,out] E
00116 *> \verbatim
00117 *>          E is REAL array, dimension (N)
00118 *>          On entry, the first (N-1) entries contain the subdiagonal
00119 *>          elements of the tridiagonal matrix T; E(N) need not be set.
00120 *>          On exit, E contains the subdiagonal elements of the unit
00121 *>          bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
00122 *>          1 <= I <= NSPLIT, contain the base points sigma_i on output.
00123 *> \endverbatim
00124 *>
00125 *> \param[in,out] E2
00126 *> \verbatim
00127 *>          E2 is REAL array, dimension (N)
00128 *>          On entry, the first (N-1) entries contain the SQUARES of the
00129 *>          subdiagonal elements of the tridiagonal matrix T;
00130 *>          E2(N) need not be set.
00131 *>          On exit, the entries E2( ISPLIT( I ) ),
00132 *>          1 <= I <= NSPLIT, have been set to zero
00133 *> \endverbatim
00134 *>
00135 *> \param[in] RTOL1
00136 *> \verbatim
00137 *>          RTOL1 is REAL
00138 *> \endverbatim
00139 *>
00140 *> \param[in] RTOL2
00141 *> \verbatim
00142 *>          RTOL2 is REAL
00143 *>           Parameters for bisection.
00144 *>           An interval [LEFT,RIGHT] has converged if
00145 *>           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
00146 *> \endverbatim
00147 *>
00148 *> \param[in] SPLTOL
00149 *> \verbatim
00150 *>          SPLTOL is REAL
00151 *>          The threshold for splitting.
00152 *> \endverbatim
00153 *>
00154 *> \param[out] NSPLIT
00155 *> \verbatim
00156 *>          NSPLIT is INTEGER
00157 *>          The number of blocks T splits into. 1 <= NSPLIT <= N.
00158 *> \endverbatim
00159 *>
00160 *> \param[out] ISPLIT
00161 *> \verbatim
00162 *>          ISPLIT is INTEGER array, dimension (N)
00163 *>          The splitting points, at which T breaks up into blocks.
00164 *>          The first block consists of rows/columns 1 to ISPLIT(1),
00165 *>          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
00166 *>          etc., and the NSPLIT-th consists of rows/columns
00167 *>          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
00168 *> \endverbatim
00169 *>
00170 *> \param[out] M
00171 *> \verbatim
00172 *>          M is INTEGER
00173 *>          The total number of eigenvalues (of all L_i D_i L_i^T)
00174 *>          found.
00175 *> \endverbatim
00176 *>
00177 *> \param[out] W
00178 *> \verbatim
00179 *>          W is REAL array, dimension (N)
00180 *>          The first M elements contain the eigenvalues. The
00181 *>          eigenvalues of each of the blocks, L_i D_i L_i^T, are
00182 *>          sorted in ascending order ( SLARRE may use the
00183 *>          remaining N-M elements as workspace).
00184 *> \endverbatim
00185 *>
00186 *> \param[out] WERR
00187 *> \verbatim
00188 *>          WERR is REAL array, dimension (N)
00189 *>          The error bound on the corresponding eigenvalue in W.
00190 *> \endverbatim
00191 *>
00192 *> \param[out] WGAP
00193 *> \verbatim
00194 *>          WGAP is REAL array, dimension (N)
00195 *>          The separation from the right neighbor eigenvalue in W.
00196 *>          The gap is only with respect to the eigenvalues of the same block
00197 *>          as each block has its own representation tree.
00198 *>          Exception: at the right end of a block we store the left gap
00199 *> \endverbatim
00200 *>
00201 *> \param[out] IBLOCK
00202 *> \verbatim
00203 *>          IBLOCK is INTEGER array, dimension (N)
00204 *>          The indices of the blocks (submatrices) associated with the
00205 *>          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
00206 *>          W(i) belongs to the first block from the top, =2 if W(i)
00207 *>          belongs to the second block, etc.
00208 *> \endverbatim
00209 *>
00210 *> \param[out] INDEXW
00211 *> \verbatim
00212 *>          INDEXW is INTEGER array, dimension (N)
00213 *>          The indices of the eigenvalues within each block (submatrix);
00214 *>          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
00215 *>          i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
00216 *> \endverbatim
00217 *>
00218 *> \param[out] GERS
00219 *> \verbatim
00220 *>          GERS is REAL array, dimension (2*N)
00221 *>          The N Gerschgorin intervals (the i-th Gerschgorin interval
00222 *>          is (GERS(2*i-1), GERS(2*i)).
00223 *> \endverbatim
00224 *>
00225 *> \param[out] PIVMIN
00226 *> \verbatim
00227 *>          PIVMIN is REAL
00228 *>          The minimum pivot in the Sturm sequence for T.
00229 *> \endverbatim
00230 *>
00231 *> \param[out] WORK
00232 *> \verbatim
00233 *>          WORK is REAL array, dimension (6*N)
00234 *>          Workspace.
00235 *> \endverbatim
00236 *>
00237 *> \param[out] IWORK
00238 *> \verbatim
00239 *>          IWORK is INTEGER array, dimension (5*N)
00240 *>          Workspace.
00241 *> \endverbatim
00242 *>
00243 *> \param[out] INFO
00244 *> \verbatim
00245 *>          INFO is INTEGER
00246 *>          = 0:  successful exit
00247 *>          > 0:  A problem occured in SLARRE.
00248 *>          < 0:  One of the called subroutines signaled an internal problem.
00249 *>                Needs inspection of the corresponding parameter IINFO
00250 *>                for further information.
00251 *>
00252 *>          =-1:  Problem in SLARRD.
00253 *>          = 2:  No base representation could be found in MAXTRY iterations.
00254 *>                Increasing MAXTRY and recompilation might be a remedy.
00255 *>          =-3:  Problem in SLARRB when computing the refined root
00256 *>                representation for SLASQ2.
00257 *>          =-4:  Problem in SLARRB when preforming bisection on the
00258 *>                desired part of the spectrum.
00259 *>          =-5:  Problem in SLASQ2.
00260 *>          =-6:  Problem in SLASQ2.
00261 *> \endverbatim
00262 *
00263 *  Authors:
00264 *  ========
00265 *
00266 *> \author Univ. of Tennessee 
00267 *> \author Univ. of California Berkeley 
00268 *> \author Univ. of Colorado Denver 
00269 *> \author NAG Ltd. 
00270 *
00271 *> \date November 2011
00272 *
00273 *> \ingroup auxOTHERauxiliary
00274 *
00275 *> \par Further Details:
00276 *  =====================
00277 *>
00278 *> \verbatim
00279 *>
00280 *>  The base representations are required to suffer very little
00281 *>  element growth and consequently define all their eigenvalues to
00282 *>  high relative accuracy.
00283 *> \endverbatim
00284 *
00285 *> \par Contributors:
00286 *  ==================
00287 *>
00288 *>     Beresford Parlett, University of California, Berkeley, USA \n
00289 *>     Jim Demmel, University of California, Berkeley, USA \n
00290 *>     Inderjit Dhillon, University of Texas, Austin, USA \n
00291 *>     Osni Marques, LBNL/NERSC, USA \n
00292 *>     Christof Voemel, University of California, Berkeley, USA \n
00293 *>
00294 *  =====================================================================
00295       SUBROUTINE SLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2,
00296      $                    RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
00297      $                    W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
00298      $                    WORK, IWORK, INFO )
00299 *
00300 *  -- LAPACK auxiliary routine (version 3.4.0) --
00301 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00302 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00303 *     November 2011
00304 *
00305 *     .. Scalar Arguments ..
00306       CHARACTER          RANGE
00307       INTEGER            IL, INFO, IU, M, N, NSPLIT
00308       REAL               PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
00309 *     ..
00310 *     .. Array Arguments ..
00311       INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * ),
00312      $                   INDEXW( * )
00313       REAL               D( * ), E( * ), E2( * ), GERS( * ),
00314      $                   W( * ),WERR( * ), WGAP( * ), WORK( * )
00315 *     ..
00316 *
00317 *  =====================================================================
00318 *
00319 *     .. Parameters ..
00320       REAL               FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
00321      $                   MAXGROWTH, ONE, PERT, TWO, ZERO
00322       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0,
00323      $                     TWO = 2.0E0, FOUR=4.0E0,
00324      $                     HNDRD = 100.0E0,
00325      $                     PERT = 4.0E0,
00326      $                     HALF = ONE/TWO, FOURTH = ONE/FOUR, FAC= HALF,
00327      $                     MAXGROWTH = 64.0E0, FUDGE = 2.0E0 )
00328       INTEGER            MAXTRY, ALLRNG, INDRNG, VALRNG
00329       PARAMETER          ( MAXTRY = 6, ALLRNG = 1, INDRNG = 2,
00330      $                     VALRNG = 3 )
00331 *     ..
00332 *     .. Local Scalars ..
00333       LOGICAL            FORCEB, NOREP, USEDQD
00334       INTEGER            CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO,
00335      $                   IN, INDL, INDU, IRANGE, J, JBLK, MB, MM,
00336      $                   WBEGIN, WEND
00337       REAL               AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
00338      $                   EMAX, EOLD, EPS, GL, GU, ISLEFT, ISRGHT, RTL,
00339      $                   RTOL, S1, S2, SAFMIN, SGNDEF, SIGMA, SPDIAM,
00340      $                   TAU, TMP, TMP1
00341 
00342 
00343 *     ..
00344 *     .. Local Arrays ..
00345       INTEGER            ISEED( 4 )
00346 *     ..
00347 *     .. External Functions ..
00348       LOGICAL            LSAME
00349       REAL                        SLAMCH
00350       EXTERNAL           SLAMCH, LSAME
00351 
00352 *     ..
00353 *     .. External Subroutines ..
00354       EXTERNAL           SCOPY, SLARNV, SLARRA, SLARRB, SLARRC, SLARRD,
00355      $                   SLASQ2
00356 *     ..
00357 *     .. Intrinsic Functions ..
00358       INTRINSIC          ABS, MAX, MIN
00359 
00360 *     ..
00361 *     .. Executable Statements ..
00362 *
00363 
00364       INFO = 0
00365 
00366 *
00367 *     Decode RANGE
00368 *
00369       IF( LSAME( RANGE, 'A' ) ) THEN
00370          IRANGE = ALLRNG
00371       ELSE IF( LSAME( RANGE, 'V' ) ) THEN
00372          IRANGE = VALRNG
00373       ELSE IF( LSAME( RANGE, 'I' ) ) THEN
00374          IRANGE = INDRNG
00375       END IF
00376 
00377       M = 0
00378 
00379 *     Get machine constants
00380       SAFMIN = SLAMCH( 'S' )
00381       EPS = SLAMCH( 'P' )
00382 
00383 *     Set parameters
00384       RTL = HNDRD*EPS
00385 *     If one were ever to ask for less initial precision in BSRTOL,
00386 *     one should keep in mind that for the subset case, the extremal
00387 *     eigenvalues must be at least as accurate as the current setting
00388 *     (eigenvalues in the middle need not as much accuracy)
00389       BSRTOL = SQRT(EPS)*(0.5E-3)
00390 
00391 *     Treat case of 1x1 matrix for quick return
00392       IF( N.EQ.1 ) THEN
00393          IF( (IRANGE.EQ.ALLRNG).OR.
00394      $       ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
00395      $       ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
00396             M = 1
00397             W(1) = D(1)
00398 *           The computation error of the eigenvalue is zero
00399             WERR(1) = ZERO
00400             WGAP(1) = ZERO
00401             IBLOCK( 1 ) = 1
00402             INDEXW( 1 ) = 1
00403             GERS(1) = D( 1 )
00404             GERS(2) = D( 1 )
00405          ENDIF
00406 *        store the shift for the initial RRR, which is zero in this case
00407          E(1) = ZERO
00408          RETURN
00409       END IF
00410 
00411 *     General case: tridiagonal matrix of order > 1
00412 *
00413 *     Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
00414 *     Compute maximum off-diagonal entry and pivmin.
00415       GL = D(1)
00416       GU = D(1)
00417       EOLD = ZERO
00418       EMAX = ZERO
00419       E(N) = ZERO
00420       DO 5 I = 1,N
00421          WERR(I) = ZERO
00422          WGAP(I) = ZERO
00423          EABS = ABS( E(I) )
00424          IF( EABS .GE. EMAX ) THEN
00425             EMAX = EABS
00426          END IF
00427          TMP1 = EABS + EOLD
00428          GERS( 2*I-1) = D(I) - TMP1
00429          GL =  MIN( GL, GERS( 2*I - 1))
00430          GERS( 2*I ) = D(I) + TMP1
00431          GU = MAX( GU, GERS(2*I) )
00432          EOLD  = EABS
00433  5    CONTINUE
00434 *     The minimum pivot allowed in the Sturm sequence for T
00435       PIVMIN = SAFMIN * MAX( ONE, EMAX**2 )
00436 *     Compute spectral diameter. The Gerschgorin bounds give an
00437 *     estimate that is wrong by at most a factor of SQRT(2)
00438       SPDIAM = GU - GL
00439 
00440 *     Compute splitting points
00441       CALL SLARRA( N, D, E, E2, SPLTOL, SPDIAM,
00442      $                    NSPLIT, ISPLIT, IINFO )
00443 
00444 *     Can force use of bisection instead of faster DQDS.
00445 *     Option left in the code for future multisection work.
00446       FORCEB = .FALSE.
00447 
00448 *     Initialize USEDQD, DQDS should be used for ALLRNG unless someone
00449 *     explicitly wants bisection.
00450       USEDQD = (( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB))
00451 
00452       IF( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ) THEN
00453 *        Set interval [VL,VU] that contains all eigenvalues
00454          VL = GL
00455          VU = GU
00456       ELSE
00457 *        We call SLARRD to find crude approximations to the eigenvalues
00458 *        in the desired range. In case IRANGE = INDRNG, we also obtain the
00459 *        interval (VL,VU] that contains all the wanted eigenvalues.
00460 *        An interval [LEFT,RIGHT] has converged if
00461 *        RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
00462 *        SLARRD needs a WORK of size 4*N, IWORK of size 3*N
00463          CALL SLARRD( RANGE, 'B', N, VL, VU, IL, IU, GERS,
00464      $                    BSRTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
00465      $                    MM, W, WERR, VL, VU, IBLOCK, INDEXW,
00466      $                    WORK, IWORK, IINFO )
00467          IF( IINFO.NE.0 ) THEN
00468             INFO = -1
00469             RETURN
00470          ENDIF
00471 *        Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
00472          DO 14 I = MM+1,N
00473             W( I ) = ZERO
00474             WERR( I ) = ZERO
00475             IBLOCK( I ) = 0
00476             INDEXW( I ) = 0
00477  14      CONTINUE
00478       END IF
00479 
00480 
00481 ***
00482 *     Loop over unreduced blocks
00483       IBEGIN = 1
00484       WBEGIN = 1
00485       DO 170 JBLK = 1, NSPLIT
00486          IEND = ISPLIT( JBLK )
00487          IN = IEND - IBEGIN + 1
00488 
00489 *        1 X 1 block
00490          IF( IN.EQ.1 ) THEN
00491             IF( (IRANGE.EQ.ALLRNG).OR.( (IRANGE.EQ.VALRNG).AND.
00492      $         ( D( IBEGIN ).GT.VL ).AND.( D( IBEGIN ).LE.VU ) )
00493      $        .OR. ( (IRANGE.EQ.INDRNG).AND.(IBLOCK(WBEGIN).EQ.JBLK))
00494      $        ) THEN
00495                M = M + 1
00496                W( M ) = D( IBEGIN )
00497                WERR(M) = ZERO
00498 *              The gap for a single block doesn't matter for the later
00499 *              algorithm and is assigned an arbitrary large value
00500                WGAP(M) = ZERO
00501                IBLOCK( M ) = JBLK
00502                INDEXW( M ) = 1
00503                WBEGIN = WBEGIN + 1
00504             ENDIF
00505 *           E( IEND ) holds the shift for the initial RRR
00506             E( IEND ) = ZERO
00507             IBEGIN = IEND + 1
00508             GO TO 170
00509          END IF
00510 *
00511 *        Blocks of size larger than 1x1
00512 *
00513 *        E( IEND ) will hold the shift for the initial RRR, for now set it =0
00514          E( IEND ) = ZERO
00515 *
00516 *        Find local outer bounds GL,GU for the block
00517          GL = D(IBEGIN)
00518          GU = D(IBEGIN)
00519          DO 15 I = IBEGIN , IEND
00520             GL = MIN( GERS( 2*I-1 ), GL )
00521             GU = MAX( GERS( 2*I ), GU )
00522  15      CONTINUE
00523          SPDIAM = GU - GL
00524 
00525          IF(.NOT. ((IRANGE.EQ.ALLRNG).AND.(.NOT.FORCEB)) ) THEN
00526 *           Count the number of eigenvalues in the current block.
00527             MB = 0
00528             DO 20 I = WBEGIN,MM
00529                IF( IBLOCK(I).EQ.JBLK ) THEN
00530                   MB = MB+1
00531                ELSE
00532                   GOTO 21
00533                ENDIF
00534  20         CONTINUE
00535  21         CONTINUE
00536 
00537             IF( MB.EQ.0) THEN
00538 *              No eigenvalue in the current block lies in the desired range
00539 *              E( IEND ) holds the shift for the initial RRR
00540                E( IEND ) = ZERO
00541                IBEGIN = IEND + 1
00542                GO TO 170
00543             ELSE
00544 
00545 *              Decide whether dqds or bisection is more efficient
00546                USEDQD = ( (MB .GT. FAC*IN) .AND. (.NOT.FORCEB) )
00547                WEND = WBEGIN + MB - 1
00548 *              Calculate gaps for the current block
00549 *              In later stages, when representations for individual
00550 *              eigenvalues are different, we use SIGMA = E( IEND ).
00551                SIGMA = ZERO
00552                DO 30 I = WBEGIN, WEND - 1
00553                   WGAP( I ) = MAX( ZERO,
00554      $                        W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
00555  30            CONTINUE
00556                WGAP( WEND ) = MAX( ZERO,
00557      $                     VU - SIGMA - (W( WEND )+WERR( WEND )))
00558 *              Find local index of the first and last desired evalue.
00559                INDL = INDEXW(WBEGIN)
00560                INDU = INDEXW( WEND )
00561             ENDIF
00562          ENDIF
00563          IF(( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ).OR.USEDQD) THEN
00564 *           Case of DQDS
00565 *           Find approximations to the extremal eigenvalues of the block
00566             CALL SLARRK( IN, 1, GL, GU, D(IBEGIN),
00567      $               E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
00568             IF( IINFO.NE.0 ) THEN
00569                INFO = -1
00570                RETURN
00571             ENDIF
00572             ISLEFT = MAX(GL, TMP - TMP1
00573      $               - HNDRD * EPS* ABS(TMP - TMP1))
00574 
00575             CALL SLARRK( IN, IN, GL, GU, D(IBEGIN),
00576      $               E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
00577             IF( IINFO.NE.0 ) THEN
00578                INFO = -1
00579                RETURN
00580             ENDIF
00581             ISRGHT = MIN(GU, TMP + TMP1
00582      $                 + HNDRD * EPS * ABS(TMP + TMP1))
00583 *           Improve the estimate of the spectral diameter
00584             SPDIAM = ISRGHT - ISLEFT
00585          ELSE
00586 *           Case of bisection
00587 *           Find approximations to the wanted extremal eigenvalues
00588             ISLEFT = MAX(GL, W(WBEGIN) - WERR(WBEGIN)
00589      $                  - HNDRD * EPS*ABS(W(WBEGIN)- WERR(WBEGIN) ))
00590             ISRGHT = MIN(GU,W(WEND) + WERR(WEND)
00591      $                  + HNDRD * EPS * ABS(W(WEND)+ WERR(WEND)))
00592          ENDIF
00593 
00594 
00595 *        Decide whether the base representation for the current block
00596 *        L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
00597 *        should be on the left or the right end of the current block.
00598 *        The strategy is to shift to the end which is "more populated"
00599 *        Furthermore, decide whether to use DQDS for the computation of
00600 *        the eigenvalue approximations at the end of SLARRE or bisection.
00601 *        dqds is chosen if all eigenvalues are desired or the number of
00602 *        eigenvalues to be computed is large compared to the blocksize.
00603          IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
00604 *           If all the eigenvalues have to be computed, we use dqd
00605             USEDQD = .TRUE.
00606 *           INDL is the local index of the first eigenvalue to compute
00607             INDL = 1
00608             INDU = IN
00609 *           MB =  number of eigenvalues to compute
00610             MB = IN
00611             WEND = WBEGIN + MB - 1
00612 *           Define 1/4 and 3/4 points of the spectrum
00613             S1 = ISLEFT + FOURTH * SPDIAM
00614             S2 = ISRGHT - FOURTH * SPDIAM
00615          ELSE
00616 *           SLARRD has computed IBLOCK and INDEXW for each eigenvalue
00617 *           approximation.
00618 *           choose sigma
00619             IF( USEDQD ) THEN
00620                S1 = ISLEFT + FOURTH * SPDIAM
00621                S2 = ISRGHT - FOURTH * SPDIAM
00622             ELSE
00623                TMP = MIN(ISRGHT,VU) -  MAX(ISLEFT,VL)
00624                S1 =  MAX(ISLEFT,VL) + FOURTH * TMP
00625                S2 =  MIN(ISRGHT,VU) - FOURTH * TMP
00626             ENDIF
00627          ENDIF
00628 
00629 *        Compute the negcount at the 1/4 and 3/4 points
00630          IF(MB.GT.1) THEN
00631             CALL SLARRC( 'T', IN, S1, S2, D(IBEGIN),
00632      $                    E(IBEGIN), PIVMIN, CNT, CNT1, CNT2, IINFO)
00633          ENDIF
00634 
00635          IF(MB.EQ.1) THEN
00636             SIGMA = GL
00637             SGNDEF = ONE
00638          ELSEIF( CNT1 - INDL .GE. INDU - CNT2 ) THEN
00639             IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
00640                SIGMA = MAX(ISLEFT,GL)
00641             ELSEIF( USEDQD ) THEN
00642 *              use Gerschgorin bound as shift to get pos def matrix
00643 *              for dqds
00644                SIGMA = ISLEFT
00645             ELSE
00646 *              use approximation of the first desired eigenvalue of the
00647 *              block as shift
00648                SIGMA = MAX(ISLEFT,VL)
00649             ENDIF
00650             SGNDEF = ONE
00651          ELSE
00652             IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
00653                SIGMA = MIN(ISRGHT,GU)
00654             ELSEIF( USEDQD ) THEN
00655 *              use Gerschgorin bound as shift to get neg def matrix
00656 *              for dqds
00657                SIGMA = ISRGHT
00658             ELSE
00659 *              use approximation of the first desired eigenvalue of the
00660 *              block as shift
00661                SIGMA = MIN(ISRGHT,VU)
00662             ENDIF
00663             SGNDEF = -ONE
00664          ENDIF
00665 
00666 
00667 *        An initial SIGMA has been chosen that will be used for computing
00668 *        T - SIGMA I = L D L^T
00669 *        Define the increment TAU of the shift in case the initial shift
00670 *        needs to be refined to obtain a factorization with not too much
00671 *        element growth.
00672          IF( USEDQD ) THEN
00673 *           The initial SIGMA was to the outer end of the spectrum
00674 *           the matrix is definite and we need not retreat.
00675             TAU = SPDIAM*EPS*N + TWO*PIVMIN
00676             TAU = MAX( TAU,TWO*EPS*ABS(SIGMA) )
00677          ELSE
00678             IF(MB.GT.1) THEN
00679                CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN)
00680                AVGAP = ABS(CLWDTH / REAL(WEND-WBEGIN))
00681                IF( SGNDEF.EQ.ONE ) THEN
00682                   TAU = HALF*MAX(WGAP(WBEGIN),AVGAP)
00683                   TAU = MAX(TAU,WERR(WBEGIN))
00684                ELSE
00685                   TAU = HALF*MAX(WGAP(WEND-1),AVGAP)
00686                   TAU = MAX(TAU,WERR(WEND))
00687                ENDIF
00688             ELSE
00689                TAU = WERR(WBEGIN)
00690             ENDIF
00691          ENDIF
00692 *
00693          DO 80 IDUM = 1, MAXTRY
00694 *           Compute L D L^T factorization of tridiagonal matrix T - sigma I.
00695 *           Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
00696 *           pivots in WORK(2*IN+1:3*IN)
00697             DPIVOT = D( IBEGIN ) - SIGMA
00698             WORK( 1 ) = DPIVOT
00699             DMAX = ABS( WORK(1) )
00700             J = IBEGIN
00701             DO 70 I = 1, IN - 1
00702                WORK( 2*IN+I ) = ONE / WORK( I )
00703                TMP = E( J )*WORK( 2*IN+I )
00704                WORK( IN+I ) = TMP
00705                DPIVOT = ( D( J+1 )-SIGMA ) - TMP*E( J )
00706                WORK( I+1 ) = DPIVOT
00707                DMAX = MAX( DMAX, ABS(DPIVOT) )
00708                J = J + 1
00709  70         CONTINUE
00710 *           check for element growth
00711             IF( DMAX .GT. MAXGROWTH*SPDIAM ) THEN
00712                NOREP = .TRUE.
00713             ELSE
00714                NOREP = .FALSE.
00715             ENDIF
00716             IF( USEDQD .AND. .NOT.NOREP ) THEN
00717 *              Ensure the definiteness of the representation
00718 *              All entries of D (of L D L^T) must have the same sign
00719                DO 71 I = 1, IN
00720                   TMP = SGNDEF*WORK( I )
00721                   IF( TMP.LT.ZERO ) NOREP = .TRUE.
00722  71            CONTINUE
00723             ENDIF
00724             IF(NOREP) THEN
00725 *              Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
00726 *              shift which makes the matrix definite. So we should end up
00727 *              here really only in the case of IRANGE = VALRNG or INDRNG.
00728                IF( IDUM.EQ.MAXTRY-1 ) THEN
00729                   IF( SGNDEF.EQ.ONE ) THEN
00730 *                    The fudged Gerschgorin shift should succeed
00731                      SIGMA =
00732      $                    GL - FUDGE*SPDIAM*EPS*N - FUDGE*TWO*PIVMIN
00733                   ELSE
00734                      SIGMA =
00735      $                    GU + FUDGE*SPDIAM*EPS*N + FUDGE*TWO*PIVMIN
00736                   END IF
00737                ELSE
00738                   SIGMA = SIGMA - SGNDEF * TAU
00739                   TAU = TWO * TAU
00740                END IF
00741             ELSE
00742 *              an initial RRR is found
00743                GO TO 83
00744             END IF
00745  80      CONTINUE
00746 *        if the program reaches this point, no base representation could be
00747 *        found in MAXTRY iterations.
00748          INFO = 2
00749          RETURN
00750 
00751  83      CONTINUE
00752 *        At this point, we have found an initial base representation
00753 *        T - SIGMA I = L D L^T with not too much element growth.
00754 *        Store the shift.
00755          E( IEND ) = SIGMA
00756 *        Store D and L.
00757          CALL SCOPY( IN, WORK, 1, D( IBEGIN ), 1 )
00758          CALL SCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 )
00759 
00760 
00761          IF(MB.GT.1 ) THEN
00762 *
00763 *           Perturb each entry of the base representation by a small
00764 *           (but random) relative amount to overcome difficulties with
00765 *           glued matrices.
00766 *
00767             DO 122 I = 1, 4
00768                ISEED( I ) = 1
00769  122        CONTINUE
00770 
00771             CALL SLARNV(2, ISEED, 2*IN-1, WORK(1))
00772             DO 125 I = 1,IN-1
00773                D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(I))
00774                E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(IN+I))
00775  125        CONTINUE
00776             D(IEND) = D(IEND)*(ONE+EPS*FOUR*WORK(IN))
00777 *
00778          ENDIF
00779 *
00780 *        Don't update the Gerschgorin intervals because keeping track
00781 *        of the updates would be too much work in SLARRV.
00782 *        We update W instead and use it to locate the proper Gerschgorin
00783 *        intervals.
00784 
00785 *        Compute the required eigenvalues of L D L' by bisection or dqds
00786          IF ( .NOT.USEDQD ) THEN
00787 *           If SLARRD has been used, shift the eigenvalue approximations
00788 *           according to their representation. This is necessary for
00789 *           a uniform SLARRV since dqds computes eigenvalues of the
00790 *           shifted representation. In SLARRV, W will always hold the
00791 *           UNshifted eigenvalue approximation.
00792             DO 134 J=WBEGIN,WEND
00793                W(J) = W(J) - SIGMA
00794                WERR(J) = WERR(J) + ABS(W(J)) * EPS
00795  134        CONTINUE
00796 *           call SLARRB to reduce eigenvalue error of the approximations
00797 *           from SLARRD
00798             DO 135 I = IBEGIN, IEND-1
00799                WORK( I ) = D( I ) * E( I )**2
00800  135        CONTINUE
00801 *           use bisection to find EV from INDL to INDU
00802             CALL SLARRB(IN, D(IBEGIN), WORK(IBEGIN),
00803      $                  INDL, INDU, RTOL1, RTOL2, INDL-1,
00804      $                  W(WBEGIN), WGAP(WBEGIN), WERR(WBEGIN),
00805      $                  WORK( 2*N+1 ), IWORK, PIVMIN, SPDIAM,
00806      $                  IN, IINFO )
00807             IF( IINFO .NE. 0 ) THEN
00808                INFO = -4
00809                RETURN
00810             END IF
00811 *           SLARRB computes all gaps correctly except for the last one
00812 *           Record distance to VU/GU
00813             WGAP( WEND ) = MAX( ZERO,
00814      $           ( VU-SIGMA ) - ( W( WEND ) + WERR( WEND ) ) )
00815             DO 138 I = INDL, INDU
00816                M = M + 1
00817                IBLOCK(M) = JBLK
00818                INDEXW(M) = I
00819  138        CONTINUE
00820          ELSE
00821 *           Call dqds to get all eigs (and then possibly delete unwanted
00822 *           eigenvalues).
00823 *           Note that dqds finds the eigenvalues of the L D L^T representation
00824 *           of T to high relative accuracy. High relative accuracy
00825 *           might be lost when the shift of the RRR is subtracted to obtain
00826 *           the eigenvalues of T. However, T is not guaranteed to define its
00827 *           eigenvalues to high relative accuracy anyway.
00828 *           Set RTOL to the order of the tolerance used in SLASQ2
00829 *           This is an ESTIMATED error, the worst case bound is 4*N*EPS
00830 *           which is usually too large and requires unnecessary work to be
00831 *           done by bisection when computing the eigenvectors
00832             RTOL = LOG(REAL(IN)) * FOUR * EPS
00833             J = IBEGIN
00834             DO 140 I = 1, IN - 1
00835                WORK( 2*I-1 ) = ABS( D( J ) )
00836                WORK( 2*I ) = E( J )*E( J )*WORK( 2*I-1 )
00837                J = J + 1
00838   140       CONTINUE
00839             WORK( 2*IN-1 ) = ABS( D( IEND ) )
00840             WORK( 2*IN ) = ZERO
00841             CALL SLASQ2( IN, WORK, IINFO )
00842             IF( IINFO .NE. 0 ) THEN
00843 *              If IINFO = -5 then an index is part of a tight cluster
00844 *              and should be changed. The index is in IWORK(1) and the
00845 *              gap is in WORK(N+1)
00846                INFO = -5
00847                RETURN
00848             ELSE
00849 *              Test that all eigenvalues are positive as expected
00850                DO 149 I = 1, IN
00851                   IF( WORK( I ).LT.ZERO ) THEN
00852                      INFO = -6
00853                      RETURN
00854                   ENDIF
00855  149           CONTINUE
00856             END IF
00857             IF( SGNDEF.GT.ZERO ) THEN
00858                DO 150 I = INDL, INDU
00859                   M = M + 1
00860                   W( M ) = WORK( IN-I+1 )
00861                   IBLOCK( M ) = JBLK
00862                   INDEXW( M ) = I
00863  150           CONTINUE
00864             ELSE
00865                DO 160 I = INDL, INDU
00866                   M = M + 1
00867                   W( M ) = -WORK( I )
00868                   IBLOCK( M ) = JBLK
00869                   INDEXW( M ) = I
00870  160           CONTINUE
00871             END IF
00872 
00873             DO 165 I = M - MB + 1, M
00874 *              the value of RTOL below should be the tolerance in SLASQ2
00875                WERR( I ) = RTOL * ABS( W(I) )
00876  165        CONTINUE
00877             DO 166 I = M - MB + 1, M - 1
00878 *              compute the right gap between the intervals
00879                WGAP( I ) = MAX( ZERO,
00880      $                          W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
00881  166        CONTINUE
00882             WGAP( M ) = MAX( ZERO,
00883      $           ( VU-SIGMA ) - ( W( M ) + WERR( M ) ) )
00884          END IF
00885 *        proceed with next block
00886          IBEGIN = IEND + 1
00887          WBEGIN = WEND + 1
00888  170  CONTINUE
00889 *
00890 
00891       RETURN
00892 *
00893 *     end of SLARRE
00894 *
00895       END
 All Files Functions