 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLARRE( 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.1) -- 00008 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00009 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00010 * -- April 2011 -- 00011 * 00012 * .. Scalar Arguments .. 00013 CHARACTER RANGE 00014 INTEGER IL, INFO, IU, M, N, NSPLIT 00015 DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU 00016 * .. 00017 * .. Array Arguments .. 00018 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ), 00019 $ INDEXW( * ) 00020 DOUBLE PRECISION 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, DLARRE 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 * DSTEMR 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 DLASQ2) to 00037 * conpute all and then discard any unwanted one. 00038 * As an added benefit, DLARRE 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) DOUBLE PRECISION 00055 * VU (input/output) DOUBLE PRECISION 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', DLARRE 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 00089 * RTOL2 (input) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 ( DLARRE may use the 00115 * remaining N-M elements as workspace). 00116 * 00117 * WERR (output) DOUBLE PRECISION array, dimension (N) 00118 * The error bound on the corresponding eigenvalue in W. 00119 * 00120 * WGAP (output) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 00142 * The minimum pivot in the Sturm sequence for T. 00143 * 00144 * WORK (workspace) DOUBLE PRECISION 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 DLARRE. 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 DLARRD. 00158 * = 2: No base representation could be found in MAXTRY iterations. 00159 * Increasing MAXTRY and recompilation might be a remedy. 00160 * =-3: Problem in DLARRB when computing the refined root 00161 * representation for DLASQ2. 00162 * =-4: Problem in DLARRB when preforming bisection on the 00163 * desired part of the spectrum. 00164 * =-5: Problem in DLASQ2. 00165 * =-6: Problem in DLASQ2. 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 DOUBLE PRECISION FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD, 00184 $ MAXGROWTH, ONE, PERT, TWO, ZERO 00185 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, 00186 $ TWO = 2.0D0, FOUR=4.0D0, 00187 $ HNDRD = 100.0D0, 00188 $ PERT = 8.0D0, 00189 $ HALF = ONE/TWO, FOURTH = ONE/FOUR, FAC= HALF, 00190 $ MAXGROWTH = 64.0D0, FUDGE = 2.0D0 ) 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 DOUBLE PRECISION 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 DOUBLE PRECISION DLAMCH 00213 EXTERNAL DLAMCH, LSAME 00214 00215 * .. 00216 * .. External Subroutines .. 00217 EXTERNAL DCOPY, DLARNV, DLARRA, DLARRB, DLARRC, DLARRD, 00218 $ DLASQ2 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 = DLAMCH( 'S' ) 00244 EPS = DLAMCH( 'P' ) 00245 00246 * Set parameters 00247 RTL = SQRT(EPS) 00248 BSRTOL = SQRT(EPS) 00249 00250 * Treat case of 1x1 matrix for quick return 00251 IF( N.EQ.1 ) THEN 00252 IF( (IRANGE.EQ.ALLRNG).OR. 00253 $ ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR. 00254 $ ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN 00255 M = 1 00256 W(1) = D(1) 00257 * The computation error of the eigenvalue is zero 00258 WERR(1) = ZERO 00259 WGAP(1) = ZERO 00260 IBLOCK( 1 ) = 1 00261 INDEXW( 1 ) = 1 00262 GERS(1) = D( 1 ) 00263 GERS(2) = D( 1 ) 00264 ENDIF 00265 * store the shift for the initial RRR, which is zero in this case 00266 E(1) = ZERO 00267 RETURN 00268 END IF 00269 00270 * General case: tridiagonal matrix of order > 1 00271 * 00272 * Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter. 00273 * Compute maximum off-diagonal entry and pivmin. 00274 GL = D(1) 00275 GU = D(1) 00276 EOLD = ZERO 00277 EMAX = ZERO 00278 E(N) = ZERO 00279 DO 5 I = 1,N 00280 WERR(I) = ZERO 00281 WGAP(I) = ZERO 00282 EABS = ABS( E(I) ) 00283 IF( EABS .GE. EMAX ) THEN 00284 EMAX = EABS 00285 END IF 00286 TMP1 = EABS + EOLD 00287 GERS( 2*I-1) = D(I) - TMP1 00288 GL = MIN( GL, GERS( 2*I - 1)) 00289 GERS( 2*I ) = D(I) + TMP1 00290 GU = MAX( GU, GERS(2*I) ) 00291 EOLD = EABS 00292 5 CONTINUE 00293 * The minimum pivot allowed in the Sturm sequence for T 00294 PIVMIN = SAFMIN * MAX( ONE, EMAX**2 ) 00295 * Compute spectral diameter. The Gerschgorin bounds give an 00296 * estimate that is wrong by at most a factor of SQRT(2) 00297 SPDIAM = GU - GL 00298 00299 * Compute splitting points 00300 CALL DLARRA( N, D, E, E2, SPLTOL, SPDIAM, 00301 $ NSPLIT, ISPLIT, IINFO ) 00302 00303 * Can force use of bisection instead of faster DQDS. 00304 * Option left in the code for future multisection work. 00305 FORCEB = .FALSE. 00306 00307 * Initialize USEDQD, DQDS should be used for ALLRNG unless someone 00308 * explicitly wants bisection. 00309 USEDQD = (( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB)) 00310 00311 IF( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ) THEN 00312 * Set interval [VL,VU] that contains all eigenvalues 00313 VL = GL 00314 VU = GU 00315 ELSE 00316 * We call DLARRD to find crude approximations to the eigenvalues 00317 * in the desired range. In case IRANGE = INDRNG, we also obtain the 00318 * interval (VL,VU] that contains all the wanted eigenvalues. 00319 * An interval [LEFT,RIGHT] has converged if 00320 * RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT)) 00321 * DLARRD needs a WORK of size 4*N, IWORK of size 3*N 00322 CALL DLARRD( RANGE, 'B', N, VL, VU, IL, IU, GERS, 00323 $ BSRTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, 00324 $ MM, W, WERR, VL, VU, IBLOCK, INDEXW, 00325 $ WORK, IWORK, IINFO ) 00326 IF( IINFO.NE.0 ) THEN 00327 INFO = -1 00328 RETURN 00329 ENDIF 00330 * Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0 00331 DO 14 I = MM+1,N 00332 W( I ) = ZERO 00333 WERR( I ) = ZERO 00334 IBLOCK( I ) = 0 00335 INDEXW( I ) = 0 00336 14 CONTINUE 00337 END IF 00338 00339 00340 *** 00341 * Loop over unreduced blocks 00342 IBEGIN = 1 00343 WBEGIN = 1 00344 DO 170 JBLK = 1, NSPLIT 00345 IEND = ISPLIT( JBLK ) 00346 IN = IEND - IBEGIN + 1 00347 00348 * 1 X 1 block 00349 IF( IN.EQ.1 ) THEN 00350 IF( (IRANGE.EQ.ALLRNG).OR.( (IRANGE.EQ.VALRNG).AND. 00351 $ ( D( IBEGIN ).GT.VL ).AND.( D( IBEGIN ).LE.VU ) ) 00352 $ .OR. ( (IRANGE.EQ.INDRNG).AND.(IBLOCK(WBEGIN).EQ.JBLK)) 00353 $ ) THEN 00354 M = M + 1 00355 W( M ) = D( IBEGIN ) 00356 WERR(M) = ZERO 00357 * The gap for a single block doesn't matter for the later 00358 * algorithm and is assigned an arbitrary large value 00359 WGAP(M) = ZERO 00360 IBLOCK( M ) = JBLK 00361 INDEXW( M ) = 1 00362 WBEGIN = WBEGIN + 1 00363 ENDIF 00364 * E( IEND ) holds the shift for the initial RRR 00365 E( IEND ) = ZERO 00366 IBEGIN = IEND + 1 00367 GO TO 170 00368 END IF 00369 * 00370 * Blocks of size larger than 1x1 00371 * 00372 * E( IEND ) will hold the shift for the initial RRR, for now set it =0 00373 E( IEND ) = ZERO 00374 * 00375 * Find local outer bounds GL,GU for the block 00376 GL = D(IBEGIN) 00377 GU = D(IBEGIN) 00378 DO 15 I = IBEGIN , IEND 00379 GL = MIN( GERS( 2*I-1 ), GL ) 00380 GU = MAX( GERS( 2*I ), GU ) 00381 15 CONTINUE 00382 SPDIAM = GU - GL 00383 00384 IF(.NOT. ((IRANGE.EQ.ALLRNG).AND.(.NOT.FORCEB)) ) THEN 00385 * Count the number of eigenvalues in the current block. 00386 MB = 0 00387 DO 20 I = WBEGIN,MM 00388 IF( IBLOCK(I).EQ.JBLK ) THEN 00389 MB = MB+1 00390 ELSE 00391 GOTO 21 00392 ENDIF 00393 20 CONTINUE 00394 21 CONTINUE 00395 00396 IF( MB.EQ.0) THEN 00397 * No eigenvalue in the current block lies in the desired range 00398 * E( IEND ) holds the shift for the initial RRR 00399 E( IEND ) = ZERO 00400 IBEGIN = IEND + 1 00401 GO TO 170 00402 ELSE 00403 00404 * Decide whether dqds or bisection is more efficient 00405 USEDQD = ( (MB .GT. FAC*IN) .AND. (.NOT.FORCEB) ) 00406 WEND = WBEGIN + MB - 1 00407 * Calculate gaps for the current block 00408 * In later stages, when representations for individual 00409 * eigenvalues are different, we use SIGMA = E( IEND ). 00410 SIGMA = ZERO 00411 DO 30 I = WBEGIN, WEND - 1 00412 WGAP( I ) = MAX( ZERO, 00413 $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) ) 00414 30 CONTINUE 00415 WGAP( WEND ) = MAX( ZERO, 00416 $ VU - SIGMA - (W( WEND )+WERR( WEND ))) 00417 * Find local index of the first and last desired evalue. 00418 INDL = INDEXW(WBEGIN) 00419 INDU = INDEXW( WEND ) 00420 ENDIF 00421 ENDIF 00422 IF(( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ).OR.USEDQD) THEN 00423 * Case of DQDS 00424 * Find approximations to the extremal eigenvalues of the block 00425 CALL DLARRK( IN, 1, GL, GU, D(IBEGIN), 00426 $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO ) 00427 IF( IINFO.NE.0 ) THEN 00428 INFO = -1 00429 RETURN 00430 ENDIF 00431 ISLEFT = MAX(GL, TMP - TMP1 00432 $ - HNDRD * EPS* ABS(TMP - TMP1)) 00433 00434 CALL DLARRK( IN, IN, GL, GU, D(IBEGIN), 00435 $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO ) 00436 IF( IINFO.NE.0 ) THEN 00437 INFO = -1 00438 RETURN 00439 ENDIF 00440 ISRGHT = MIN(GU, TMP + TMP1 00441 $ + HNDRD * EPS * ABS(TMP + TMP1)) 00442 * Improve the estimate of the spectral diameter 00443 SPDIAM = ISRGHT - ISLEFT 00444 ELSE 00445 * Case of bisection 00446 * Find approximations to the wanted extremal eigenvalues 00447 ISLEFT = MAX(GL, W(WBEGIN) - WERR(WBEGIN) 00448 $ - HNDRD * EPS*ABS(W(WBEGIN)- WERR(WBEGIN) )) 00449 ISRGHT = MIN(GU,W(WEND) + WERR(WEND) 00450 $ + HNDRD * EPS * ABS(W(WEND)+ WERR(WEND))) 00451 ENDIF 00452 00453 00454 * Decide whether the base representation for the current block 00455 * L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I 00456 * should be on the left or the right end of the current block. 00457 * The strategy is to shift to the end which is "more populated" 00458 * Furthermore, decide whether to use DQDS for the computation of 00459 * the eigenvalue approximations at the end of DLARRE or bisection. 00460 * dqds is chosen if all eigenvalues are desired or the number of 00461 * eigenvalues to be computed is large compared to the blocksize. 00462 IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN 00463 * If all the eigenvalues have to be computed, we use dqd 00464 USEDQD = .TRUE. 00465 * INDL is the local index of the first eigenvalue to compute 00466 INDL = 1 00467 INDU = IN 00468 * MB = number of eigenvalues to compute 00469 MB = IN 00470 WEND = WBEGIN + MB - 1 00471 * Define 1/4 and 3/4 points of the spectrum 00472 S1 = ISLEFT + FOURTH * SPDIAM 00473 S2 = ISRGHT - FOURTH * SPDIAM 00474 ELSE 00475 * DLARRD has computed IBLOCK and INDEXW for each eigenvalue 00476 * approximation. 00477 * choose sigma 00478 IF( USEDQD ) THEN 00479 S1 = ISLEFT + FOURTH * SPDIAM 00480 S2 = ISRGHT - FOURTH * SPDIAM 00481 ELSE 00482 TMP = MIN(ISRGHT,VU) - MAX(ISLEFT,VL) 00483 S1 = MAX(ISLEFT,VL) + FOURTH * TMP 00484 S2 = MIN(ISRGHT,VU) - FOURTH * TMP 00485 ENDIF 00486 ENDIF 00487 00488 * Compute the negcount at the 1/4 and 3/4 points 00489 IF(MB.GT.1) THEN 00490 CALL DLARRC( 'T', IN, S1, S2, D(IBEGIN), 00491 $ E(IBEGIN), PIVMIN, CNT, CNT1, CNT2, IINFO) 00492 ENDIF 00493 00494 IF(MB.EQ.1) THEN 00495 SIGMA = GL 00496 SGNDEF = ONE 00497 ELSEIF( CNT1 - INDL .GE. INDU - CNT2 ) THEN 00498 IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN 00499 SIGMA = MAX(ISLEFT,GL) 00500 ELSEIF( USEDQD ) THEN 00501 * use Gerschgorin bound as shift to get pos def matrix 00502 * for dqds 00503 SIGMA = ISLEFT 00504 ELSE 00505 * use approximation of the first desired eigenvalue of the 00506 * block as shift 00507 SIGMA = MAX(ISLEFT,VL) 00508 ENDIF 00509 SGNDEF = ONE 00510 ELSE 00511 IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN 00512 SIGMA = MIN(ISRGHT,GU) 00513 ELSEIF( USEDQD ) THEN 00514 * use Gerschgorin bound as shift to get neg def matrix 00515 * for dqds 00516 SIGMA = ISRGHT 00517 ELSE 00518 * use approximation of the first desired eigenvalue of the 00519 * block as shift 00520 SIGMA = MIN(ISRGHT,VU) 00521 ENDIF 00522 SGNDEF = -ONE 00523 ENDIF 00524 00525 00526 * An initial SIGMA has been chosen that will be used for computing 00527 * T - SIGMA I = L D L^T 00528 * Define the increment TAU of the shift in case the initial shift 00529 * needs to be refined to obtain a factorization with not too much 00530 * element growth. 00531 IF( USEDQD ) THEN 00532 * The initial SIGMA was to the outer end of the spectrum 00533 * the matrix is definite and we need not retreat. 00534 TAU = SPDIAM*EPS*N + TWO*PIVMIN 00535 TAU = MAX( TAU,TWO*EPS*ABS(SIGMA) ) 00536 ELSE 00537 IF(MB.GT.1) THEN 00538 CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN) 00539 AVGAP = ABS(CLWDTH / DBLE(WEND-WBEGIN)) 00540 IF( SGNDEF.EQ.ONE ) THEN 00541 TAU = HALF*MAX(WGAP(WBEGIN),AVGAP) 00542 TAU = MAX(TAU,WERR(WBEGIN)) 00543 ELSE 00544 TAU = HALF*MAX(WGAP(WEND-1),AVGAP) 00545 TAU = MAX(TAU,WERR(WEND)) 00546 ENDIF 00547 ELSE 00548 TAU = WERR(WBEGIN) 00549 ENDIF 00550 ENDIF 00551 * 00552 DO 80 IDUM = 1, MAXTRY 00553 * Compute L D L^T factorization of tridiagonal matrix T - sigma I. 00554 * Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of 00555 * pivots in WORK(2*IN+1:3*IN) 00556 DPIVOT = D( IBEGIN ) - SIGMA 00557 WORK( 1 ) = DPIVOT 00558 DMAX = ABS( WORK(1) ) 00559 J = IBEGIN 00560 DO 70 I = 1, IN - 1 00561 WORK( 2*IN+I ) = ONE / WORK( I ) 00562 TMP = E( J )*WORK( 2*IN+I ) 00563 WORK( IN+I ) = TMP 00564 DPIVOT = ( D( J+1 )-SIGMA ) - TMP*E( J ) 00565 WORK( I+1 ) = DPIVOT 00566 DMAX = MAX( DMAX, ABS(DPIVOT) ) 00567 J = J + 1 00568 70 CONTINUE 00569 * check for element growth 00570 IF( DMAX .GT. MAXGROWTH*SPDIAM ) THEN 00571 NOREP = .TRUE. 00572 ELSE 00573 NOREP = .FALSE. 00574 ENDIF 00575 IF( USEDQD .AND. .NOT.NOREP ) THEN 00576 * Ensure the definiteness of the representation 00577 * All entries of D (of L D L^T) must have the same sign 00578 DO 71 I = 1, IN 00579 TMP = SGNDEF*WORK( I ) 00580 IF( TMP.LT.ZERO ) NOREP = .TRUE. 00581 71 CONTINUE 00582 ENDIF 00583 IF(NOREP) THEN 00584 * Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin 00585 * shift which makes the matrix definite. So we should end up 00586 * here really only in the case of IRANGE = VALRNG or INDRNG. 00587 IF( IDUM.EQ.MAXTRY-1 ) THEN 00588 IF( SGNDEF.EQ.ONE ) THEN 00589 * The fudged Gerschgorin shift should succeed 00590 SIGMA = 00591 $ GL - FUDGE*SPDIAM*EPS*N - FUDGE*TWO*PIVMIN 00592 ELSE 00593 SIGMA = 00594 $ GU + FUDGE*SPDIAM*EPS*N + FUDGE*TWO*PIVMIN 00595 END IF 00596 ELSE 00597 SIGMA = SIGMA - SGNDEF * TAU 00598 TAU = TWO * TAU 00599 END IF 00600 ELSE 00601 * an initial RRR is found 00602 GO TO 83 00603 END IF 00604 80 CONTINUE 00605 * if the program reaches this point, no base representation could be 00606 * found in MAXTRY iterations. 00607 INFO = 2 00608 RETURN 00609 00610 83 CONTINUE 00611 * At this point, we have found an initial base representation 00612 * T - SIGMA I = L D L^T with not too much element growth. 00613 * Store the shift. 00614 E( IEND ) = SIGMA 00615 * Store D and L. 00616 CALL DCOPY( IN, WORK, 1, D( IBEGIN ), 1 ) 00617 CALL DCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 ) 00618 00619 00620 IF(MB.GT.1 ) THEN 00621 * 00622 * Perturb each entry of the base representation by a small 00623 * (but random) relative amount to overcome difficulties with 00624 * glued matrices. 00625 * 00626 DO 122 I = 1, 4 00627 ISEED( I ) = 1 00628 122 CONTINUE 00629 00630 CALL DLARNV(2, ISEED, 2*IN-1, WORK(1)) 00631 DO 125 I = 1,IN-1 00632 D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(I)) 00633 E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(IN+I)) 00634 125 CONTINUE 00635 D(IEND) = D(IEND)*(ONE+EPS*FOUR*WORK(IN)) 00636 * 00637 ENDIF 00638 * 00639 * Don't update the Gerschgorin intervals because keeping track 00640 * of the updates would be too much work in DLARRV. 00641 * We update W instead and use it to locate the proper Gerschgorin 00642 * intervals. 00643 00644 * Compute the required eigenvalues of L D L' by bisection or dqds 00645 IF ( .NOT.USEDQD ) THEN 00646 * If DLARRD has been used, shift the eigenvalue approximations 00647 * according to their representation. This is necessary for 00648 * a uniform DLARRV since dqds computes eigenvalues of the 00649 * shifted representation. In DLARRV, W will always hold the 00650 * UNshifted eigenvalue approximation. 00651 DO 134 J=WBEGIN,WEND 00652 W(J) = W(J) - SIGMA 00653 WERR(J) = WERR(J) + ABS(W(J)) * EPS 00654 134 CONTINUE 00655 * call DLARRB to reduce eigenvalue error of the approximations 00656 * from DLARRD 00657 DO 135 I = IBEGIN, IEND-1 00658 WORK( I ) = D( I ) * E( I )**2 00659 135 CONTINUE 00660 * use bisection to find EV from INDL to INDU 00661 CALL DLARRB(IN, D(IBEGIN), WORK(IBEGIN), 00662 $ INDL, INDU, RTOL1, RTOL2, INDL-1, 00663 $ W(WBEGIN), WGAP(WBEGIN), WERR(WBEGIN), 00664 $ WORK( 2*N+1 ), IWORK, PIVMIN, SPDIAM, 00665 $ IN, IINFO ) 00666 IF( IINFO .NE. 0 ) THEN 00667 INFO = -4 00668 RETURN 00669 END IF 00670 * DLARRB computes all gaps correctly except for the last one 00671 * Record distance to VU/GU 00672 WGAP( WEND ) = MAX( ZERO, 00673 $ ( VU-SIGMA ) - ( W( WEND ) + WERR( WEND ) ) ) 00674 DO 138 I = INDL, INDU 00675 M = M + 1 00676 IBLOCK(M) = JBLK 00677 INDEXW(M) = I 00678 138 CONTINUE 00679 ELSE 00680 * Call dqds to get all eigs (and then possibly delete unwanted 00681 * eigenvalues). 00682 * Note that dqds finds the eigenvalues of the L D L^T representation 00683 * of T to high relative accuracy. High relative accuracy 00684 * might be lost when the shift of the RRR is subtracted to obtain 00685 * the eigenvalues of T. However, T is not guaranteed to define its 00686 * eigenvalues to high relative accuracy anyway. 00687 * Set RTOL to the order of the tolerance used in DLASQ2 00688 * This is an ESTIMATED error, the worst case bound is 4*N*EPS 00689 * which is usually too large and requires unnecessary work to be 00690 * done by bisection when computing the eigenvectors 00691 RTOL = LOG(DBLE(IN)) * FOUR * EPS 00692 J = IBEGIN 00693 DO 140 I = 1, IN - 1 00694 WORK( 2*I-1 ) = ABS( D( J ) ) 00695 WORK( 2*I ) = E( J )*E( J )*WORK( 2*I-1 ) 00696 J = J + 1 00697 140 CONTINUE 00698 WORK( 2*IN-1 ) = ABS( D( IEND ) ) 00699 WORK( 2*IN ) = ZERO 00700 CALL DLASQ2( IN, WORK, IINFO ) 00701 IF( IINFO .NE. 0 ) THEN 00702 * If IINFO = -5 then an index is part of a tight cluster 00703 * and should be changed. The index is in IWORK(1) and the 00704 * gap is in WORK(N+1) 00705 INFO = -5 00706 RETURN 00707 ELSE 00708 * Test that all eigenvalues are positive as expected 00709 DO 149 I = 1, IN 00710 IF( WORK( I ).LT.ZERO ) THEN 00711 INFO = -6 00712 RETURN 00713 ENDIF 00714 149 CONTINUE 00715 END IF 00716 IF( SGNDEF.GT.ZERO ) THEN 00717 DO 150 I = INDL, INDU 00718 M = M + 1 00719 W( M ) = WORK( IN-I+1 ) 00720 IBLOCK( M ) = JBLK 00721 INDEXW( M ) = I 00722 150 CONTINUE 00723 ELSE 00724 DO 160 I = INDL, INDU 00725 M = M + 1 00726 W( M ) = -WORK( I ) 00727 IBLOCK( M ) = JBLK 00728 INDEXW( M ) = I 00729 160 CONTINUE 00730 END IF 00731 00732 DO 165 I = M - MB + 1, M 00733 * the value of RTOL below should be the tolerance in DLASQ2 00734 WERR( I ) = RTOL * ABS( W(I) ) 00735 165 CONTINUE 00736 DO 166 I = M - MB + 1, M - 1 00737 * compute the right gap between the intervals 00738 WGAP( I ) = MAX( ZERO, 00739 $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) ) 00740 166 CONTINUE 00741 WGAP( M ) = MAX( ZERO, 00742 $ ( VU-SIGMA ) - ( W( M ) + WERR( M ) ) ) 00743 END IF 00744 * proceed with next block 00745 IBEGIN = IEND + 1 00746 WBEGIN = WEND + 1 00747 170 CONTINUE 00748 * 00749 00750 RETURN 00751 * 00752 * end of DLARRE 00753 * 00754 END
1.7.3 
