LAPACK: SRC/slarre.f Source File

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