|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
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
1.7.3 
