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