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