3 /* Table of constant values */
5 static integer c__1 = 1;
7 /* Subroutine */ int slarrf_(integer *n, real *d__, real *l, real *ld,
8 integer *clstrt, integer *clend, real *w, real *wgap, real *werr,
9 real *spdiam, real *clgapl, real *clgapr, real *pivmin, real *sigma,
10 real *dplus, real *lplus, real *work, integer *info)
12 /* System generated locals */
14 real r__1, r__2, r__3;
16 /* Builtin functions */
17 double sqrt(doublereal);
21 real s, bestshift, smlgrowth, eps, tmp, max1, max2, rrr1, rrr2, znm2,
22 growthbound, fail, fact, oldp;
26 real fail2, avgap, ldmax, rdmax;
28 extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
32 extern doublereal slamch_(char *);
34 real mingap, lsigma, rdelta;
37 extern logical sisnan_(real *);
38 logical sawnan1, sawnan2, tryrrr1;
41 /* -- LAPACK auxiliary routine (version 3.1) -- */
42 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
45 /* .. Scalar Arguments .. */
47 /* .. Array Arguments .. */
53 /* Given the initial representation L D L^T and its cluster of close */
54 /* eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ... */
55 /* W( CLEND ), SLARRF finds a new relatively robust representation */
56 /* L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the */
57 /* eigenvalues of L(+) D(+) L(+)^T is relatively isolated. */
62 /* N (input) INTEGER */
63 /* The order of the matrix (subblock, if the matrix splitted). */
65 /* D (input) REAL array, dimension (N) */
66 /* The N diagonal elements of the diagonal matrix D. */
68 /* L (input) REAL array, dimension (N-1) */
69 /* The (N-1) subdiagonal elements of the unit bidiagonal */
72 /* LD (input) REAL array, dimension (N-1) */
73 /* The (N-1) elements L(i)*D(i). */
75 /* CLSTRT (input) INTEGER */
76 /* The index of the first eigenvalue in the cluster. */
78 /* CLEND (input) INTEGER */
79 /* The index of the last eigenvalue in the cluster. */
81 /* W (input) REAL array, dimension >= (CLEND-CLSTRT+1) */
82 /* The eigenvalue APPROXIMATIONS of L D L^T in ascending order. */
83 /* W( CLSTRT ) through W( CLEND ) form the cluster of relatively */
84 /* close eigenalues. */
86 /* WGAP (input/output) REAL array, dimension >= (CLEND-CLSTRT+1) */
87 /* The separation from the right neighbor eigenvalue in W. */
89 /* WERR (input) REAL array, dimension >= (CLEND-CLSTRT+1) */
90 /* WERR contain the semiwidth of the uncertainty */
91 /* interval of the corresponding eigenvalue APPROXIMATION in W */
93 /* SPDIAM (input) estimate of the spectral diameter obtained from the */
94 /* Gerschgorin intervals */
96 /* CLGAPL, CLGAPR (input) absolute gap on each end of the cluster. */
97 /* Set by the calling routine to protect against shifts too close */
98 /* to eigenvalues outside the cluster. */
100 /* PIVMIN (input) DOUBLE PRECISION */
101 /* The minimum pivot allowed in the Sturm sequence. */
103 /* SIGMA (output) REAL */
104 /* The shift used to form L(+) D(+) L(+)^T. */
106 /* DPLUS (output) REAL array, dimension (N) */
107 /* The N diagonal elements of the diagonal matrix D(+). */
109 /* LPLUS (output) REAL array, dimension (N-1) */
110 /* The first (N-1) elements of LPLUS contain the subdiagonal */
111 /* elements of the unit bidiagonal matrix L(+). */
113 /* WORK (workspace) REAL array, dimension (2*N) */
116 /* Further Details */
117 /* =============== */
119 /* Based on contributions by */
120 /* Beresford Parlett, University of California, Berkeley, USA */
121 /* Jim Demmel, University of California, Berkeley, USA */
122 /* Inderjit Dhillon, University of Texas, Austin, USA */
123 /* Osni Marques, LBNL/NERSC, USA */
124 /* Christof Voemel, University of California, Berkeley, USA */
126 /* ===================================================================== */
128 /* .. Parameters .. */
130 /* .. Local Scalars .. */
132 /* .. External Functions .. */
134 /* .. External Subroutines .. */
136 /* .. Intrinsic Functions .. */
138 /* .. Executable Statements .. */
140 /* Parameter adjustments */
154 eps = slamch_("Precision");
157 /* Note that we cannot guarantee that for any of the shifts tried, */
158 /* the factorization has a small or even moderate element growth. */
159 /* There could be Ritz values at both ends of the cluster and despite */
160 /* backing off, there are examples where all factorizations tried */
161 /* (in IEEE mode, allowing zero pivots & infinities) have INFINITE */
162 /* element growth. */
163 /* For this reason, we should use PIVMIN in this subroutine so that at */
164 /* least the L D L^T factorization exists. It can be checked afterwards */
165 /* whether the element growth caused bad residuals/orthogonality. */
166 /* Decide whether the code should accept the best among all */
167 /* representations despite large element growth or signal INFO=1 */
170 /* Compute the average gap length of the cluster */
171 clwdth = (r__1 = w[*clend] - w[*clstrt], dabs(r__1)) + werr[*clend] +
173 avgap = clwdth / (real) (*clend - *clstrt);
174 mingap = dmin(*clgapl,*clgapr);
175 /* Initial values for shifts to both ends of cluster */
177 r__1 = w[*clstrt], r__2 = w[*clend];
178 lsigma = dmin(r__1,r__2) - werr[*clstrt];
180 r__1 = w[*clstrt], r__2 = w[*clend];
181 rsigma = dmax(r__1,r__2) + werr[*clend];
182 /* Use a small fudge to make sure that we really shift to the outside */
183 lsigma -= dabs(lsigma) * 2.f * eps;
184 rsigma += dabs(rsigma) * 2.f * eps;
185 /* Compute upper bounds for how much to back off the initial shifts */
186 ldmax = mingap * .25f + *pivmin * 2.f;
187 rdmax = mingap * .25f + *pivmin * 2.f;
189 r__1 = avgap, r__2 = wgap[*clstrt];
190 ldelta = dmax(r__1,r__2) / fact;
192 r__1 = avgap, r__2 = wgap[*clend - 1];
193 rdelta = dmax(r__1,r__2) / fact;
195 /* Initialize the record of the best representation found */
199 fail = (real) (*n - 1) * mingap / (*spdiam * eps);
200 fail2 = (real) (*n - 1) * mingap / (*spdiam * sqrt(eps));
203 /* while (KTRY <= KTRYMAX) */
205 growthbound = *spdiam * 8.f;
209 /* Ensure that we do not back off too much of the initial shifts */
210 ldelta = dmin(ldmax,ldelta);
211 rdelta = dmin(rdmax,rdelta);
212 /* Compute the element growth when shifting to both ends of the cluster */
213 /* accept the shift if there is no element growth at one of the two ends */
216 dplus[1] = d__[1] + s;
217 if (dabs(dplus[1]) < *pivmin) {
218 dplus[1] = -(*pivmin);
219 /* Need to set SAWNAN1 because refined RRR test should not be used */
223 max1 = dabs(dplus[1]);
225 for (i__ = 1; i__ <= i__1; ++i__) {
226 lplus[i__] = ld[i__] / dplus[i__];
227 s = s * lplus[i__] * l[i__] - lsigma;
228 dplus[i__ + 1] = d__[i__ + 1] + s;
229 if ((r__1 = dplus[i__ + 1], dabs(r__1)) < *pivmin) {
230 dplus[i__ + 1] = -(*pivmin);
231 /* Need to set SAWNAN1 because refined RRR test should not be used */
236 r__2 = max1, r__3 = (r__1 = dplus[i__ + 1], dabs(r__1));
237 max1 = dmax(r__2,r__3);
240 sawnan1 = sawnan1 || sisnan_(&max1);
241 if (forcer || max1 <= growthbound && ! sawnan1) {
248 work[1] = d__[1] + s;
249 if (dabs(work[1]) < *pivmin) {
250 work[1] = -(*pivmin);
251 /* Need to set SAWNAN2 because refined RRR test should not be used */
255 max2 = dabs(work[1]);
257 for (i__ = 1; i__ <= i__1; ++i__) {
258 work[*n + i__] = ld[i__] / work[i__];
259 s = s * work[*n + i__] * l[i__] - rsigma;
260 work[i__ + 1] = d__[i__ + 1] + s;
261 if ((r__1 = work[i__ + 1], dabs(r__1)) < *pivmin) {
262 work[i__ + 1] = -(*pivmin);
263 /* Need to set SAWNAN2 because refined RRR test should not be used */
268 r__2 = max2, r__3 = (r__1 = work[i__ + 1], dabs(r__1));
269 max2 = dmax(r__2,r__3);
272 sawnan2 = sawnan2 || sisnan_(&max2);
273 if (forcer || max2 <= growthbound && ! sawnan2) {
278 /* If we are at this point, both shifts led to too much element growth */
279 /* Record the better of the two shifts (provided it didn't lead to NaN) */
280 if (sawnan1 && sawnan2) {
281 /* both MAX1 and MAX2 are NaN */
286 if (max1 <= smlgrowth) {
292 if (sawnan1 || max2 <= max1) {
295 if (max2 <= smlgrowth) {
301 /* If we are here, both the left and the right shift led to */
302 /* element growth. If the element growth is moderate, then */
303 /* we may still accept the representation, if it passes a */
304 /* refined test for RRR. This test supposes that no NaN occurred. */
305 /* Moreover, we use the refined RRR test only for isolated clusters. */
306 if (clwdth < mingap / 128.f && dmin(max1,max2) < fail2 && ! sawnan1 && !
313 if (tryrrr1 && dorrr1) {
315 tmp = (r__1 = dplus[*n], dabs(r__1));
319 for (i__ = *n - 1; i__ >= 1; --i__) {
321 prod = dplus[i__ + 1] * work[*n + i__ + 1] / (dplus[i__] *
322 work[*n + i__]) * oldp;
324 prod *= (r__1 = work[*n + i__], dabs(r__1));
327 /* Computing 2nd power */
331 r__2 = tmp, r__3 = (r__1 = dplus[i__] * prod, dabs(r__1));
332 tmp = dmax(r__2,r__3);
335 rrr1 = tmp / (*spdiam * sqrt(znm2));
341 } else if (indx == 2) {
342 tmp = (r__1 = work[*n], dabs(r__1));
346 for (i__ = *n - 1; i__ >= 1; --i__) {
348 prod = work[i__ + 1] * lplus[i__ + 1] / (work[i__] *
351 prod *= (r__1 = lplus[i__], dabs(r__1));
354 /* Computing 2nd power */
358 r__2 = tmp, r__3 = (r__1 = work[i__] * prod, dabs(r__1));
359 tmp = dmax(r__2,r__3);
362 rrr2 = tmp / (*spdiam * sqrt(znm2));
372 /* If we are here, both shifts failed also the RRR test. */
373 /* Back off to the outside */
375 r__1 = lsigma - ldelta, r__2 = lsigma - ldmax;
376 lsigma = dmax(r__1,r__2);
378 r__1 = rsigma + rdelta, r__2 = rsigma + rdmax;
379 rsigma = dmin(r__1,r__2);
385 /* None of the representations investigated satisfied our */
386 /* criteria. Take the best one we found. */
387 if (smlgrowth < fail || nofail) {
399 } else if (shift == 2) {
400 /* store new L and D back into DPLUS, LPLUS */
401 scopy_(n, &work[1], &c__1, &dplus[1], &c__1);
403 scopy_(&i__1, &work[*n + 1], &c__1, &lplus[1], &c__1);