![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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