3 /* Subroutine */ int slar1v_(integer *n, integer *b1, integer *bn, real *
4 lambda, real *d__, real *l, real *ld, real *lld, real *pivmin, real *
5 gaptol, real *z__, logical *wantnc, integer *negcnt, real *ztz, real *
6 mingma, integer *r__, integer *isuppz, real *nrminv, real *resid,
7 real *rqcorr, real *work)
9 /* System generated locals */
11 real r__1, r__2, r__3;
13 /* Builtin functions */
14 double sqrt(doublereal);
21 integer neg1, neg2, indp, inds;
23 extern doublereal slamch_(char *);
24 integer indlpl, indumn;
25 extern logical sisnan_(real *);
27 logical sawnan1, sawnan2;
30 /* -- LAPACK auxiliary routine (version 3.1) -- */
31 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
34 /* .. Scalar Arguments .. */
36 /* .. Array Arguments .. */
42 /* SLAR1V computes the (scaled) r-th column of the inverse of */
43 /* the sumbmatrix in rows B1 through BN of the tridiagonal matrix */
44 /* L D L^T - sigma I. When sigma is close to an eigenvalue, the */
45 /* computed vector is an accurate eigenvector. Usually, r corresponds */
46 /* to the index where the eigenvector is largest in magnitude. */
47 /* The following steps accomplish this computation : */
48 /* (a) Stationary qd transform, L D L^T - sigma I = L(+) D(+) L(+)^T, */
49 /* (b) Progressive qd transform, L D L^T - sigma I = U(-) D(-) U(-)^T, */
50 /* (c) Computation of the diagonal elements of the inverse of */
51 /* L D L^T - sigma I by combining the above transforms, and choosing */
52 /* r as the index where the diagonal of the inverse is (one of the) */
53 /* largest in magnitude. */
54 /* (d) Computation of the (scaled) r-th column of the inverse using the */
55 /* twisted factorization obtained by combining the top part of the */
56 /* the stationary and the bottom part of the progressive transform. */
61 /* N (input) INTEGER */
62 /* The order of the matrix L D L^T. */
64 /* B1 (input) INTEGER */
65 /* First index of the submatrix of L D L^T. */
67 /* BN (input) INTEGER */
68 /* Last index of the submatrix of L D L^T. */
70 /* LAMBDA (input) REAL */
71 /* The shift. In order to compute an accurate eigenvector, */
72 /* LAMBDA should be a good approximation to an eigenvalue */
75 /* L (input) REAL array, dimension (N-1) */
76 /* The (n-1) subdiagonal elements of the unit bidiagonal matrix */
77 /* L, in elements 1 to N-1. */
79 /* D (input) REAL array, dimension (N) */
80 /* The n diagonal elements of the diagonal matrix D. */
82 /* LD (input) REAL array, dimension (N-1) */
83 /* The n-1 elements L(i)*D(i). */
85 /* LLD (input) REAL array, dimension (N-1) */
86 /* The n-1 elements L(i)*L(i)*D(i). */
88 /* PIVMIN (input) REAL */
89 /* The minimum pivot in the Sturm sequence. */
91 /* GAPTOL (input) REAL */
92 /* Tolerance that indicates when eigenvector entries are negligible */
93 /* w.r.t. their contribution to the residual. */
95 /* Z (input/output) REAL array, dimension (N) */
96 /* On input, all entries of Z must be set to 0. */
97 /* On output, Z contains the (scaled) r-th column of the */
98 /* inverse. The scaling is such that Z(R) equals 1. */
100 /* WANTNC (input) LOGICAL */
101 /* Specifies whether NEGCNT has to be computed. */
103 /* NEGCNT (output) INTEGER */
104 /* If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin */
105 /* in the matrix factorization L D L^T, and NEGCNT = -1 otherwise. */
107 /* ZTZ (output) REAL */
108 /* The square of the 2-norm of Z. */
110 /* MINGMA (output) REAL */
111 /* The reciprocal of the largest (in magnitude) diagonal */
112 /* element of the inverse of L D L^T - sigma I. */
114 /* R (input/output) INTEGER */
115 /* The twist index for the twisted factorization used to */
117 /* On input, 0 <= R <= N. If R is input as 0, R is set to */
118 /* the index where (L D L^T - sigma I)^{-1} is largest */
119 /* in magnitude. If 1 <= R <= N, R is unchanged. */
120 /* On output, R contains the twist index used to compute Z. */
121 /* Ideally, R designates the position of the maximum entry in the */
124 /* ISUPPZ (output) INTEGER array, dimension (2) */
125 /* The support of the vector in Z, i.e., the vector Z is */
126 /* nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ). */
128 /* NRMINV (output) REAL */
129 /* NRMINV = 1/SQRT( ZTZ ) */
131 /* RESID (output) REAL */
132 /* The residual of the FP vector. */
133 /* RESID = ABS( MINGMA )/SQRT( ZTZ ) */
135 /* RQCORR (output) REAL */
136 /* The Rayleigh Quotient correction to LAMBDA. */
137 /* RQCORR = MINGMA*TMP */
139 /* WORK (workspace) REAL array, dimension (4*N) */
141 /* Further Details */
142 /* =============== */
144 /* Based on contributions by */
145 /* Beresford Parlett, University of California, Berkeley, USA */
146 /* Jim Demmel, University of California, Berkeley, USA */
147 /* Inderjit Dhillon, University of Texas, Austin, USA */
148 /* Osni Marques, LBNL/NERSC, USA */
149 /* Christof Voemel, University of California, Berkeley, USA */
151 /* ===================================================================== */
153 /* .. Parameters .. */
155 /* .. Local Scalars .. */
157 /* .. External Functions .. */
159 /* .. Intrinsic Functions .. */
161 /* .. Executable Statements .. */
163 /* Parameter adjustments */
173 eps = slamch_("Precision");
181 /* Storage for LPLUS */
183 /* Storage for UMINUS */
185 inds = (*n << 1) + 1;
190 work[inds + *b1 - 1] = lld[*b1 - 1];
193 /* Compute the stationary transform (using the differential form) */
194 /* until the index R2. */
198 s = work[inds + *b1 - 1] - *lambda;
200 for (i__ = *b1; i__ <= i__1; ++i__) {
201 dplus = d__[i__] + s;
202 work[indlpl + i__] = ld[i__] / dplus;
206 work[inds + i__] = s * work[indlpl + i__] * l[i__];
207 s = work[inds + i__] - *lambda;
210 sawnan1 = sisnan_(&s);
215 for (i__ = r1; i__ <= i__1; ++i__) {
216 dplus = d__[i__] + s;
217 work[indlpl + i__] = ld[i__] / dplus;
218 work[inds + i__] = s * work[indlpl + i__] * l[i__];
219 s = work[inds + i__] - *lambda;
222 sawnan1 = sisnan_(&s);
226 /* Runs a slower version of the above loop if a NaN is detected */
228 s = work[inds + *b1 - 1] - *lambda;
230 for (i__ = *b1; i__ <= i__1; ++i__) {
231 dplus = d__[i__] + s;
232 if (dabs(dplus) < *pivmin) {
235 work[indlpl + i__] = ld[i__] / dplus;
239 work[inds + i__] = s * work[indlpl + i__] * l[i__];
240 if (work[indlpl + i__] == 0.f) {
241 work[inds + i__] = lld[i__];
243 s = work[inds + i__] - *lambda;
247 for (i__ = r1; i__ <= i__1; ++i__) {
248 dplus = d__[i__] + s;
249 if (dabs(dplus) < *pivmin) {
252 work[indlpl + i__] = ld[i__] / dplus;
253 work[inds + i__] = s * work[indlpl + i__] * l[i__];
254 if (work[indlpl + i__] == 0.f) {
255 work[inds + i__] = lld[i__];
257 s = work[inds + i__] - *lambda;
262 /* Compute the progressive transform (using the differential form) */
263 /* until the index R1 */
267 work[indp + *bn - 1] = d__[*bn] - *lambda;
269 for (i__ = *bn - 1; i__ >= i__1; --i__) {
270 dminus = lld[i__] + work[indp + i__];
271 tmp = d__[i__] / dminus;
275 work[indumn + i__] = l[i__] * tmp;
276 work[indp + i__ - 1] = work[indp + i__] * tmp - *lambda;
279 tmp = work[indp + r1 - 1];
280 sawnan2 = sisnan_(&tmp);
282 /* Runs a slower version of the above loop if a NaN is detected */
285 for (i__ = *bn - 1; i__ >= i__1; --i__) {
286 dminus = lld[i__] + work[indp + i__];
287 if (dabs(dminus) < *pivmin) {
290 tmp = d__[i__] / dminus;
294 work[indumn + i__] = l[i__] * tmp;
295 work[indp + i__ - 1] = work[indp + i__] * tmp - *lambda;
297 work[indp + i__ - 1] = d__[i__] - *lambda;
303 /* Find the index (from R1 to R2) of the largest (in magnitude) */
304 /* diagonal element of the inverse */
306 *mingma = work[inds + r1 - 1] + work[indp + r1 - 1];
311 *negcnt = neg1 + neg2;
315 if (dabs(*mingma) == 0.f) {
316 *mingma = eps * work[inds + r1 - 1];
320 for (i__ = r1; i__ <= i__1; ++i__) {
321 tmp = work[inds + i__] + work[indp + i__];
323 tmp = eps * work[inds + i__];
325 if (dabs(tmp) <= dabs(*mingma)) {
332 /* Compute the FP vector: solve N^T v = e_r */
339 /* Compute the FP vector upwards from R */
341 if (! sawnan1 && ! sawnan2) {
343 for (i__ = *r__ - 1; i__ >= i__1; --i__) {
344 z__[i__] = -(work[indlpl + i__] * z__[i__ + 1]);
345 if (((r__1 = z__[i__], dabs(r__1)) + (r__2 = z__[i__ + 1], dabs(
346 r__2))) * (r__3 = ld[i__], dabs(r__3)) < *gaptol) {
351 *ztz += z__[i__] * z__[i__];
357 /* Run slower loop if NaN occurred. */
359 for (i__ = *r__ - 1; i__ >= i__1; --i__) {
360 if (z__[i__ + 1] == 0.f) {
361 z__[i__] = -(ld[i__ + 1] / ld[i__]) * z__[i__ + 2];
363 z__[i__] = -(work[indlpl + i__] * z__[i__ + 1]);
365 if (((r__1 = z__[i__], dabs(r__1)) + (r__2 = z__[i__ + 1], dabs(
366 r__2))) * (r__3 = ld[i__], dabs(r__3)) < *gaptol) {
371 *ztz += z__[i__] * z__[i__];
377 /* Compute the FP vector downwards from R in blocks of size BLKSIZ */
378 if (! sawnan1 && ! sawnan2) {
380 for (i__ = *r__; i__ <= i__1; ++i__) {
381 z__[i__ + 1] = -(work[indumn + i__] * z__[i__]);
382 if (((r__1 = z__[i__], dabs(r__1)) + (r__2 = z__[i__ + 1], dabs(
383 r__2))) * (r__3 = ld[i__], dabs(r__3)) < *gaptol) {
388 *ztz += z__[i__ + 1] * z__[i__ + 1];
394 /* Run slower loop if NaN occurred. */
396 for (i__ = *r__; i__ <= i__1; ++i__) {
397 if (z__[i__] == 0.f) {
398 z__[i__ + 1] = -(ld[i__ - 1] / ld[i__]) * z__[i__ - 1];
400 z__[i__ + 1] = -(work[indumn + i__] * z__[i__]);
402 if (((r__1 = z__[i__], dabs(r__1)) + (r__2 = z__[i__ + 1], dabs(
403 r__2))) * (r__3 = ld[i__], dabs(r__3)) < *gaptol) {
408 *ztz += z__[i__ + 1] * z__[i__ + 1];
415 /* Compute quantities for convergence test */
419 *resid = dabs(*mingma) * *nrminv;
420 *rqcorr = *mingma * tmp;