|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DLARRD( RANGE, ORDER, N, VL, VU, IL, IU, GERS, 00002 $ RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, 00003 $ M, W, WERR, WL, WU, IBLOCK, INDEXW, 00004 $ WORK, IWORK, INFO ) 00005 * 00006 * -- LAPACK auxiliary routine (version 3.3.0) -- 00007 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00008 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00009 * November 2010 00010 * 00011 * .. Scalar Arguments .. 00012 CHARACTER ORDER, RANGE 00013 INTEGER IL, INFO, IU, M, N, NSPLIT 00014 DOUBLE PRECISION PIVMIN, RELTOL, VL, VU, WL, WU 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER IBLOCK( * ), INDEXW( * ), 00018 $ ISPLIT( * ), IWORK( * ) 00019 DOUBLE PRECISION D( * ), E( * ), E2( * ), 00020 $ GERS( * ), W( * ), WERR( * ), WORK( * ) 00021 * .. 00022 * 00023 * Purpose 00024 * ======= 00025 * 00026 * DLARRD computes the eigenvalues of a symmetric tridiagonal 00027 * matrix T to suitable accuracy. This is an auxiliary code to be 00028 * called from DSTEMR. 00029 * The user may ask for all eigenvalues, all eigenvalues 00030 * in the half-open interval (VL, VU], or the IL-th through IU-th 00031 * eigenvalues. 00032 * 00033 * To avoid overflow, the matrix must be scaled so that its 00034 * largest element is no greater than overflow**(1/2) * 00035 * underflow**(1/4) in absolute value, and for greatest 00036 * accuracy, it should not be much smaller than that. 00037 * 00038 * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal 00039 * Matrix", Report CS41, Computer Science Dept., Stanford 00040 * University, July 21, 1966. 00041 * 00042 * Arguments 00043 * ========= 00044 * 00045 * RANGE (input) CHARACTER*1 00046 * = 'A': ("All") all eigenvalues will be found. 00047 * = 'V': ("Value") all eigenvalues in the half-open interval 00048 * (VL, VU] will be found. 00049 * = 'I': ("Index") the IL-th through IU-th eigenvalues (of the 00050 * entire matrix) will be found. 00051 * 00052 * ORDER (input) CHARACTER*1 00053 * = 'B': ("By Block") the eigenvalues will be grouped by 00054 * split-off block (see IBLOCK, ISPLIT) and 00055 * ordered from smallest to largest within 00056 * the block. 00057 * = 'E': ("Entire matrix") 00058 * the eigenvalues for the entire matrix 00059 * will be ordered from smallest to 00060 * largest. 00061 * 00062 * N (input) INTEGER 00063 * The order of the tridiagonal matrix T. N >= 0. 00064 * 00065 * VL (input) DOUBLE PRECISION 00066 * VU (input) DOUBLE PRECISION 00067 * If RANGE='V', the lower and upper bounds of the interval to 00068 * be searched for eigenvalues. Eigenvalues less than or equal 00069 * to VL, or greater than VU, will not be returned. VL < VU. 00070 * Not referenced if RANGE = 'A' or 'I'. 00071 * 00072 * IL (input) INTEGER 00073 * IU (input) INTEGER 00074 * If RANGE='I', the indices (in ascending order) of the 00075 * smallest and largest eigenvalues to be returned. 00076 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 00077 * Not referenced if RANGE = 'A' or 'V'. 00078 * 00079 * GERS (input) DOUBLE PRECISION array, dimension (2*N) 00080 * The N Gerschgorin intervals (the i-th Gerschgorin interval 00081 * is (GERS(2*i-1), GERS(2*i)). 00082 * 00083 * RELTOL (input) DOUBLE PRECISION 00084 * The minimum relative width of an interval. When an interval 00085 * is narrower than RELTOL times the larger (in 00086 * magnitude) endpoint, then it is considered to be 00087 * sufficiently small, i.e., converged. Note: this should 00088 * always be at least radix*machine epsilon. 00089 * 00090 * D (input) DOUBLE PRECISION array, dimension (N) 00091 * The n diagonal elements of the tridiagonal matrix T. 00092 * 00093 * E (input) DOUBLE PRECISION array, dimension (N-1) 00094 * The (n-1) off-diagonal elements of the tridiagonal matrix T. 00095 * 00096 * E2 (input) DOUBLE PRECISION array, dimension (N-1) 00097 * The (n-1) squared off-diagonal elements of the tridiagonal matrix T. 00098 * 00099 * PIVMIN (input) DOUBLE PRECISION 00100 * The minimum pivot allowed in the Sturm sequence for T. 00101 * 00102 * NSPLIT (input) INTEGER 00103 * The number of diagonal blocks in the matrix T. 00104 * 1 <= NSPLIT <= N. 00105 * 00106 * ISPLIT (input) INTEGER array, dimension (N) 00107 * The splitting points, at which T breaks up into submatrices. 00108 * The first submatrix consists of rows/columns 1 to ISPLIT(1), 00109 * the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), 00110 * etc., and the NSPLIT-th consists of rows/columns 00111 * ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. 00112 * (Only the first NSPLIT elements will actually be used, but 00113 * since the user cannot know a priori what value NSPLIT will 00114 * have, N words must be reserved for ISPLIT.) 00115 * 00116 * M (output) INTEGER 00117 * The actual number of eigenvalues found. 0 <= M <= N. 00118 * (See also the description of INFO=2,3.) 00119 * 00120 * W (output) DOUBLE PRECISION array, dimension (N) 00121 * On exit, the first M elements of W will contain the 00122 * eigenvalue approximations. DLARRD computes an interval 00123 * I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue 00124 * approximation is given as the interval midpoint 00125 * W(j)= ( a_j + b_j)/2. The corresponding error is bounded by 00126 * WERR(j) = abs( a_j - b_j)/2 00127 * 00128 * WERR (output) DOUBLE PRECISION array, dimension (N) 00129 * The error bound on the corresponding eigenvalue approximation 00130 * in W. 00131 * 00132 * WL (output) DOUBLE PRECISION 00133 * WU (output) DOUBLE PRECISION 00134 * The interval (WL, WU] contains all the wanted eigenvalues. 00135 * If RANGE='V', then WL=VL and WU=VU. 00136 * If RANGE='A', then WL and WU are the global Gerschgorin bounds 00137 * on the spectrum. 00138 * If RANGE='I', then WL and WU are computed by DLAEBZ from the 00139 * index range specified. 00140 * 00141 * IBLOCK (output) INTEGER array, dimension (N) 00142 * At each row/column j where E(j) is zero or small, the 00143 * matrix T is considered to split into a block diagonal 00144 * matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which 00145 * block (from 1 to the number of blocks) the eigenvalue W(i) 00146 * belongs. (DLARRD may use the remaining N-M elements as 00147 * workspace.) 00148 * 00149 * INDEXW (output) INTEGER array, dimension (N) 00150 * The indices of the eigenvalues within each block (submatrix); 00151 * for example, INDEXW(i)= j and IBLOCK(i)=k imply that the 00152 * i-th eigenvalue W(i) is the j-th eigenvalue in block k. 00153 * 00154 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N) 00155 * 00156 * IWORK (workspace) INTEGER array, dimension (3*N) 00157 * 00158 * INFO (output) INTEGER 00159 * = 0: successful exit 00160 * < 0: if INFO = -i, the i-th argument had an illegal value 00161 * > 0: some or all of the eigenvalues failed to converge or 00162 * were not computed: 00163 * =1 or 3: Bisection failed to converge for some 00164 * eigenvalues; these eigenvalues are flagged by a 00165 * negative block number. The effect is that the 00166 * eigenvalues may not be as accurate as the 00167 * absolute and relative tolerances. This is 00168 * generally caused by unexpectedly inaccurate 00169 * arithmetic. 00170 * =2 or 3: RANGE='I' only: Not all of the eigenvalues 00171 * IL:IU were found. 00172 * Effect: M < IU+1-IL 00173 * Cause: non-monotonic arithmetic, causing the 00174 * Sturm sequence to be non-monotonic. 00175 * Cure: recalculate, using RANGE='A', and pick 00176 * out eigenvalues IL:IU. In some cases, 00177 * increasing the PARAMETER "FUDGE" may 00178 * make things work. 00179 * = 4: RANGE='I', and the Gershgorin interval 00180 * initially used was too small. No eigenvalues 00181 * were computed. 00182 * Probable cause: your machine has sloppy 00183 * floating-point arithmetic. 00184 * Cure: Increase the PARAMETER "FUDGE", 00185 * recompile, and try again. 00186 * 00187 * Internal Parameters 00188 * =================== 00189 * 00190 * FUDGE DOUBLE PRECISION, default = 2 00191 * A "fudge factor" to widen the Gershgorin intervals. Ideally, 00192 * a value of 1 should work, but on machines with sloppy 00193 * arithmetic, this needs to be larger. The default for 00194 * publicly released versions should be large enough to handle 00195 * the worst machine around. Note that this has no effect 00196 * on accuracy of the solution. 00197 * 00198 * Based on contributions by 00199 * W. Kahan, University of California, Berkeley, USA 00200 * Beresford Parlett, University of California, Berkeley, USA 00201 * Jim Demmel, University of California, Berkeley, USA 00202 * Inderjit Dhillon, University of Texas, Austin, USA 00203 * Osni Marques, LBNL/NERSC, USA 00204 * Christof Voemel, University of California, Berkeley, USA 00205 * 00206 * ===================================================================== 00207 * 00208 * .. Parameters .. 00209 DOUBLE PRECISION ZERO, ONE, TWO, HALF, FUDGE 00210 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, 00211 $ TWO = 2.0D0, HALF = ONE/TWO, 00212 $ FUDGE = TWO ) 00213 INTEGER ALLRNG, VALRNG, INDRNG 00214 PARAMETER ( ALLRNG = 1, VALRNG = 2, INDRNG = 3 ) 00215 * .. 00216 * .. Local Scalars .. 00217 LOGICAL NCNVRG, TOOFEW 00218 INTEGER I, IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO, 00219 $ IM, IN, IOFF, IOUT, IRANGE, ITMAX, ITMP1, 00220 $ ITMP2, IW, IWOFF, J, JBLK, JDISC, JE, JEE, NB, 00221 $ NWL, NWU 00222 DOUBLE PRECISION ATOLI, EPS, GL, GU, RTOLI, TMP1, TMP2, 00223 $ TNORM, UFLOW, WKILL, WLU, WUL 00224 00225 * .. 00226 * .. Local Arrays .. 00227 INTEGER IDUMMA( 1 ) 00228 * .. 00229 * .. External Functions .. 00230 LOGICAL LSAME 00231 INTEGER ILAENV 00232 DOUBLE PRECISION DLAMCH 00233 EXTERNAL LSAME, ILAENV, DLAMCH 00234 * .. 00235 * .. External Subroutines .. 00236 EXTERNAL DLAEBZ 00237 * .. 00238 * .. Intrinsic Functions .. 00239 INTRINSIC ABS, INT, LOG, MAX, MIN 00240 * .. 00241 * .. Executable Statements .. 00242 * 00243 INFO = 0 00244 * 00245 * Decode RANGE 00246 * 00247 IF( LSAME( RANGE, 'A' ) ) THEN 00248 IRANGE = ALLRNG 00249 ELSE IF( LSAME( RANGE, 'V' ) ) THEN 00250 IRANGE = VALRNG 00251 ELSE IF( LSAME( RANGE, 'I' ) ) THEN 00252 IRANGE = INDRNG 00253 ELSE 00254 IRANGE = 0 00255 END IF 00256 * 00257 * Check for Errors 00258 * 00259 IF( IRANGE.LE.0 ) THEN 00260 INFO = -1 00261 ELSE IF( .NOT.(LSAME(ORDER,'B').OR.LSAME(ORDER,'E')) ) THEN 00262 INFO = -2 00263 ELSE IF( N.LT.0 ) THEN 00264 INFO = -3 00265 ELSE IF( IRANGE.EQ.VALRNG ) THEN 00266 IF( VL.GE.VU ) 00267 $ INFO = -5 00268 ELSE IF( IRANGE.EQ.INDRNG .AND. 00269 $ ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) THEN 00270 INFO = -6 00271 ELSE IF( IRANGE.EQ.INDRNG .AND. 00272 $ ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN 00273 INFO = -7 00274 END IF 00275 * 00276 IF( INFO.NE.0 ) THEN 00277 RETURN 00278 END IF 00279 00280 * Initialize error flags 00281 INFO = 0 00282 NCNVRG = .FALSE. 00283 TOOFEW = .FALSE. 00284 00285 * Quick return if possible 00286 M = 0 00287 IF( N.EQ.0 ) RETURN 00288 00289 * Simplification: 00290 IF( IRANGE.EQ.INDRNG .AND. IL.EQ.1 .AND. IU.EQ.N ) IRANGE = 1 00291 00292 * Get machine constants 00293 EPS = DLAMCH( 'P' ) 00294 UFLOW = DLAMCH( 'U' ) 00295 00296 00297 * Special Case when N=1 00298 * Treat case of 1x1 matrix for quick return 00299 IF( N.EQ.1 ) THEN 00300 IF( (IRANGE.EQ.ALLRNG).OR. 00301 $ ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR. 00302 $ ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN 00303 M = 1 00304 W(1) = D(1) 00305 * The computation error of the eigenvalue is zero 00306 WERR(1) = ZERO 00307 IBLOCK( 1 ) = 1 00308 INDEXW( 1 ) = 1 00309 ENDIF 00310 RETURN 00311 END IF 00312 00313 * NB is the minimum vector length for vector bisection, or 0 00314 * if only scalar is to be done. 00315 NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 ) 00316 IF( NB.LE.1 ) NB = 0 00317 00318 * Find global spectral radius 00319 GL = D(1) 00320 GU = D(1) 00321 DO 5 I = 1,N 00322 GL = MIN( GL, GERS( 2*I - 1)) 00323 GU = MAX( GU, GERS(2*I) ) 00324 5 CONTINUE 00325 * Compute global Gerschgorin bounds and spectral diameter 00326 TNORM = MAX( ABS( GL ), ABS( GU ) ) 00327 GL = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN 00328 GU = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN 00329 * [JAN/28/2009] remove the line below since SPDIAM variable not use 00330 * SPDIAM = GU - GL 00331 * Input arguments for DLAEBZ: 00332 * The relative tolerance. An interval (a,b] lies within 00333 * "relative tolerance" if b-a < RELTOL*max(|a|,|b|), 00334 RTOLI = RELTOL 00335 * Set the absolute tolerance for interval convergence to zero to force 00336 * interval convergence based on relative size of the interval. 00337 * This is dangerous because intervals might not converge when RELTOL is 00338 * small. But at least a very small number should be selected so that for 00339 * strongly graded matrices, the code can get relatively accurate 00340 * eigenvalues. 00341 ATOLI = FUDGE*TWO*UFLOW + FUDGE*TWO*PIVMIN 00342 00343 IF( IRANGE.EQ.INDRNG ) THEN 00344 00345 * RANGE='I': Compute an interval containing eigenvalues 00346 * IL through IU. The initial interval [GL,GU] from the global 00347 * Gerschgorin bounds GL and GU is refined by DLAEBZ. 00348 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) / 00349 $ LOG( TWO ) ) + 2 00350 WORK( N+1 ) = GL 00351 WORK( N+2 ) = GL 00352 WORK( N+3 ) = GU 00353 WORK( N+4 ) = GU 00354 WORK( N+5 ) = GL 00355 WORK( N+6 ) = GU 00356 IWORK( 1 ) = -1 00357 IWORK( 2 ) = -1 00358 IWORK( 3 ) = N + 1 00359 IWORK( 4 ) = N + 1 00360 IWORK( 5 ) = IL - 1 00361 IWORK( 6 ) = IU 00362 * 00363 CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, 00364 $ D, E, E2, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT, 00365 $ IWORK, W, IBLOCK, IINFO ) 00366 IF( IINFO .NE. 0 ) THEN 00367 INFO = IINFO 00368 RETURN 00369 END IF 00370 * On exit, output intervals may not be ordered by ascending negcount 00371 IF( IWORK( 6 ).EQ.IU ) THEN 00372 WL = WORK( N+1 ) 00373 WLU = WORK( N+3 ) 00374 NWL = IWORK( 1 ) 00375 WU = WORK( N+4 ) 00376 WUL = WORK( N+2 ) 00377 NWU = IWORK( 4 ) 00378 ELSE 00379 WL = WORK( N+2 ) 00380 WLU = WORK( N+4 ) 00381 NWL = IWORK( 2 ) 00382 WU = WORK( N+3 ) 00383 WUL = WORK( N+1 ) 00384 NWU = IWORK( 3 ) 00385 END IF 00386 * On exit, the interval [WL, WLU] contains a value with negcount NWL, 00387 * and [WUL, WU] contains a value with negcount NWU. 00388 IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN 00389 INFO = 4 00390 RETURN 00391 END IF 00392 00393 ELSEIF( IRANGE.EQ.VALRNG ) THEN 00394 WL = VL 00395 WU = VU 00396 00397 ELSEIF( IRANGE.EQ.ALLRNG ) THEN 00398 WL = GL 00399 WU = GU 00400 ENDIF 00401 00402 00403 00404 * Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU. 00405 * NWL accumulates the number of eigenvalues .le. WL, 00406 * NWU accumulates the number of eigenvalues .le. WU 00407 M = 0 00408 IEND = 0 00409 INFO = 0 00410 NWL = 0 00411 NWU = 0 00412 * 00413 DO 70 JBLK = 1, NSPLIT 00414 IOFF = IEND 00415 IBEGIN = IOFF + 1 00416 IEND = ISPLIT( JBLK ) 00417 IN = IEND - IOFF 00418 * 00419 IF( IN.EQ.1 ) THEN 00420 * 1x1 block 00421 IF( WL.GE.D( IBEGIN )-PIVMIN ) 00422 $ NWL = NWL + 1 00423 IF( WU.GE.D( IBEGIN )-PIVMIN ) 00424 $ NWU = NWU + 1 00425 IF( IRANGE.EQ.ALLRNG .OR. 00426 $ ( WL.LT.D( IBEGIN )-PIVMIN 00427 $ .AND. WU.GE. D( IBEGIN )-PIVMIN ) ) THEN 00428 M = M + 1 00429 W( M ) = D( IBEGIN ) 00430 WERR(M) = ZERO 00431 * The gap for a single block doesn't matter for the later 00432 * algorithm and is assigned an arbitrary large value 00433 IBLOCK( M ) = JBLK 00434 INDEXW( M ) = 1 00435 END IF 00436 00437 * Disabled 2x2 case because of a failure on the following matrix 00438 * RANGE = 'I', IL = IU = 4 00439 * Original Tridiagonal, d = [ 00440 * -0.150102010615740E+00 00441 * -0.849897989384260E+00 00442 * -0.128208148052635E-15 00443 * 0.128257718286320E-15 00444 * ]; 00445 * e = [ 00446 * -0.357171383266986E+00 00447 * -0.180411241501588E-15 00448 * -0.175152352710251E-15 00449 * ]; 00450 * 00451 * ELSE IF( IN.EQ.2 ) THEN 00452 ** 2x2 block 00453 * DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 ) 00454 * TMP1 = HALF*(D(IBEGIN)+D(IEND)) 00455 * L1 = TMP1 - DISC 00456 * IF( WL.GE. L1-PIVMIN ) 00457 * $ NWL = NWL + 1 00458 * IF( WU.GE. L1-PIVMIN ) 00459 * $ NWU = NWU + 1 00460 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE. 00461 * $ L1-PIVMIN ) ) THEN 00462 * M = M + 1 00463 * W( M ) = L1 00464 ** The uncertainty of eigenvalues of a 2x2 matrix is very small 00465 * WERR( M ) = EPS * ABS( W( M ) ) * TWO 00466 * IBLOCK( M ) = JBLK 00467 * INDEXW( M ) = 1 00468 * ENDIF 00469 * L2 = TMP1 + DISC 00470 * IF( WL.GE. L2-PIVMIN ) 00471 * $ NWL = NWL + 1 00472 * IF( WU.GE. L2-PIVMIN ) 00473 * $ NWU = NWU + 1 00474 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE. 00475 * $ L2-PIVMIN ) ) THEN 00476 * M = M + 1 00477 * W( M ) = L2 00478 ** The uncertainty of eigenvalues of a 2x2 matrix is very small 00479 * WERR( M ) = EPS * ABS( W( M ) ) * TWO 00480 * IBLOCK( M ) = JBLK 00481 * INDEXW( M ) = 2 00482 * ENDIF 00483 ELSE 00484 * General Case - block of size IN >= 2 00485 * Compute local Gerschgorin interval and use it as the initial 00486 * interval for DLAEBZ 00487 GU = D( IBEGIN ) 00488 GL = D( IBEGIN ) 00489 TMP1 = ZERO 00490 00491 DO 40 J = IBEGIN, IEND 00492 GL = MIN( GL, GERS( 2*J - 1)) 00493 GU = MAX( GU, GERS(2*J) ) 00494 40 CONTINUE 00495 * [JAN/28/2009] 00496 * change SPDIAM by TNORM in lines 2 and 3 thereafter 00497 * line 1: remove computation of SPDIAM (not useful anymore) 00498 * SPDIAM = GU - GL 00499 * GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN 00500 * GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN 00501 GL = GL - FUDGE*TNORM*EPS*IN - FUDGE*PIVMIN 00502 GU = GU + FUDGE*TNORM*EPS*IN + FUDGE*PIVMIN 00503 * 00504 IF( IRANGE.GT.1 ) THEN 00505 IF( GU.LT.WL ) THEN 00506 * the local block contains none of the wanted eigenvalues 00507 NWL = NWL + IN 00508 NWU = NWU + IN 00509 GO TO 70 00510 END IF 00511 * refine search interval if possible, only range (WL,WU] matters 00512 GL = MAX( GL, WL ) 00513 GU = MIN( GU, WU ) 00514 IF( GL.GE.GU ) 00515 $ GO TO 70 00516 END IF 00517 00518 * Find negcount of initial interval boundaries GL and GU 00519 WORK( N+1 ) = GL 00520 WORK( N+IN+1 ) = GU 00521 CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 00522 $ D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ), 00523 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM, 00524 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 00525 IF( IINFO .NE. 0 ) THEN 00526 INFO = IINFO 00527 RETURN 00528 END IF 00529 * 00530 NWL = NWL + IWORK( 1 ) 00531 NWU = NWU + IWORK( IN+1 ) 00532 IWOFF = M - IWORK( 1 ) 00533 00534 * Compute Eigenvalues 00535 ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) / 00536 $ LOG( TWO ) ) + 2 00537 CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 00538 $ D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ), 00539 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT, 00540 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 00541 IF( IINFO .NE. 0 ) THEN 00542 INFO = IINFO 00543 RETURN 00544 END IF 00545 * 00546 * Copy eigenvalues into W and IBLOCK 00547 * Use -JBLK for block number for unconverged eigenvalues. 00548 * Loop over the number of output intervals from DLAEBZ 00549 DO 60 J = 1, IOUT 00550 * eigenvalue approximation is middle point of interval 00551 TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) ) 00552 * semi length of error interval 00553 TMP2 = HALF*ABS( WORK( J+N )-WORK( J+IN+N ) ) 00554 IF( J.GT.IOUT-IINFO ) THEN 00555 * Flag non-convergence. 00556 NCNVRG = .TRUE. 00557 IB = -JBLK 00558 ELSE 00559 IB = JBLK 00560 END IF 00561 DO 50 JE = IWORK( J ) + 1 + IWOFF, 00562 $ IWORK( J+IN ) + IWOFF 00563 W( JE ) = TMP1 00564 WERR( JE ) = TMP2 00565 INDEXW( JE ) = JE - IWOFF 00566 IBLOCK( JE ) = IB 00567 50 CONTINUE 00568 60 CONTINUE 00569 * 00570 M = M + IM 00571 END IF 00572 70 CONTINUE 00573 00574 * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU 00575 * If NWL+1 < IL or NWU > IU, discard extra eigenvalues. 00576 IF( IRANGE.EQ.INDRNG ) THEN 00577 IDISCL = IL - 1 - NWL 00578 IDISCU = NWU - IU 00579 * 00580 IF( IDISCL.GT.0 ) THEN 00581 IM = 0 00582 DO 80 JE = 1, M 00583 * Remove some of the smallest eigenvalues from the left so that 00584 * at the end IDISCL =0. Move all eigenvalues up to the left. 00585 IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN 00586 IDISCL = IDISCL - 1 00587 ELSE 00588 IM = IM + 1 00589 W( IM ) = W( JE ) 00590 WERR( IM ) = WERR( JE ) 00591 INDEXW( IM ) = INDEXW( JE ) 00592 IBLOCK( IM ) = IBLOCK( JE ) 00593 END IF 00594 80 CONTINUE 00595 M = IM 00596 END IF 00597 IF( IDISCU.GT.0 ) THEN 00598 * Remove some of the largest eigenvalues from the right so that 00599 * at the end IDISCU =0. Move all eigenvalues up to the left. 00600 IM=M+1 00601 DO 81 JE = M, 1, -1 00602 IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN 00603 IDISCU = IDISCU - 1 00604 ELSE 00605 IM = IM - 1 00606 W( IM ) = W( JE ) 00607 WERR( IM ) = WERR( JE ) 00608 INDEXW( IM ) = INDEXW( JE ) 00609 IBLOCK( IM ) = IBLOCK( JE ) 00610 END IF 00611 81 CONTINUE 00612 JEE = 0 00613 DO 82 JE = IM, M 00614 JEE = JEE + 1 00615 W( JEE ) = W( JE ) 00616 WERR( JEE ) = WERR( JE ) 00617 INDEXW( JEE ) = INDEXW( JE ) 00618 IBLOCK( JEE ) = IBLOCK( JE ) 00619 82 CONTINUE 00620 M = M-IM+1 00621 END IF 00622 00623 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN 00624 * Code to deal with effects of bad arithmetic. (If N(w) is 00625 * monotone non-decreasing, this should never happen.) 00626 * Some low eigenvalues to be discarded are not in (WL,WLU], 00627 * or high eigenvalues to be discarded are not in (WUL,WU] 00628 * so just kill off the smallest IDISCL/largest IDISCU 00629 * eigenvalues, by marking the corresponding IBLOCK = 0 00630 IF( IDISCL.GT.0 ) THEN 00631 WKILL = WU 00632 DO 100 JDISC = 1, IDISCL 00633 IW = 0 00634 DO 90 JE = 1, M 00635 IF( IBLOCK( JE ).NE.0 .AND. 00636 $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN 00637 IW = JE 00638 WKILL = W( JE ) 00639 END IF 00640 90 CONTINUE 00641 IBLOCK( IW ) = 0 00642 100 CONTINUE 00643 END IF 00644 IF( IDISCU.GT.0 ) THEN 00645 WKILL = WL 00646 DO 120 JDISC = 1, IDISCU 00647 IW = 0 00648 DO 110 JE = 1, M 00649 IF( IBLOCK( JE ).NE.0 .AND. 00650 $ ( W( JE ).GE.WKILL .OR. IW.EQ.0 ) ) THEN 00651 IW = JE 00652 WKILL = W( JE ) 00653 END IF 00654 110 CONTINUE 00655 IBLOCK( IW ) = 0 00656 120 CONTINUE 00657 END IF 00658 * Now erase all eigenvalues with IBLOCK set to zero 00659 IM = 0 00660 DO 130 JE = 1, M 00661 IF( IBLOCK( JE ).NE.0 ) THEN 00662 IM = IM + 1 00663 W( IM ) = W( JE ) 00664 WERR( IM ) = WERR( JE ) 00665 INDEXW( IM ) = INDEXW( JE ) 00666 IBLOCK( IM ) = IBLOCK( JE ) 00667 END IF 00668 130 CONTINUE 00669 M = IM 00670 END IF 00671 IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN 00672 TOOFEW = .TRUE. 00673 END IF 00674 END IF 00675 * 00676 IF(( IRANGE.EQ.ALLRNG .AND. M.NE.N ).OR. 00677 $ ( IRANGE.EQ.INDRNG .AND. M.NE.IU-IL+1 ) ) THEN 00678 TOOFEW = .TRUE. 00679 END IF 00680 00681 * If ORDER='B', do nothing the eigenvalues are already sorted by 00682 * block. 00683 * If ORDER='E', sort the eigenvalues from smallest to largest 00684 00685 IF( LSAME(ORDER,'E') .AND. NSPLIT.GT.1 ) THEN 00686 DO 150 JE = 1, M - 1 00687 IE = 0 00688 TMP1 = W( JE ) 00689 DO 140 J = JE + 1, M 00690 IF( W( J ).LT.TMP1 ) THEN 00691 IE = J 00692 TMP1 = W( J ) 00693 END IF 00694 140 CONTINUE 00695 IF( IE.NE.0 ) THEN 00696 TMP2 = WERR( IE ) 00697 ITMP1 = IBLOCK( IE ) 00698 ITMP2 = INDEXW( IE ) 00699 W( IE ) = W( JE ) 00700 WERR( IE ) = WERR( JE ) 00701 IBLOCK( IE ) = IBLOCK( JE ) 00702 INDEXW( IE ) = INDEXW( JE ) 00703 W( JE ) = TMP1 00704 WERR( JE ) = TMP2 00705 IBLOCK( JE ) = ITMP1 00706 INDEXW( JE ) = ITMP2 00707 END IF 00708 150 CONTINUE 00709 END IF 00710 * 00711 INFO = 0 00712 IF( NCNVRG ) 00713 $ INFO = INFO + 1 00714 IF( TOOFEW ) 00715 $ INFO = INFO + 2 00716 RETURN 00717 * 00718 * End of DLARRD 00719 * 00720 END
1.7.3 
