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