00001 /* dlarrb.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Subroutine */ int dlarrb_(integer *n, doublereal *d__, doublereal *lld, 00017 integer *ifirst, integer *ilast, doublereal *rtol1, doublereal *rtol2, 00018 integer *offset, doublereal *w, doublereal *wgap, doublereal *werr, 00019 doublereal *work, integer *iwork, doublereal *pivmin, doublereal * 00020 spdiam, integer *twist, integer *info) 00021 { 00022 /* System generated locals */ 00023 integer i__1; 00024 doublereal d__1, d__2; 00025 00026 /* Builtin functions */ 00027 double log(doublereal); 00028 00029 /* Local variables */ 00030 integer i__, k, r__, i1, ii, ip; 00031 doublereal gap, mid, tmp, back, lgap, rgap, left; 00032 integer iter, nint, prev, next; 00033 doublereal cvrgd, right, width; 00034 extern integer dlaneg_(integer *, doublereal *, doublereal *, doublereal * 00035 , doublereal *, integer *); 00036 integer negcnt; 00037 doublereal mnwdth; 00038 integer olnint, maxitr; 00039 00040 00041 /* -- LAPACK auxiliary routine (version 3.2) -- */ 00042 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00043 /* November 2006 */ 00044 00045 /* .. Scalar Arguments .. */ 00046 /* .. */ 00047 /* .. Array Arguments .. */ 00048 /* .. */ 00049 00050 /* Purpose */ 00051 /* ======= */ 00052 00053 /* Given the relatively robust representation(RRR) L D L^T, DLARRB */ 00054 /* does "limited" bisection to refine the eigenvalues of L D L^T, */ 00055 /* W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial */ 00056 /* guesses for these eigenvalues are input in W, the corresponding estimate */ 00057 /* of the error in these guesses and their gaps are input in WERR */ 00058 /* and WGAP, respectively. During bisection, intervals */ 00059 /* [left, right] are maintained by storing their mid-points and */ 00060 /* semi-widths in the arrays W and WERR respectively. */ 00061 00062 /* Arguments */ 00063 /* ========= */ 00064 00065 /* N (input) INTEGER */ 00066 /* The order of the matrix. */ 00067 00068 /* D (input) DOUBLE PRECISION array, dimension (N) */ 00069 /* The N diagonal elements of the diagonal matrix D. */ 00070 00071 /* LLD (input) DOUBLE PRECISION array, dimension (N-1) */ 00072 /* The (N-1) elements L(i)*L(i)*D(i). */ 00073 00074 /* IFIRST (input) INTEGER */ 00075 /* The index of the first eigenvalue to be computed. */ 00076 00077 /* ILAST (input) INTEGER */ 00078 /* The index of the last eigenvalue to be computed. */ 00079 00080 /* RTOL1 (input) DOUBLE PRECISION */ 00081 /* RTOL2 (input) DOUBLE PRECISION */ 00082 /* Tolerance for the convergence of the bisection intervals. */ 00083 /* An interval [LEFT,RIGHT] has converged if */ 00084 /* RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) */ 00085 /* where GAP is the (estimated) distance to the nearest */ 00086 /* eigenvalue. */ 00087 00088 /* OFFSET (input) INTEGER */ 00089 /* Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET */ 00090 /* through ILAST-OFFSET elements of these arrays are to be used. */ 00091 00092 /* W (input/output) DOUBLE PRECISION array, dimension (N) */ 00093 /* On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are */ 00094 /* estimates of the eigenvalues of L D L^T indexed IFIRST throug */ 00095 /* ILAST. */ 00096 /* On output, these estimates are refined. */ 00097 00098 /* WGAP (input/output) DOUBLE PRECISION array, dimension (N-1) */ 00099 /* On input, the (estimated) gaps between consecutive */ 00100 /* eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between */ 00101 /* eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST */ 00102 /* then WGAP(IFIRST-OFFSET) must be set to ZERO. */ 00103 /* On output, these gaps are refined. */ 00104 00105 /* WERR (input/output) DOUBLE PRECISION array, dimension (N) */ 00106 /* On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are */ 00107 /* the errors in the estimates of the corresponding elements in W. */ 00108 /* On output, these errors are refined. */ 00109 00110 /* WORK (workspace) DOUBLE PRECISION array, dimension (2*N) */ 00111 /* Workspace. */ 00112 00113 /* IWORK (workspace) INTEGER array, dimension (2*N) */ 00114 /* Workspace. */ 00115 00116 /* PIVMIN (input) DOUBLE PRECISION */ 00117 /* The minimum pivot in the Sturm sequence. */ 00118 00119 /* SPDIAM (input) DOUBLE PRECISION */ 00120 /* The spectral diameter of the matrix. */ 00121 00122 /* TWIST (input) INTEGER */ 00123 /* The twist index for the twisted factorization that is used */ 00124 /* for the negcount. */ 00125 /* TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T */ 00126 /* TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T */ 00127 /* TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r) */ 00128 00129 /* INFO (output) INTEGER */ 00130 /* Error flag. */ 00131 00132 /* Further Details */ 00133 /* =============== */ 00134 00135 /* Based on contributions by */ 00136 /* Beresford Parlett, University of California, Berkeley, USA */ 00137 /* Jim Demmel, University of California, Berkeley, USA */ 00138 /* Inderjit Dhillon, University of Texas, Austin, USA */ 00139 /* Osni Marques, LBNL/NERSC, USA */ 00140 /* Christof Voemel, University of California, Berkeley, USA */ 00141 00142 /* ===================================================================== */ 00143 00144 /* .. Parameters .. */ 00145 /* .. */ 00146 /* .. Local Scalars .. */ 00147 /* .. */ 00148 /* .. External Functions .. */ 00149 00150 /* .. */ 00151 /* .. Intrinsic Functions .. */ 00152 /* .. */ 00153 /* .. Executable Statements .. */ 00154 00155 /* Parameter adjustments */ 00156 --iwork; 00157 --work; 00158 --werr; 00159 --wgap; 00160 --w; 00161 --lld; 00162 --d__; 00163 00164 /* Function Body */ 00165 *info = 0; 00166 00167 maxitr = (integer) ((log(*spdiam + *pivmin) - log(*pivmin)) / log(2.)) + 00168 2; 00169 mnwdth = *pivmin * 2.; 00170 00171 r__ = *twist; 00172 if (r__ < 1 || r__ > *n) { 00173 r__ = *n; 00174 } 00175 00176 /* Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ]. */ 00177 /* The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while */ 00178 /* Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 ) */ 00179 /* for an unconverged interval is set to the index of the next unconverged */ 00180 /* interval, and is -1 or 0 for a converged interval. Thus a linked */ 00181 /* list of unconverged intervals is set up. */ 00182 00183 i1 = *ifirst; 00184 /* The number of unconverged intervals */ 00185 nint = 0; 00186 /* The last unconverged interval found */ 00187 prev = 0; 00188 rgap = wgap[i1 - *offset]; 00189 i__1 = *ilast; 00190 for (i__ = i1; i__ <= i__1; ++i__) { 00191 k = i__ << 1; 00192 ii = i__ - *offset; 00193 left = w[ii] - werr[ii]; 00194 right = w[ii] + werr[ii]; 00195 lgap = rgap; 00196 rgap = wgap[ii]; 00197 gap = min(lgap,rgap); 00198 /* Make sure that [LEFT,RIGHT] contains the desired eigenvalue */ 00199 /* Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT */ 00200 00201 /* Do while( NEGCNT(LEFT).GT.I-1 ) */ 00202 00203 back = werr[ii]; 00204 L20: 00205 negcnt = dlaneg_(n, &d__[1], &lld[1], &left, pivmin, &r__); 00206 if (negcnt > i__ - 1) { 00207 left -= back; 00208 back *= 2.; 00209 goto L20; 00210 } 00211 00212 /* Do while( NEGCNT(RIGHT).LT.I ) */ 00213 /* Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT */ 00214 00215 back = werr[ii]; 00216 L50: 00217 negcnt = dlaneg_(n, &d__[1], &lld[1], &right, pivmin, &r__); 00218 if (negcnt < i__) { 00219 right += back; 00220 back *= 2.; 00221 goto L50; 00222 } 00223 width = (d__1 = left - right, abs(d__1)) * .5; 00224 /* Computing MAX */ 00225 d__1 = abs(left), d__2 = abs(right); 00226 tmp = max(d__1,d__2); 00227 /* Computing MAX */ 00228 d__1 = *rtol1 * gap, d__2 = *rtol2 * tmp; 00229 cvrgd = max(d__1,d__2); 00230 if (width <= cvrgd || width <= mnwdth) { 00231 /* This interval has already converged and does not need refinement. */ 00232 /* (Note that the gaps might change through refining the */ 00233 /* eigenvalues, however, they can only get bigger.) */ 00234 /* Remove it from the list. */ 00235 iwork[k - 1] = -1; 00236 /* Make sure that I1 always points to the first unconverged interval */ 00237 if (i__ == i1 && i__ < *ilast) { 00238 i1 = i__ + 1; 00239 } 00240 if (prev >= i1 && i__ <= *ilast) { 00241 iwork[(prev << 1) - 1] = i__ + 1; 00242 } 00243 } else { 00244 /* unconverged interval found */ 00245 prev = i__; 00246 ++nint; 00247 iwork[k - 1] = i__ + 1; 00248 iwork[k] = negcnt; 00249 } 00250 work[k - 1] = left; 00251 work[k] = right; 00252 /* L75: */ 00253 } 00254 00255 /* Do while( NINT.GT.0 ), i.e. there are still unconverged intervals */ 00256 /* and while (ITER.LT.MAXITR) */ 00257 00258 iter = 0; 00259 L80: 00260 prev = i1 - 1; 00261 i__ = i1; 00262 olnint = nint; 00263 i__1 = olnint; 00264 for (ip = 1; ip <= i__1; ++ip) { 00265 k = i__ << 1; 00266 ii = i__ - *offset; 00267 rgap = wgap[ii]; 00268 lgap = rgap; 00269 if (ii > 1) { 00270 lgap = wgap[ii - 1]; 00271 } 00272 gap = min(lgap,rgap); 00273 next = iwork[k - 1]; 00274 left = work[k - 1]; 00275 right = work[k]; 00276 mid = (left + right) * .5; 00277 /* semiwidth of interval */ 00278 width = right - mid; 00279 /* Computing MAX */ 00280 d__1 = abs(left), d__2 = abs(right); 00281 tmp = max(d__1,d__2); 00282 /* Computing MAX */ 00283 d__1 = *rtol1 * gap, d__2 = *rtol2 * tmp; 00284 cvrgd = max(d__1,d__2); 00285 if (width <= cvrgd || width <= mnwdth || iter == maxitr) { 00286 /* reduce number of unconverged intervals */ 00287 --nint; 00288 /* Mark interval as converged. */ 00289 iwork[k - 1] = 0; 00290 if (i1 == i__) { 00291 i1 = next; 00292 } else { 00293 /* Prev holds the last unconverged interval previously examined */ 00294 if (prev >= i1) { 00295 iwork[(prev << 1) - 1] = next; 00296 } 00297 } 00298 i__ = next; 00299 goto L100; 00300 } 00301 prev = i__; 00302 00303 /* Perform one bisection step */ 00304 00305 negcnt = dlaneg_(n, &d__[1], &lld[1], &mid, pivmin, &r__); 00306 if (negcnt <= i__ - 1) { 00307 work[k - 1] = mid; 00308 } else { 00309 work[k] = mid; 00310 } 00311 i__ = next; 00312 L100: 00313 ; 00314 } 00315 ++iter; 00316 /* do another loop if there are still unconverged intervals */ 00317 /* However, in the last iteration, all intervals are accepted */ 00318 /* since this is the best we can do. */ 00319 if (nint > 0 && iter <= maxitr) { 00320 goto L80; 00321 } 00322 00323 00324 /* At this point, all the intervals have converged */ 00325 i__1 = *ilast; 00326 for (i__ = *ifirst; i__ <= i__1; ++i__) { 00327 k = i__ << 1; 00328 ii = i__ - *offset; 00329 /* All intervals marked by '0' have been refined. */ 00330 if (iwork[k - 1] == 0) { 00331 w[ii] = (work[k - 1] + work[k]) * .5; 00332 werr[ii] = work[k] - w[ii]; 00333 } 00334 /* L110: */ 00335 } 00336 00337 i__1 = *ilast; 00338 for (i__ = *ifirst + 1; i__ <= i__1; ++i__) { 00339 k = i__ << 1; 00340 ii = i__ - *offset; 00341 /* Computing MAX */ 00342 d__1 = 0., d__2 = w[ii] - werr[ii] - w[ii - 1] - werr[ii - 1]; 00343 wgap[ii - 1] = max(d__1,d__2); 00344 /* L111: */ 00345 } 00346 return 0; 00347 00348 /* End of DLARRB */ 00349 00350 } /* dlarrb_ */