1 /* ./src_f77/dlarrv.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static doublereal c_b6 = 0.;
11 static integer c__1 = 1;
12
dlarrv_(integer * n,doublereal * d__,doublereal * l,integer * isplit,integer * m,doublereal * w,integer * iblock,doublereal * gersch,doublereal * tol,doublereal * z__,integer * ldz,integer * isuppz,doublereal * work,integer * iwork,integer * info)13 /* Subroutine */ int dlarrv_(integer *n, doublereal *d__, doublereal *l,
14 integer *isplit, integer *m, doublereal *w, integer *iblock,
15 doublereal *gersch, doublereal *tol, doublereal *z__, integer *ldz,
16 integer *isuppz, doublereal *work, integer *iwork, integer *info)
17 {
18 /* System generated locals */
19 integer z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5, i__6;
20 doublereal d__1, d__2;
21
22 /* Builtin functions */
23 double sqrt(doublereal);
24
25 /* Local variables */
26 static integer i__, j, k, p, q, im, in;
27 static doublereal gap, eps, ztz, tmp1;
28 static integer iend, jblk;
29 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
30 integer *);
31 static integer iter, temp[1], ktot;
32 extern doublereal dnrm2_(integer *, doublereal *, integer *);
33 static integer itmp1, itmp2;
34 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
35 integer *);
36 static integer indld;
37 static doublereal sigma;
38 static integer ndone, iinfo, iindr;
39 static doublereal resid;
40 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
41 doublereal *, integer *);
42 static integer nclus;
43 extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *,
44 integer *, doublereal *, integer *);
45 static integer iindc1, iindc2;
46 extern /* Subroutine */ int dlar1v_(integer *, integer *, integer *,
47 doublereal *, doublereal *, doublereal *, doublereal *,
48 doublereal *, doublereal *, doublereal *, doublereal *,
49 doublereal *, integer *, integer *, doublereal *);
50 static doublereal lambda;
51 extern doublereal dlamch_(char *, ftnlen);
52 static integer ibegin, indgap, indlld;
53 extern /* Subroutine */ int dlarrb_(integer *, doublereal *, doublereal *,
54 doublereal *, doublereal *, integer *, integer *, doublereal *,
55 doublereal *, doublereal *, doublereal *, doublereal *,
56 doublereal *, integer *, integer *);
57 static doublereal mingma;
58 static integer oldien, oldncl;
59 static doublereal relgap;
60 extern /* Subroutine */ int dlarrf_(integer *, doublereal *, doublereal *,
61 doublereal *, doublereal *, integer *, integer *, doublereal *,
62 doublereal *, doublereal *, doublereal *, integer *, integer *),
63 dlaset_(char *, integer *, integer *, doublereal *, doublereal *,
64 doublereal *, integer *, ftnlen);
65 static integer oldcls, ndepth, inderr, iindwk;
66 extern /* Subroutine */ int dstein_(integer *, doublereal *, doublereal *,
67 integer *, doublereal *, integer *, integer *, doublereal *,
68 integer *, doublereal *, integer *, integer *, integer *);
69 static logical mgscls;
70 static integer lsbdpt, newcls, oldfst;
71 static doublereal minrgp;
72 static integer indwrk, oldlst;
73 static doublereal reltol;
74 static integer maxitr, newfrs, newftt;
75 static doublereal mgstol;
76 static integer nsplit;
77 static doublereal nrminv, rqcorr;
78 static integer newlst, newsiz;
79
80
81 /* -- LAPACK auxiliary routine (version 3.0) -- */
82 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
83 /* Courant Institute, Argonne National Lab, and Rice University */
84 /* June 30, 1999 */
85
86 /* .. Scalar Arguments .. */
87 /* .. */
88 /* .. Array Arguments .. */
89 /* .. */
90
91 /* Purpose */
92 /* ======= */
93
94 /* DLARRV computes the eigenvectors of the tridiagonal matrix */
95 /* T = L D L^T given L, D and the eigenvalues of L D L^T. */
96 /* The input eigenvalues should have high relative accuracy with */
97 /* respect to the entries of L and D. The desired accuracy of the */
98 /* output can be specified by the input parameter TOL. */
99
100 /* Arguments */
101 /* ========= */
102
103 /* N (input) INTEGER */
104 /* The order of the matrix. N >= 0. */
105
106 /* D (input/output) DOUBLE PRECISION array, dimension (N) */
107 /* On entry, the n diagonal elements of the diagonal matrix D. */
108 /* On exit, D may be overwritten. */
109
110 /* L (input/output) DOUBLE PRECISION array, dimension (N-1) */
111 /* On entry, the (n-1) subdiagonal elements of the unit */
112 /* bidiagonal matrix L in elements 1 to N-1 of L. L(N) need */
113 /* not be set. On exit, L is overwritten. */
114
115 /* ISPLIT (input) INTEGER array, dimension (N) */
116 /* The splitting points, at which T breaks up into submatrices. */
117 /* The first submatrix consists of rows/columns 1 to */
118 /* ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 */
119 /* through ISPLIT( 2 ), etc. */
120
121 /* TOL (input) DOUBLE PRECISION */
122 /* The absolute error tolerance for the */
123 /* eigenvalues/eigenvectors. */
124 /* Errors in the input eigenvalues must be bounded by TOL. */
125 /* The eigenvectors output have residual norms */
126 /* bounded by TOL, and the dot products between different */
127 /* eigenvectors are bounded by TOL. TOL must be at least */
128 /* N*EPS*|T|, where EPS is the machine precision and |T| is */
129 /* the 1-norm of the tridiagonal matrix. */
130
131 /* M (input) INTEGER */
132 /* The total number of eigenvalues found. 0 <= M <= N. */
133 /* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. */
134
135 /* W (input) DOUBLE PRECISION array, dimension (N) */
136 /* The first M elements of W contain the eigenvalues for */
137 /* which eigenvectors are to be computed. The eigenvalues */
138 /* should be grouped by split-off block and ordered from */
139 /* smallest to largest within the block ( The output array */
140 /* W from DLARRE is expected here ). */
141 /* Errors in W must be bounded by TOL (see above). */
142
143 /* IBLOCK (input) INTEGER array, dimension (N) */
144 /* The submatrix indices associated with the corresponding */
145 /* eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to */
146 /* the first submatrix from the top, =2 if W(i) belongs to */
147 /* the second submatrix, etc. */
148
149 /* Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) ) */
150 /* If JOBZ = 'V', then if INFO = 0, the first M columns of Z */
151 /* contain the orthonormal eigenvectors of the matrix T */
152 /* corresponding to the selected eigenvalues, with the i-th */
153 /* column of Z holding the eigenvector associated with W(i). */
154 /* If JOBZ = 'N', then Z is not referenced. */
155 /* Note: the user must ensure that at least max(1,M) columns are */
156 /* supplied in the array Z; if RANGE = 'V', the exact value of M */
157 /* is not known in advance and an upper bound must be used. */
158
159 /* LDZ (input) INTEGER */
160 /* The leading dimension of the array Z. LDZ >= 1, and if */
161 /* JOBZ = 'V', LDZ >= max(1,N). */
162
163 /* ISUPPZ (output) INTEGER ARRAY, dimension ( 2*max(1,M) ) */
164 /* The support of the eigenvectors in Z, i.e., the indices */
165 /* indicating the nonzero elements in Z. The i-th eigenvector */
166 /* is nonzero only in elements ISUPPZ( 2*i-1 ) through */
167 /* ISUPPZ( 2*i ). */
168
169 /* WORK (workspace) DOUBLE PRECISION array, dimension (13*N) */
170
171 /* IWORK (workspace) INTEGER array, dimension (6*N) */
172
173 /* INFO (output) INTEGER */
174 /* = 0: successful exit */
175 /* < 0: if INFO = -i, the i-th argument had an illegal value */
176 /* > 0: if INFO = 1, internal error in DLARRB */
177 /* if INFO = 2, internal error in DSTEIN */
178
179 /* Further Details */
180 /* =============== */
181
182 /* Based on contributions by */
183 /* Inderjit Dhillon, IBM Almaden, USA */
184 /* Osni Marques, LBNL/NERSC, USA */
185
186 /* ===================================================================== */
187
188 /* .. Parameters .. */
189 /* .. */
190 /* .. Local Scalars .. */
191 /* .. */
192 /* .. External Functions .. */
193 /* .. */
194 /* .. External Subroutines .. */
195 /* .. */
196 /* .. Intrinsic Functions .. */
197 /* .. */
198 /* .. Local Arrays .. */
199 /* .. */
200 /* .. Executable Statements .. */
201
202 /* Test the input parameters. */
203
204 /* Parameter adjustments */
205 --d__;
206 --l;
207 --isplit;
208 --w;
209 --iblock;
210 --gersch;
211 z_dim1 = *ldz;
212 z_offset = 1 + z_dim1;
213 z__ -= z_offset;
214 --isuppz;
215 --work;
216 --iwork;
217
218 /* Function Body */
219 inderr = *n + 1;
220 indld = *n << 1;
221 indlld = *n * 3;
222 indgap = *n << 2;
223 indwrk = *n * 5 + 1;
224
225 iindr = *n;
226 iindc1 = *n << 1;
227 iindc2 = *n * 3;
228 iindwk = (*n << 2) + 1;
229
230 eps = dlamch_("Precision", (ftnlen)9);
231
232 i__1 = *n << 1;
233 for (i__ = 1; i__ <= i__1; ++i__) {
234 iwork[i__] = 0;
235 /* L10: */
236 }
237 i__1 = *m;
238 for (i__ = 1; i__ <= i__1; ++i__) {
239 work[inderr + i__ - 1] = eps * (d__1 = w[i__], abs(d__1));
240 /* L20: */
241 }
242 dlaset_("Full", n, n, &c_b6, &c_b6, &z__[z_offset], ldz, (ftnlen)4);
243 mgstol = eps * 5.;
244
245 nsplit = iblock[*m];
246 ibegin = 1;
247 i__1 = nsplit;
248 for (jblk = 1; jblk <= i__1; ++jblk) {
249 iend = isplit[jblk];
250
251 /* Find the eigenvectors of the submatrix indexed IBEGIN */
252 /* through IEND. */
253
254 if (ibegin == iend) {
255 z__[ibegin + ibegin * z_dim1] = 1.;
256 isuppz[(ibegin << 1) - 1] = ibegin;
257 isuppz[ibegin * 2] = ibegin;
258 ibegin = iend + 1;
259 goto L170;
260 }
261 oldien = ibegin - 1;
262 in = iend - oldien;
263 /* Computing MIN */
264 d__1 = .01, d__2 = 1. / (doublereal) in;
265 reltol = min(d__1,d__2);
266 im = in;
267 dcopy_(&im, &w[ibegin], &c__1, &work[1], &c__1);
268 i__2 = in - 1;
269 for (i__ = 1; i__ <= i__2; ++i__) {
270 work[indgap + i__] = work[i__ + 1] - work[i__];
271 /* L30: */
272 }
273 /* Computing MAX */
274 d__2 = (d__1 = work[in], abs(d__1));
275 work[indgap + in] = max(d__2,eps);
276 ndone = 0;
277
278 ndepth = 0;
279 lsbdpt = 1;
280 nclus = 1;
281 iwork[iindc1 + 1] = 1;
282 iwork[iindc1 + 2] = in;
283
284 /* While( NDONE.LT.IM ) do */
285
286 L40:
287 if (ndone < im) {
288 oldncl = nclus;
289 nclus = 0;
290 lsbdpt = 1 - lsbdpt;
291 i__2 = oldncl;
292 for (i__ = 1; i__ <= i__2; ++i__) {
293 if (lsbdpt == 0) {
294 oldcls = iindc1;
295 newcls = iindc2;
296 } else {
297 oldcls = iindc2;
298 newcls = iindc1;
299 }
300
301 /* If NDEPTH > 1, retrieve the relatively robust */
302 /* representation (RRR) and perform limited bisection */
303 /* (if necessary) to get approximate eigenvalues. */
304
305 j = oldcls + (i__ << 1);
306 oldfst = iwork[j - 1];
307 oldlst = iwork[j];
308 if (ndepth > 0) {
309 j = oldien + oldfst;
310 dcopy_(&in, &z__[ibegin + j * z_dim1], &c__1, &d__[ibegin]
311 , &c__1);
312 dcopy_(&in, &z__[ibegin + (j + 1) * z_dim1], &c__1, &l[
313 ibegin], &c__1);
314 sigma = l[iend];
315 }
316 k = ibegin;
317 i__3 = in - 1;
318 for (j = 1; j <= i__3; ++j) {
319 work[indld + j] = d__[k] * l[k];
320 work[indlld + j] = work[indld + j] * l[k];
321 ++k;
322 /* L50: */
323 }
324 if (ndepth > 0) {
325 dlarrb_(&in, &d__[ibegin], &l[ibegin], &work[indld + 1], &
326 work[indlld + 1], &oldfst, &oldlst, &sigma, &
327 reltol, &work[1], &work[indgap + 1], &work[inderr]
328 , &work[indwrk], &iwork[iindwk], &iinfo);
329 if (iinfo != 0) {
330 *info = 1;
331 return 0;
332 }
333 }
334
335 /* Classify eigenvalues of the current representation (RRR) */
336 /* as (i) isolated, (ii) loosely clustered or (iii) tightly */
337 /* clustered */
338
339 newfrs = oldfst;
340 i__3 = oldlst;
341 for (j = oldfst; j <= i__3; ++j) {
342 if (j == oldlst || work[indgap + j] >= reltol * (d__1 =
343 work[j], abs(d__1))) {
344 newlst = j;
345 } else {
346
347 /* continue (to the next loop) */
348
349 relgap = work[indgap + j] / (d__1 = work[j], abs(d__1)
350 );
351 if (j == newfrs) {
352 minrgp = relgap;
353 } else {
354 minrgp = min(minrgp,relgap);
355 }
356 goto L140;
357 }
358 newsiz = newlst - newfrs + 1;
359 maxitr = 10;
360 newftt = oldien + newfrs;
361 if (newsiz > 1) {
362 mgscls = newsiz <= 20 && minrgp >= mgstol;
363 if (! mgscls) {
364 dlarrf_(&in, &d__[ibegin], &l[ibegin], &work[
365 indld + 1], &work[indlld + 1], &newfrs, &
366 newlst, &work[1], &z__[ibegin + newftt *
367 z_dim1], &z__[ibegin + (newftt + 1) *
368 z_dim1], &work[indwrk], &iwork[iindwk],
369 info);
370 if (*info == 0) {
371 ++nclus;
372 k = newcls + (nclus << 1);
373 iwork[k - 1] = newfrs;
374 iwork[k] = newlst;
375 } else {
376 *info = 0;
377 if (minrgp >= mgstol) {
378 mgscls = TRUE_;
379 } else {
380
381 /* Call DSTEIN to process this tight cluster. */
382 /* This happens only if MINRGP <= MGSTOL */
383 /* and DLARRF returns INFO = 1. The latter */
384 /* means that a new RRR to "break" the */
385 /* cluster could not be found. */
386
387 work[indwrk] = d__[ibegin];
388 i__4 = in - 1;
389 for (k = 1; k <= i__4; ++k) {
390 work[indwrk + k] = d__[ibegin + k] +
391 work[indlld + k];
392 /* L60: */
393 }
394 i__4 = newsiz;
395 for (k = 1; k <= i__4; ++k) {
396 iwork[iindwk + k - 1] = 1;
397 /* L70: */
398 }
399 i__4 = newlst;
400 for (k = newfrs; k <= i__4; ++k) {
401 isuppz[(ibegin + k << 1) - 3] = 1;
402 isuppz[(ibegin + k << 1) - 2] = in;
403 /* L80: */
404 }
405 temp[0] = in;
406 dstein_(&in, &work[indwrk], &work[indld +
407 1], &newsiz, &work[newfrs], &
408 iwork[iindwk], temp, &z__[ibegin
409 + newftt * z_dim1], ldz, &work[
410 indwrk + in], &iwork[iindwk + in],
411 &iwork[iindwk + (in << 1)], &
412 iinfo);
413 if (iinfo != 0) {
414 *info = 2;
415 return 0;
416 }
417 ndone += newsiz;
418 }
419 }
420 }
421 } else {
422 mgscls = FALSE_;
423 }
424 if (newsiz == 1 || mgscls) {
425 ktot = newftt;
426 i__4 = newlst;
427 for (k = newfrs; k <= i__4; ++k) {
428 iter = 0;
429 L90:
430 lambda = work[k];
431 dlar1v_(&in, &c__1, &in, &lambda, &d__[ibegin], &
432 l[ibegin], &work[indld + 1], &work[indlld
433 + 1], &gersch[(oldien << 1) + 1], &z__[
434 ibegin + ktot * z_dim1], &ztz, &mingma, &
435 iwork[iindr + ktot], &isuppz[(ktot << 1)
436 - 1], &work[indwrk]);
437 tmp1 = 1. / ztz;
438 nrminv = sqrt(tmp1);
439 resid = abs(mingma) * nrminv;
440 rqcorr = mingma * tmp1;
441 if (k == in) {
442 gap = work[indgap + k - 1];
443 } else if (k == 1) {
444 gap = work[indgap + k];
445 } else {
446 /* Computing MIN */
447 d__1 = work[indgap + k - 1], d__2 = work[
448 indgap + k];
449 gap = min(d__1,d__2);
450 }
451 ++iter;
452 if (resid > *tol * gap && abs(rqcorr) > eps * 4. *
453 abs(lambda)) {
454 work[k] = lambda + rqcorr;
455 if (iter < maxitr) {
456 goto L90;
457 }
458 }
459 iwork[ktot] = 1;
460 if (newsiz == 1) {
461 ++ndone;
462 }
463 dscal_(&in, &nrminv, &z__[ibegin + ktot * z_dim1],
464 &c__1);
465 ++ktot;
466 /* L100: */
467 }
468 if (newsiz > 1) {
469 itmp1 = isuppz[(newftt << 1) - 1];
470 itmp2 = isuppz[newftt * 2];
471 ktot = oldien + newlst;
472 i__4 = ktot;
473 for (p = newftt + 1; p <= i__4; ++p) {
474 i__5 = p - 1;
475 for (q = newftt; q <= i__5; ++q) {
476 tmp1 = -ddot_(&in, &z__[ibegin + p *
477 z_dim1], &c__1, &z__[ibegin + q *
478 z_dim1], &c__1);
479 daxpy_(&in, &tmp1, &z__[ibegin + q *
480 z_dim1], &c__1, &z__[ibegin + p *
481 z_dim1], &c__1);
482 /* L110: */
483 }
484 tmp1 = 1. / dnrm2_(&in, &z__[ibegin + p *
485 z_dim1], &c__1);
486 dscal_(&in, &tmp1, &z__[ibegin + p * z_dim1],
487 &c__1);
488 /* Computing MIN */
489 i__5 = itmp1, i__6 = isuppz[(p << 1) - 1];
490 itmp1 = min(i__5,i__6);
491 /* Computing MAX */
492 i__5 = itmp2, i__6 = isuppz[p * 2];
493 itmp2 = max(i__5,i__6);
494 /* L120: */
495 }
496 i__4 = ktot;
497 for (p = newftt; p <= i__4; ++p) {
498 isuppz[(p << 1) - 1] = itmp1;
499 isuppz[p * 2] = itmp2;
500 /* L130: */
501 }
502 ndone += newsiz;
503 }
504 }
505 newfrs = j + 1;
506 L140:
507 ;
508 }
509 /* L150: */
510 }
511 ++ndepth;
512 goto L40;
513 }
514 j = ibegin << 1;
515 i__2 = iend;
516 for (i__ = ibegin; i__ <= i__2; ++i__) {
517 isuppz[j - 1] += oldien;
518 isuppz[j] += oldien;
519 j += 2;
520 /* L160: */
521 }
522 ibegin = iend + 1;
523 L170:
524 ;
525 }
526
527 return 0;
528
529 /* End of DLARRV */
530
531 } /* dlarrv_ */
532
533