|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST, 00002 $ RTOL, OFFSET, W, WERR, WORK, IWORK, 00003 $ PIVMIN, SPDIAM, INFO ) 00004 * 00005 * -- LAPACK auxiliary routine (version 3.2.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * June 2010 00009 * 00010 * .. Scalar Arguments .. 00011 INTEGER IFIRST, ILAST, INFO, N, OFFSET 00012 DOUBLE PRECISION PIVMIN, RTOL, SPDIAM 00013 * .. 00014 * .. Array Arguments .. 00015 INTEGER IWORK( * ) 00016 DOUBLE PRECISION D( * ), E2( * ), W( * ), 00017 $ WERR( * ), WORK( * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * Given the initial eigenvalue approximations of T, DLARRJ 00024 * does bisection to refine the eigenvalues of T, 00025 * W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial 00026 * guesses for these eigenvalues are input in W, the corresponding estimate 00027 * of the error in these guesses in WERR. During bisection, intervals 00028 * [left, right] are maintained by storing their mid-points and 00029 * semi-widths in the arrays W and WERR respectively. 00030 * 00031 * Arguments 00032 * ========= 00033 * 00034 * N (input) INTEGER 00035 * The order of the matrix. 00036 * 00037 * D (input) DOUBLE PRECISION array, dimension (N) 00038 * The N diagonal elements of T. 00039 * 00040 * E2 (input) DOUBLE PRECISION array, dimension (N-1) 00041 * The Squares of the (N-1) subdiagonal elements of T. 00042 * 00043 * IFIRST (input) INTEGER 00044 * The index of the first eigenvalue to be computed. 00045 * 00046 * ILAST (input) INTEGER 00047 * The index of the last eigenvalue to be computed. 00048 * 00049 * RTOL (input) DOUBLE PRECISION 00050 * Tolerance for the convergence of the bisection intervals. 00051 * An interval [LEFT,RIGHT] has converged if 00052 * RIGHT-LEFT.LT.RTOL*MAX(|LEFT|,|RIGHT|). 00053 * 00054 * OFFSET (input) INTEGER 00055 * Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET 00056 * through ILAST-OFFSET elements of these arrays are to be used. 00057 * 00058 * W (input/output) DOUBLE PRECISION array, dimension (N) 00059 * On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are 00060 * estimates of the eigenvalues of L D L^T indexed IFIRST through 00061 * ILAST. 00062 * On output, these estimates are refined. 00063 * 00064 * WERR (input/output) DOUBLE PRECISION array, dimension (N) 00065 * On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are 00066 * the errors in the estimates of the corresponding elements in W. 00067 * On output, these errors are refined. 00068 * 00069 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N) 00070 * Workspace. 00071 * 00072 * IWORK (workspace) INTEGER array, dimension (2*N) 00073 * Workspace. 00074 * 00075 * PIVMIN (input) DOUBLE PRECISION 00076 * The minimum pivot in the Sturm sequence for T. 00077 * 00078 * SPDIAM (input) DOUBLE PRECISION 00079 * The spectral diameter of T. 00080 * 00081 * INFO (output) INTEGER 00082 * Error flag. 00083 * 00084 * Further Details 00085 * =============== 00086 * 00087 * Based on contributions by 00088 * Beresford Parlett, University of California, Berkeley, USA 00089 * Jim Demmel, University of California, Berkeley, USA 00090 * Inderjit Dhillon, University of Texas, Austin, USA 00091 * Osni Marques, LBNL/NERSC, USA 00092 * Christof Voemel, University of California, Berkeley, USA 00093 * 00094 * ===================================================================== 00095 * 00096 * .. Parameters .. 00097 DOUBLE PRECISION ZERO, ONE, TWO, HALF 00098 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 00099 $ HALF = 0.5D0 ) 00100 INTEGER MAXITR 00101 * .. 00102 * .. Local Scalars .. 00103 INTEGER CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT, 00104 $ OLNINT, P, PREV, SAVI1 00105 DOUBLE PRECISION DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH 00106 * 00107 * .. 00108 * .. Intrinsic Functions .. 00109 INTRINSIC ABS, MAX 00110 * .. 00111 * .. Executable Statements .. 00112 * 00113 INFO = 0 00114 * 00115 MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) / 00116 $ LOG( TWO ) ) + 2 00117 * 00118 * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ]. 00119 * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while 00120 * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 ) 00121 * for an unconverged interval is set to the index of the next unconverged 00122 * interval, and is -1 or 0 for a converged interval. Thus a linked 00123 * list of unconverged intervals is set up. 00124 * 00125 00126 I1 = IFIRST 00127 I2 = ILAST 00128 * The number of unconverged intervals 00129 NINT = 0 00130 * The last unconverged interval found 00131 PREV = 0 00132 DO 75 I = I1, I2 00133 K = 2*I 00134 II = I - OFFSET 00135 LEFT = W( II ) - WERR( II ) 00136 MID = W(II) 00137 RIGHT = W( II ) + WERR( II ) 00138 WIDTH = RIGHT - MID 00139 TMP = MAX( ABS( LEFT ), ABS( RIGHT ) ) 00140 00141 * The following test prevents the test of converged intervals 00142 IF( WIDTH.LT.RTOL*TMP ) THEN 00143 * This interval has already converged and does not need refinement. 00144 * (Note that the gaps might change through refining the 00145 * eigenvalues, however, they can only get bigger.) 00146 * Remove it from the list. 00147 IWORK( K-1 ) = -1 00148 * Make sure that I1 always points to the first unconverged interval 00149 IF((I.EQ.I1).AND.(I.LT.I2)) I1 = I + 1 00150 IF((PREV.GE.I1).AND.(I.LE.I2)) IWORK( 2*PREV-1 ) = I + 1 00151 ELSE 00152 * unconverged interval found 00153 PREV = I 00154 * Make sure that [LEFT,RIGHT] contains the desired eigenvalue 00155 * 00156 * Do while( CNT(LEFT).GT.I-1 ) 00157 * 00158 FAC = ONE 00159 20 CONTINUE 00160 CNT = 0 00161 S = LEFT 00162 DPLUS = D( 1 ) - S 00163 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 00164 DO 30 J = 2, N 00165 DPLUS = D( J ) - S - E2( J-1 )/DPLUS 00166 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 00167 30 CONTINUE 00168 IF( CNT.GT.I-1 ) THEN 00169 LEFT = LEFT - WERR( II )*FAC 00170 FAC = TWO*FAC 00171 GO TO 20 00172 END IF 00173 * 00174 * Do while( CNT(RIGHT).LT.I ) 00175 * 00176 FAC = ONE 00177 50 CONTINUE 00178 CNT = 0 00179 S = RIGHT 00180 DPLUS = D( 1 ) - S 00181 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 00182 DO 60 J = 2, N 00183 DPLUS = D( J ) - S - E2( J-1 )/DPLUS 00184 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 00185 60 CONTINUE 00186 IF( CNT.LT.I ) THEN 00187 RIGHT = RIGHT + WERR( II )*FAC 00188 FAC = TWO*FAC 00189 GO TO 50 00190 END IF 00191 NINT = NINT + 1 00192 IWORK( K-1 ) = I + 1 00193 IWORK( K ) = CNT 00194 END IF 00195 WORK( K-1 ) = LEFT 00196 WORK( K ) = RIGHT 00197 75 CONTINUE 00198 00199 00200 SAVI1 = I1 00201 * 00202 * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals 00203 * and while (ITER.LT.MAXITR) 00204 * 00205 ITER = 0 00206 80 CONTINUE 00207 PREV = I1 - 1 00208 I = I1 00209 OLNINT = NINT 00210 00211 DO 100 P = 1, OLNINT 00212 K = 2*I 00213 II = I - OFFSET 00214 NEXT = IWORK( K-1 ) 00215 LEFT = WORK( K-1 ) 00216 RIGHT = WORK( K ) 00217 MID = HALF*( LEFT + RIGHT ) 00218 00219 * semiwidth of interval 00220 WIDTH = RIGHT - MID 00221 TMP = MAX( ABS( LEFT ), ABS( RIGHT ) ) 00222 00223 IF( ( WIDTH.LT.RTOL*TMP ) .OR. 00224 $ (ITER.EQ.MAXITR) )THEN 00225 * reduce number of unconverged intervals 00226 NINT = NINT - 1 00227 * Mark interval as converged. 00228 IWORK( K-1 ) = 0 00229 IF( I1.EQ.I ) THEN 00230 I1 = NEXT 00231 ELSE 00232 * Prev holds the last unconverged interval previously examined 00233 IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT 00234 END IF 00235 I = NEXT 00236 GO TO 100 00237 END IF 00238 PREV = I 00239 * 00240 * Perform one bisection step 00241 * 00242 CNT = 0 00243 S = MID 00244 DPLUS = D( 1 ) - S 00245 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 00246 DO 90 J = 2, N 00247 DPLUS = D( J ) - S - E2( J-1 )/DPLUS 00248 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 00249 90 CONTINUE 00250 IF( CNT.LE.I-1 ) THEN 00251 WORK( K-1 ) = MID 00252 ELSE 00253 WORK( K ) = MID 00254 END IF 00255 I = NEXT 00256 00257 100 CONTINUE 00258 ITER = ITER + 1 00259 * do another loop if there are still unconverged intervals 00260 * However, in the last iteration, all intervals are accepted 00261 * since this is the best we can do. 00262 IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80 00263 * 00264 * 00265 * At this point, all the intervals have converged 00266 DO 110 I = SAVI1, ILAST 00267 K = 2*I 00268 II = I - OFFSET 00269 * All intervals marked by '0' have been refined. 00270 IF( IWORK( K-1 ).EQ.0 ) THEN 00271 W( II ) = HALF*( WORK( K-1 )+WORK( K ) ) 00272 WERR( II ) = WORK( K ) - W( II ) 00273 END IF 00274 110 CONTINUE 00275 * 00276 00277 RETURN 00278 * 00279 * End of DLARRJ 00280 * 00281 END
1.7.3 
