3 /* Table of constant values */
5 static integer c__1 = 1;
7 /* Subroutine */ int dlarrf_(integer *n, doublereal *d__, doublereal *l,
8 doublereal *ld, integer *clstrt, integer *clend, doublereal *w,
9 doublereal *wgap, doublereal *werr, doublereal *spdiam, doublereal *
10 clgapl, doublereal *clgapr, doublereal *pivmin, doublereal *sigma,
11 doublereal *dplus, doublereal *lplus, doublereal *work, integer *info)
13 /* System generated locals */
15 doublereal d__1, d__2, d__3;
17 /* Builtin functions */
18 double sqrt(doublereal);
22 doublereal s, bestshift, smlgrowth, eps, tmp, max1, max2, rrr1, rrr2,
23 znm2, growthbound, fail, fact, oldp;
27 doublereal fail2, avgap, ldmax, rdmax;
29 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
30 doublereal *, integer *);
32 extern doublereal dlamch_(char *);
35 doublereal mingap, lsigma, rdelta;
36 extern logical disnan_(doublereal *);
38 doublereal rsigma, clwdth;
39 logical sawnan1, sawnan2, tryrrr1;
42 /* -- LAPACK auxiliary routine (version 3.1) -- */
43 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
46 /* .. Scalar Arguments .. */
48 /* .. Array Arguments .. */
54 /* Given the initial representation L D L^T and its cluster of close */
55 /* eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ... */
56 /* W( CLEND ), DLARRF finds a new relatively robust representation */
57 /* L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the */
58 /* eigenvalues of L(+) D(+) L(+)^T is relatively isolated. */
63 /* N (input) INTEGER */
64 /* The order of the matrix (subblock, if the matrix splitted). */
66 /* D (input) DOUBLE PRECISION array, dimension (N) */
67 /* The N diagonal elements of the diagonal matrix D. */
69 /* L (input) DOUBLE PRECISION array, dimension (N-1) */
70 /* The (N-1) subdiagonal elements of the unit bidiagonal */
73 /* LD (input) DOUBLE PRECISION array, dimension (N-1) */
74 /* The (N-1) elements L(i)*D(i). */
76 /* CLSTRT (input) INTEGER */
77 /* The index of the first eigenvalue in the cluster. */
79 /* CLEND (input) INTEGER */
80 /* The index of the last eigenvalue in the cluster. */
82 /* W (input) DOUBLE PRECISION array, dimension >= (CLEND-CLSTRT+1) */
83 /* The eigenvalue APPROXIMATIONS of L D L^T in ascending order. */
84 /* W( CLSTRT ) through W( CLEND ) form the cluster of relatively */
85 /* close eigenalues. */
87 /* WGAP (input/output) DOUBLE PRECISION array, dimension >= (CLEND-CLSTRT+1) */
88 /* The separation from the right neighbor eigenvalue in W. */
90 /* WERR (input) DOUBLE PRECISION array, dimension >= (CLEND-CLSTRT+1) */
91 /* WERR contain the semiwidth of the uncertainty */
92 /* interval of the corresponding eigenvalue APPROXIMATION in W */
94 /* SPDIAM (input) estimate of the spectral diameter obtained from the */
95 /* Gerschgorin intervals */
97 /* CLGAPL, CLGAPR (input) absolute gap on each end of the cluster. */
98 /* Set by the calling routine to protect against shifts too close */
99 /* to eigenvalues outside the cluster. */
101 /* PIVMIN (input) DOUBLE PRECISION */
102 /* The minimum pivot allowed in the Sturm sequence. */
104 /* SIGMA (output) DOUBLE PRECISION */
105 /* The shift used to form L(+) D(+) L(+)^T. */
107 /* DPLUS (output) DOUBLE PRECISION array, dimension (N) */
108 /* The N diagonal elements of the diagonal matrix D(+). */
110 /* LPLUS (output) DOUBLE PRECISION array, dimension (N-1) */
111 /* The first (N-1) elements of LPLUS contain the subdiagonal */
112 /* elements of the unit bidiagonal matrix L(+). */
114 /* WORK (workspace) DOUBLE PRECISION array, dimension (2*N) */
117 /* Further Details */
118 /* =============== */
120 /* Based on contributions by */
121 /* Beresford Parlett, University of California, Berkeley, USA */
122 /* Jim Demmel, University of California, Berkeley, USA */
123 /* Inderjit Dhillon, University of Texas, Austin, USA */
124 /* Osni Marques, LBNL/NERSC, USA */
125 /* Christof Voemel, University of California, Berkeley, USA */
127 /* ===================================================================== */
129 /* .. Parameters .. */
131 /* .. Local Scalars .. */
133 /* .. External Functions .. */
135 /* .. External Subroutines .. */
137 /* .. Intrinsic Functions .. */
139 /* .. Executable Statements .. */
141 /* Parameter adjustments */
155 eps = dlamch_("Precision");
158 /* Note that we cannot guarantee that for any of the shifts tried, */
159 /* the factorization has a small or even moderate element growth. */
160 /* There could be Ritz values at both ends of the cluster and despite */
161 /* backing off, there are examples where all factorizations tried */
162 /* (in IEEE mode, allowing zero pivots & infinities) have INFINITE */
163 /* element growth. */
164 /* For this reason, we should use PIVMIN in this subroutine so that at */
165 /* least the L D L^T factorization exists. It can be checked afterwards */
166 /* whether the element growth caused bad residuals/orthogonality. */
167 /* Decide whether the code should accept the best among all */
168 /* representations despite large element growth or signal INFO=1 */
171 /* Compute the average gap length of the cluster */
172 clwdth = (d__1 = w[*clend] - w[*clstrt], abs(d__1)) + werr[*clend] + werr[
174 avgap = clwdth / (doublereal) (*clend - *clstrt);
175 mingap = min(*clgapl,*clgapr);
176 /* Initial values for shifts to both ends of cluster */
178 d__1 = w[*clstrt], d__2 = w[*clend];
179 lsigma = min(d__1,d__2) - werr[*clstrt];
181 d__1 = w[*clstrt], d__2 = w[*clend];
182 rsigma = max(d__1,d__2) + werr[*clend];
183 /* Use a small fudge to make sure that we really shift to the outside */
184 lsigma -= abs(lsigma) * 4. * eps;
185 rsigma += abs(rsigma) * 4. * eps;
186 /* Compute upper bounds for how much to back off the initial shifts */
187 ldmax = mingap * .25 + *pivmin * 2.;
188 rdmax = mingap * .25 + *pivmin * 2.;
190 d__1 = avgap, d__2 = wgap[*clstrt];
191 ldelta = max(d__1,d__2) / fact;
193 d__1 = avgap, d__2 = wgap[*clend - 1];
194 rdelta = max(d__1,d__2) / fact;
196 /* Initialize the record of the best representation found */
200 fail = (doublereal) (*n - 1) * mingap / (*spdiam * eps);
201 fail2 = (doublereal) (*n - 1) * mingap / (*spdiam * sqrt(eps));
204 /* while (KTRY <= KTRYMAX) */
206 growthbound = *spdiam * 8.;
210 /* Ensure that we do not back off too much of the initial shifts */
211 ldelta = min(ldmax,ldelta);
212 rdelta = min(rdmax,rdelta);
213 /* Compute the element growth when shifting to both ends of the cluster */
214 /* accept the shift if there is no element growth at one of the two ends */
217 dplus[1] = d__[1] + s;
218 if (abs(dplus[1]) < *pivmin) {
219 dplus[1] = -(*pivmin);
220 /* Need to set SAWNAN1 because refined RRR test should not be used */
224 max1 = abs(dplus[1]);
226 for (i__ = 1; i__ <= i__1; ++i__) {
227 lplus[i__] = ld[i__] / dplus[i__];
228 s = s * lplus[i__] * l[i__] - lsigma;
229 dplus[i__ + 1] = d__[i__ + 1] + s;
230 if ((d__1 = dplus[i__ + 1], abs(d__1)) < *pivmin) {
231 dplus[i__ + 1] = -(*pivmin);
232 /* Need to set SAWNAN1 because refined RRR test should not be used */
237 d__2 = max1, d__3 = (d__1 = dplus[i__ + 1], abs(d__1));
238 max1 = max(d__2,d__3);
241 sawnan1 = sawnan1 || disnan_(&max1);
242 if (forcer || max1 <= growthbound && ! sawnan1) {
249 work[1] = d__[1] + s;
250 if (abs(work[1]) < *pivmin) {
251 work[1] = -(*pivmin);
252 /* Need to set SAWNAN2 because refined RRR test should not be used */
258 for (i__ = 1; i__ <= i__1; ++i__) {
259 work[*n + i__] = ld[i__] / work[i__];
260 s = s * work[*n + i__] * l[i__] - rsigma;
261 work[i__ + 1] = d__[i__ + 1] + s;
262 if ((d__1 = work[i__ + 1], abs(d__1)) < *pivmin) {
263 work[i__ + 1] = -(*pivmin);
264 /* Need to set SAWNAN2 because refined RRR test should not be used */
269 d__2 = max2, d__3 = (d__1 = work[i__ + 1], abs(d__1));
270 max2 = max(d__2,d__3);
273 sawnan2 = sawnan2 || disnan_(&max2);
274 if (forcer || max2 <= growthbound && ! sawnan2) {
279 /* If we are at this point, both shifts led to too much element growth */
280 /* Record the better of the two shifts (provided it didn't lead to NaN) */
281 if (sawnan1 && sawnan2) {
282 /* both MAX1 and MAX2 are NaN */
287 if (max1 <= smlgrowth) {
293 if (sawnan1 || max2 <= max1) {
296 if (max2 <= smlgrowth) {
302 /* If we are here, both the left and the right shift led to */
303 /* element growth. If the element growth is moderate, then */
304 /* we may still accept the representation, if it passes a */
305 /* refined test for RRR. This test supposes that no NaN occurred. */
306 /* Moreover, we use the refined RRR test only for isolated clusters. */
307 if (clwdth < mingap / 128. && min(max1,max2) < fail2 && ! sawnan1 && !
314 if (tryrrr1 && dorrr1) {
316 tmp = (d__1 = dplus[*n], abs(d__1));
320 for (i__ = *n - 1; i__ >= 1; --i__) {
322 prod = dplus[i__ + 1] * work[*n + i__ + 1] / (dplus[i__] *
323 work[*n + i__]) * oldp;
325 prod *= (d__1 = work[*n + i__], abs(d__1));
328 /* Computing 2nd power */
332 d__2 = tmp, d__3 = (d__1 = dplus[i__] * prod, abs(d__1));
333 tmp = max(d__2,d__3);
336 rrr1 = tmp / (*spdiam * sqrt(znm2));
342 } else if (indx == 2) {
343 tmp = (d__1 = work[*n], abs(d__1));
347 for (i__ = *n - 1; i__ >= 1; --i__) {
349 prod = work[i__ + 1] * lplus[i__ + 1] / (work[i__] *
352 prod *= (d__1 = lplus[i__], abs(d__1));
355 /* Computing 2nd power */
359 d__2 = tmp, d__3 = (d__1 = work[i__] * prod, abs(d__1));
360 tmp = max(d__2,d__3);
363 rrr2 = tmp / (*spdiam * sqrt(znm2));
373 /* If we are here, both shifts failed also the RRR test. */
374 /* Back off to the outside */
376 d__1 = lsigma - ldelta, d__2 = lsigma - ldmax;
377 lsigma = max(d__1,d__2);
379 d__1 = rsigma + rdelta, d__2 = rsigma + rdmax;
380 rsigma = min(d__1,d__2);
386 /* None of the representations investigated satisfied our */
387 /* criteria. Take the best one we found. */
388 if (smlgrowth < fail || nofail) {
400 } else if (shift == 2) {
401 /* store new L and D back into DPLUS, LPLUS */
402 dcopy_(n, &work[1], &c__1, &dplus[1], &c__1);
404 dcopy_(&i__1, &work[*n + 1], &c__1, &lplus[1], &c__1);