![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLARRE 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLARRE + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarre.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarre.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarre.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLARRE( 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 * DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU 00030 * .. 00031 * .. Array Arguments .. 00032 * INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ), 00033 * $ INDEXW( * ) 00034 * DOUBLE PRECISION 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, DLARRE 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 *> DSTEMR 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 DLASQ2) to 00054 *> conpute all and then discard any unwanted one. 00055 *> As an added benefit, DLARRE 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 DOUBLE PRECISION 00081 *> \endverbatim 00082 *> 00083 *> \param[in,out] VU 00084 *> \verbatim 00085 *> VU is DOUBLE PRECISION 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', DLARRE 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 00138 *> \endverbatim 00139 *> 00140 *> \param[in] RTOL2 00141 *> \verbatim 00142 *> RTOL2 is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 ( DLARRE may use the 00183 *> remaining N-M elements as workspace). 00184 *> \endverbatim 00185 *> 00186 *> \param[out] WERR 00187 *> \verbatim 00188 *> WERR is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 00228 *> The minimum pivot in the Sturm sequence for T. 00229 *> \endverbatim 00230 *> 00231 *> \param[out] WORK 00232 *> \verbatim 00233 *> WORK is DOUBLE PRECISION 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 DLARRE. 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 DLARRD. 00253 *> = 2: No base representation could be found in MAXTRY iterations. 00254 *> Increasing MAXTRY and recompilation might be a remedy. 00255 *> =-3: Problem in DLARRB when computing the refined root 00256 *> representation for DLASQ2. 00257 *> =-4: Problem in DLARRB when preforming bisection on the 00258 *> desired part of the spectrum. 00259 *> =-5: Problem in DLASQ2. 00260 *> =-6: Problem in DLASQ2. 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 DLARRE( 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 DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU 00309 * .. 00310 * .. Array Arguments .. 00311 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ), 00312 $ INDEXW( * ) 00313 DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ), 00314 $ W( * ),WERR( * ), WGAP( * ), WORK( * ) 00315 * .. 00316 * 00317 * ===================================================================== 00318 * 00319 * .. Parameters .. 00320 DOUBLE PRECISION FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD, 00321 $ MAXGROWTH, ONE, PERT, TWO, ZERO 00322 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, 00323 $ TWO = 2.0D0, FOUR=4.0D0, 00324 $ HNDRD = 100.0D0, 00325 $ PERT = 8.0D0, 00326 $ HALF = ONE/TWO, FOURTH = ONE/FOUR, FAC= HALF, 00327 $ MAXGROWTH = 64.0D0, FUDGE = 2.0D0 ) 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 DOUBLE PRECISION 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 DOUBLE PRECISION DLAMCH 00350 EXTERNAL DLAMCH, LSAME 00351 00352 * .. 00353 * .. External Subroutines .. 00354 EXTERNAL DCOPY, DLARNV, DLARRA, DLARRB, DLARRC, DLARRD, 00355 $ DLASQ2 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 = DLAMCH( 'S' ) 00381 EPS = DLAMCH( 'P' ) 00382 00383 * Set parameters 00384 RTL = SQRT(EPS) 00385 BSRTOL = SQRT(EPS) 00386 00387 * Treat case of 1x1 matrix for quick return 00388 IF( N.EQ.1 ) THEN 00389 IF( (IRANGE.EQ.ALLRNG).OR. 00390 $ ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR. 00391 $ ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN 00392 M = 1 00393 W(1) = D(1) 00394 * The computation error of the eigenvalue is zero 00395 WERR(1) = ZERO 00396 WGAP(1) = ZERO 00397 IBLOCK( 1 ) = 1 00398 INDEXW( 1 ) = 1 00399 GERS(1) = D( 1 ) 00400 GERS(2) = D( 1 ) 00401 ENDIF 00402 * store the shift for the initial RRR, which is zero in this case 00403 E(1) = ZERO 00404 RETURN 00405 END IF 00406 00407 * General case: tridiagonal matrix of order > 1 00408 * 00409 * Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter. 00410 * Compute maximum off-diagonal entry and pivmin. 00411 GL = D(1) 00412 GU = D(1) 00413 EOLD = ZERO 00414 EMAX = ZERO 00415 E(N) = ZERO 00416 DO 5 I = 1,N 00417 WERR(I) = ZERO 00418 WGAP(I) = ZERO 00419 EABS = ABS( E(I) ) 00420 IF( EABS .GE. EMAX ) THEN 00421 EMAX = EABS 00422 END IF 00423 TMP1 = EABS + EOLD 00424 GERS( 2*I-1) = D(I) - TMP1 00425 GL = MIN( GL, GERS( 2*I - 1)) 00426 GERS( 2*I ) = D(I) + TMP1 00427 GU = MAX( GU, GERS(2*I) ) 00428 EOLD = EABS 00429 5 CONTINUE 00430 * The minimum pivot allowed in the Sturm sequence for T 00431 PIVMIN = SAFMIN * MAX( ONE, EMAX**2 ) 00432 * Compute spectral diameter. The Gerschgorin bounds give an 00433 * estimate that is wrong by at most a factor of SQRT(2) 00434 SPDIAM = GU - GL 00435 00436 * Compute splitting points 00437 CALL DLARRA( N, D, E, E2, SPLTOL, SPDIAM, 00438 $ NSPLIT, ISPLIT, IINFO ) 00439 00440 * Can force use of bisection instead of faster DQDS. 00441 * Option left in the code for future multisection work. 00442 FORCEB = .FALSE. 00443 00444 * Initialize USEDQD, DQDS should be used for ALLRNG unless someone 00445 * explicitly wants bisection. 00446 USEDQD = (( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB)) 00447 00448 IF( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ) THEN 00449 * Set interval [VL,VU] that contains all eigenvalues 00450 VL = GL 00451 VU = GU 00452 ELSE 00453 * We call DLARRD to find crude approximations to the eigenvalues 00454 * in the desired range. In case IRANGE = INDRNG, we also obtain the 00455 * interval (VL,VU] that contains all the wanted eigenvalues. 00456 * An interval [LEFT,RIGHT] has converged if 00457 * RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT)) 00458 * DLARRD needs a WORK of size 4*N, IWORK of size 3*N 00459 CALL DLARRD( RANGE, 'B', N, VL, VU, IL, IU, GERS, 00460 $ BSRTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, 00461 $ MM, W, WERR, VL, VU, IBLOCK, INDEXW, 00462 $ WORK, IWORK, IINFO ) 00463 IF( IINFO.NE.0 ) THEN 00464 INFO = -1 00465 RETURN 00466 ENDIF 00467 * Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0 00468 DO 14 I = MM+1,N 00469 W( I ) = ZERO 00470 WERR( I ) = ZERO 00471 IBLOCK( I ) = 0 00472 INDEXW( I ) = 0 00473 14 CONTINUE 00474 END IF 00475 00476 00477 *** 00478 * Loop over unreduced blocks 00479 IBEGIN = 1 00480 WBEGIN = 1 00481 DO 170 JBLK = 1, NSPLIT 00482 IEND = ISPLIT( JBLK ) 00483 IN = IEND - IBEGIN + 1 00484 00485 * 1 X 1 block 00486 IF( IN.EQ.1 ) THEN 00487 IF( (IRANGE.EQ.ALLRNG).OR.( (IRANGE.EQ.VALRNG).AND. 00488 $ ( D( IBEGIN ).GT.VL ).AND.( D( IBEGIN ).LE.VU ) ) 00489 $ .OR. ( (IRANGE.EQ.INDRNG).AND.(IBLOCK(WBEGIN).EQ.JBLK)) 00490 $ ) THEN 00491 M = M + 1 00492 W( M ) = D( IBEGIN ) 00493 WERR(M) = ZERO 00494 * The gap for a single block doesn't matter for the later 00495 * algorithm and is assigned an arbitrary large value 00496 WGAP(M) = ZERO 00497 IBLOCK( M ) = JBLK 00498 INDEXW( M ) = 1 00499 WBEGIN = WBEGIN + 1 00500 ENDIF 00501 * E( IEND ) holds the shift for the initial RRR 00502 E( IEND ) = ZERO 00503 IBEGIN = IEND + 1 00504 GO TO 170 00505 END IF 00506 * 00507 * Blocks of size larger than 1x1 00508 * 00509 * E( IEND ) will hold the shift for the initial RRR, for now set it =0 00510 E( IEND ) = ZERO 00511 * 00512 * Find local outer bounds GL,GU for the block 00513 GL = D(IBEGIN) 00514 GU = D(IBEGIN) 00515 DO 15 I = IBEGIN , IEND 00516 GL = MIN( GERS( 2*I-1 ), GL ) 00517 GU = MAX( GERS( 2*I ), GU ) 00518 15 CONTINUE 00519 SPDIAM = GU - GL 00520 00521 IF(.NOT. ((IRANGE.EQ.ALLRNG).AND.(.NOT.FORCEB)) ) THEN 00522 * Count the number of eigenvalues in the current block. 00523 MB = 0 00524 DO 20 I = WBEGIN,MM 00525 IF( IBLOCK(I).EQ.JBLK ) THEN 00526 MB = MB+1 00527 ELSE 00528 GOTO 21 00529 ENDIF 00530 20 CONTINUE 00531 21 CONTINUE 00532 00533 IF( MB.EQ.0) THEN 00534 * No eigenvalue in the current block lies in the desired range 00535 * E( IEND ) holds the shift for the initial RRR 00536 E( IEND ) = ZERO 00537 IBEGIN = IEND + 1 00538 GO TO 170 00539 ELSE 00540 00541 * Decide whether dqds or bisection is more efficient 00542 USEDQD = ( (MB .GT. FAC*IN) .AND. (.NOT.FORCEB) ) 00543 WEND = WBEGIN + MB - 1 00544 * Calculate gaps for the current block 00545 * In later stages, when representations for individual 00546 * eigenvalues are different, we use SIGMA = E( IEND ). 00547 SIGMA = ZERO 00548 DO 30 I = WBEGIN, WEND - 1 00549 WGAP( I ) = MAX( ZERO, 00550 $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) ) 00551 30 CONTINUE 00552 WGAP( WEND ) = MAX( ZERO, 00553 $ VU - SIGMA - (W( WEND )+WERR( WEND ))) 00554 * Find local index of the first and last desired evalue. 00555 INDL = INDEXW(WBEGIN) 00556 INDU = INDEXW( WEND ) 00557 ENDIF 00558 ENDIF 00559 IF(( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ).OR.USEDQD) THEN 00560 * Case of DQDS 00561 * Find approximations to the extremal eigenvalues of the block 00562 CALL DLARRK( IN, 1, GL, GU, D(IBEGIN), 00563 $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO ) 00564 IF( IINFO.NE.0 ) THEN 00565 INFO = -1 00566 RETURN 00567 ENDIF 00568 ISLEFT = MAX(GL, TMP - TMP1 00569 $ - HNDRD * EPS* ABS(TMP - TMP1)) 00570 00571 CALL DLARRK( IN, IN, GL, GU, D(IBEGIN), 00572 $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO ) 00573 IF( IINFO.NE.0 ) THEN 00574 INFO = -1 00575 RETURN 00576 ENDIF 00577 ISRGHT = MIN(GU, TMP + TMP1 00578 $ + HNDRD * EPS * ABS(TMP + TMP1)) 00579 * Improve the estimate of the spectral diameter 00580 SPDIAM = ISRGHT - ISLEFT 00581 ELSE 00582 * Case of bisection 00583 * Find approximations to the wanted extremal eigenvalues 00584 ISLEFT = MAX(GL, W(WBEGIN) - WERR(WBEGIN) 00585 $ - HNDRD * EPS*ABS(W(WBEGIN)- WERR(WBEGIN) )) 00586 ISRGHT = MIN(GU,W(WEND) + WERR(WEND) 00587 $ + HNDRD * EPS * ABS(W(WEND)+ WERR(WEND))) 00588 ENDIF 00589 00590 00591 * Decide whether the base representation for the current block 00592 * L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I 00593 * should be on the left or the right end of the current block. 00594 * The strategy is to shift to the end which is "more populated" 00595 * Furthermore, decide whether to use DQDS for the computation of 00596 * the eigenvalue approximations at the end of DLARRE or bisection. 00597 * dqds is chosen if all eigenvalues are desired or the number of 00598 * eigenvalues to be computed is large compared to the blocksize. 00599 IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN 00600 * If all the eigenvalues have to be computed, we use dqd 00601 USEDQD = .TRUE. 00602 * INDL is the local index of the first eigenvalue to compute 00603 INDL = 1 00604 INDU = IN 00605 * MB = number of eigenvalues to compute 00606 MB = IN 00607 WEND = WBEGIN + MB - 1 00608 * Define 1/4 and 3/4 points of the spectrum 00609 S1 = ISLEFT + FOURTH * SPDIAM 00610 S2 = ISRGHT - FOURTH * SPDIAM 00611 ELSE 00612 * DLARRD has computed IBLOCK and INDEXW for each eigenvalue 00613 * approximation. 00614 * choose sigma 00615 IF( USEDQD ) THEN 00616 S1 = ISLEFT + FOURTH * SPDIAM 00617 S2 = ISRGHT - FOURTH * SPDIAM 00618 ELSE 00619 TMP = MIN(ISRGHT,VU) - MAX(ISLEFT,VL) 00620 S1 = MAX(ISLEFT,VL) + FOURTH * TMP 00621 S2 = MIN(ISRGHT,VU) - FOURTH * TMP 00622 ENDIF 00623 ENDIF 00624 00625 * Compute the negcount at the 1/4 and 3/4 points 00626 IF(MB.GT.1) THEN 00627 CALL DLARRC( 'T', IN, S1, S2, D(IBEGIN), 00628 $ E(IBEGIN), PIVMIN, CNT, CNT1, CNT2, IINFO) 00629 ENDIF 00630 00631 IF(MB.EQ.1) THEN 00632 SIGMA = GL 00633 SGNDEF = ONE 00634 ELSEIF( CNT1 - INDL .GE. INDU - CNT2 ) THEN 00635 IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN 00636 SIGMA = MAX(ISLEFT,GL) 00637 ELSEIF( USEDQD ) THEN 00638 * use Gerschgorin bound as shift to get pos def matrix 00639 * for dqds 00640 SIGMA = ISLEFT 00641 ELSE 00642 * use approximation of the first desired eigenvalue of the 00643 * block as shift 00644 SIGMA = MAX(ISLEFT,VL) 00645 ENDIF 00646 SGNDEF = ONE 00647 ELSE 00648 IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN 00649 SIGMA = MIN(ISRGHT,GU) 00650 ELSEIF( USEDQD ) THEN 00651 * use Gerschgorin bound as shift to get neg def matrix 00652 * for dqds 00653 SIGMA = ISRGHT 00654 ELSE 00655 * use approximation of the first desired eigenvalue of the 00656 * block as shift 00657 SIGMA = MIN(ISRGHT,VU) 00658 ENDIF 00659 SGNDEF = -ONE 00660 ENDIF 00661 00662 00663 * An initial SIGMA has been chosen that will be used for computing 00664 * T - SIGMA I = L D L^T 00665 * Define the increment TAU of the shift in case the initial shift 00666 * needs to be refined to obtain a factorization with not too much 00667 * element growth. 00668 IF( USEDQD ) THEN 00669 * The initial SIGMA was to the outer end of the spectrum 00670 * the matrix is definite and we need not retreat. 00671 TAU = SPDIAM*EPS*N + TWO*PIVMIN 00672 TAU = MAX( TAU,TWO*EPS*ABS(SIGMA) ) 00673 ELSE 00674 IF(MB.GT.1) THEN 00675 CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN) 00676 AVGAP = ABS(CLWDTH / DBLE(WEND-WBEGIN)) 00677 IF( SGNDEF.EQ.ONE ) THEN 00678 TAU = HALF*MAX(WGAP(WBEGIN),AVGAP) 00679 TAU = MAX(TAU,WERR(WBEGIN)) 00680 ELSE 00681 TAU = HALF*MAX(WGAP(WEND-1),AVGAP) 00682 TAU = MAX(TAU,WERR(WEND)) 00683 ENDIF 00684 ELSE 00685 TAU = WERR(WBEGIN) 00686 ENDIF 00687 ENDIF 00688 * 00689 DO 80 IDUM = 1, MAXTRY 00690 * Compute L D L^T factorization of tridiagonal matrix T - sigma I. 00691 * Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of 00692 * pivots in WORK(2*IN+1:3*IN) 00693 DPIVOT = D( IBEGIN ) - SIGMA 00694 WORK( 1 ) = DPIVOT 00695 DMAX = ABS( WORK(1) ) 00696 J = IBEGIN 00697 DO 70 I = 1, IN - 1 00698 WORK( 2*IN+I ) = ONE / WORK( I ) 00699 TMP = E( J )*WORK( 2*IN+I ) 00700 WORK( IN+I ) = TMP 00701 DPIVOT = ( D( J+1 )-SIGMA ) - TMP*E( J ) 00702 WORK( I+1 ) = DPIVOT 00703 DMAX = MAX( DMAX, ABS(DPIVOT) ) 00704 J = J + 1 00705 70 CONTINUE 00706 * check for element growth 00707 IF( DMAX .GT. MAXGROWTH*SPDIAM ) THEN 00708 NOREP = .TRUE. 00709 ELSE 00710 NOREP = .FALSE. 00711 ENDIF 00712 IF( USEDQD .AND. .NOT.NOREP ) THEN 00713 * Ensure the definiteness of the representation 00714 * All entries of D (of L D L^T) must have the same sign 00715 DO 71 I = 1, IN 00716 TMP = SGNDEF*WORK( I ) 00717 IF( TMP.LT.ZERO ) NOREP = .TRUE. 00718 71 CONTINUE 00719 ENDIF 00720 IF(NOREP) THEN 00721 * Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin 00722 * shift which makes the matrix definite. So we should end up 00723 * here really only in the case of IRANGE = VALRNG or INDRNG. 00724 IF( IDUM.EQ.MAXTRY-1 ) THEN 00725 IF( SGNDEF.EQ.ONE ) THEN 00726 * The fudged Gerschgorin shift should succeed 00727 SIGMA = 00728 $ GL - FUDGE*SPDIAM*EPS*N - FUDGE*TWO*PIVMIN 00729 ELSE 00730 SIGMA = 00731 $ GU + FUDGE*SPDIAM*EPS*N + FUDGE*TWO*PIVMIN 00732 END IF 00733 ELSE 00734 SIGMA = SIGMA - SGNDEF * TAU 00735 TAU = TWO * TAU 00736 END IF 00737 ELSE 00738 * an initial RRR is found 00739 GO TO 83 00740 END IF 00741 80 CONTINUE 00742 * if the program reaches this point, no base representation could be 00743 * found in MAXTRY iterations. 00744 INFO = 2 00745 RETURN 00746 00747 83 CONTINUE 00748 * At this point, we have found an initial base representation 00749 * T - SIGMA I = L D L^T with not too much element growth. 00750 * Store the shift. 00751 E( IEND ) = SIGMA 00752 * Store D and L. 00753 CALL DCOPY( IN, WORK, 1, D( IBEGIN ), 1 ) 00754 CALL DCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 ) 00755 00756 00757 IF(MB.GT.1 ) THEN 00758 * 00759 * Perturb each entry of the base representation by a small 00760 * (but random) relative amount to overcome difficulties with 00761 * glued matrices. 00762 * 00763 DO 122 I = 1, 4 00764 ISEED( I ) = 1 00765 122 CONTINUE 00766 00767 CALL DLARNV(2, ISEED, 2*IN-1, WORK(1)) 00768 DO 125 I = 1,IN-1 00769 D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(I)) 00770 E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(IN+I)) 00771 125 CONTINUE 00772 D(IEND) = D(IEND)*(ONE+EPS*FOUR*WORK(IN)) 00773 * 00774 ENDIF 00775 * 00776 * Don't update the Gerschgorin intervals because keeping track 00777 * of the updates would be too much work in DLARRV. 00778 * We update W instead and use it to locate the proper Gerschgorin 00779 * intervals. 00780 00781 * Compute the required eigenvalues of L D L' by bisection or dqds 00782 IF ( .NOT.USEDQD ) THEN 00783 * If DLARRD has been used, shift the eigenvalue approximations 00784 * according to their representation. This is necessary for 00785 * a uniform DLARRV since dqds computes eigenvalues of the 00786 * shifted representation. In DLARRV, W will always hold the 00787 * UNshifted eigenvalue approximation. 00788 DO 134 J=WBEGIN,WEND 00789 W(J) = W(J) - SIGMA 00790 WERR(J) = WERR(J) + ABS(W(J)) * EPS 00791 134 CONTINUE 00792 * call DLARRB to reduce eigenvalue error of the approximations 00793 * from DLARRD 00794 DO 135 I = IBEGIN, IEND-1 00795 WORK( I ) = D( I ) * E( I )**2 00796 135 CONTINUE 00797 * use bisection to find EV from INDL to INDU 00798 CALL DLARRB(IN, D(IBEGIN), WORK(IBEGIN), 00799 $ INDL, INDU, RTOL1, RTOL2, INDL-1, 00800 $ W(WBEGIN), WGAP(WBEGIN), WERR(WBEGIN), 00801 $ WORK( 2*N+1 ), IWORK, PIVMIN, SPDIAM, 00802 $ IN, IINFO ) 00803 IF( IINFO .NE. 0 ) THEN 00804 INFO = -4 00805 RETURN 00806 END IF 00807 * DLARRB computes all gaps correctly except for the last one 00808 * Record distance to VU/GU 00809 WGAP( WEND ) = MAX( ZERO, 00810 $ ( VU-SIGMA ) - ( W( WEND ) + WERR( WEND ) ) ) 00811 DO 138 I = INDL, INDU 00812 M = M + 1 00813 IBLOCK(M) = JBLK 00814 INDEXW(M) = I 00815 138 CONTINUE 00816 ELSE 00817 * Call dqds to get all eigs (and then possibly delete unwanted 00818 * eigenvalues). 00819 * Note that dqds finds the eigenvalues of the L D L^T representation 00820 * of T to high relative accuracy. High relative accuracy 00821 * might be lost when the shift of the RRR is subtracted to obtain 00822 * the eigenvalues of T. However, T is not guaranteed to define its 00823 * eigenvalues to high relative accuracy anyway. 00824 * Set RTOL to the order of the tolerance used in DLASQ2 00825 * This is an ESTIMATED error, the worst case bound is 4*N*EPS 00826 * which is usually too large and requires unnecessary work to be 00827 * done by bisection when computing the eigenvectors 00828 RTOL = LOG(DBLE(IN)) * FOUR * EPS 00829 J = IBEGIN 00830 DO 140 I = 1, IN - 1 00831 WORK( 2*I-1 ) = ABS( D( J ) ) 00832 WORK( 2*I ) = E( J )*E( J )*WORK( 2*I-1 ) 00833 J = J + 1 00834 140 CONTINUE 00835 WORK( 2*IN-1 ) = ABS( D( IEND ) ) 00836 WORK( 2*IN ) = ZERO 00837 CALL DLASQ2( IN, WORK, IINFO ) 00838 IF( IINFO .NE. 0 ) THEN 00839 * If IINFO = -5 then an index is part of a tight cluster 00840 * and should be changed. The index is in IWORK(1) and the 00841 * gap is in WORK(N+1) 00842 INFO = -5 00843 RETURN 00844 ELSE 00845 * Test that all eigenvalues are positive as expected 00846 DO 149 I = 1, IN 00847 IF( WORK( I ).LT.ZERO ) THEN 00848 INFO = -6 00849 RETURN 00850 ENDIF 00851 149 CONTINUE 00852 END IF 00853 IF( SGNDEF.GT.ZERO ) THEN 00854 DO 150 I = INDL, INDU 00855 M = M + 1 00856 W( M ) = WORK( IN-I+1 ) 00857 IBLOCK( M ) = JBLK 00858 INDEXW( M ) = I 00859 150 CONTINUE 00860 ELSE 00861 DO 160 I = INDL, INDU 00862 M = M + 1 00863 W( M ) = -WORK( I ) 00864 IBLOCK( M ) = JBLK 00865 INDEXW( M ) = I 00866 160 CONTINUE 00867 END IF 00868 00869 DO 165 I = M - MB + 1, M 00870 * the value of RTOL below should be the tolerance in DLASQ2 00871 WERR( I ) = RTOL * ABS( W(I) ) 00872 165 CONTINUE 00873 DO 166 I = M - MB + 1, M - 1 00874 * compute the right gap between the intervals 00875 WGAP( I ) = MAX( ZERO, 00876 $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) ) 00877 166 CONTINUE 00878 WGAP( M ) = MAX( ZERO, 00879 $ ( VU-SIGMA ) - ( W( M ) + WERR( M ) ) ) 00880 END IF 00881 * proceed with next block 00882 IBEGIN = IEND + 1 00883 WBEGIN = WEND + 1 00884 170 CONTINUE 00885 * 00886 00887 RETURN 00888 * 00889 * end of DLARRE 00890 * 00891 END