1 /* ../netlib/clarrv.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib;
2 on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */
3 #include "FLA_f2c.h" /* Table of constant values */
4 static complex c_b1 =
5 {
6 0.f,0.f
7 }
8 ;
9 static integer c__1 = 1;
10 static integer c__2 = 2;
11 static real c_b28 = 0.f;
12 /* > \brief \b CLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenv alues of L D LT. */
13 /* =========== DOCUMENTATION =========== */
14 /* Online html documentation available at */
15 /* http://www.netlib.org/lapack/explore-html/ */
16 /* > \htmlonly */
17 /* > Download CLARRV + dependencies */
18 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarrv. f"> */
19 /* > [TGZ]</a> */
20 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarrv. f"> */
21 /* > [ZIP]</a> */
22 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarrv. f"> */
23 /* > [TXT]</a> */
24 /* > \endhtmlonly */
25 /* Definition: */
26 /* =========== */
27 /* SUBROUTINE CLARRV( N, VL, VU, D, L, PIVMIN, */
28 /* ISPLIT, M, DOL, DOU, MINRGP, */
29 /* RTOL1, RTOL2, W, WERR, WGAP, */
30 /* IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ, */
31 /* WORK, IWORK, INFO ) */
32 /* .. Scalar Arguments .. */
33 /* INTEGER DOL, DOU, INFO, LDZ, M, N */
34 /* REAL MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU */
35 /* .. */
36 /* .. Array Arguments .. */
37 /* INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ), */
38 /* $ ISUPPZ( * ), IWORK( * ) */
39 /* REAL D( * ), GERS( * ), L( * ), W( * ), WERR( * ), */
40 /* $ WGAP( * ), WORK( * ) */
41 /* COMPLEX Z( LDZ, * ) */
42 /* .. */
43 /* > \par Purpose: */
44 /* ============= */
45 /* > */
46 /* > \verbatim */
47 /* > */
48 /* > CLARRV computes the eigenvectors of the tridiagonal matrix */
49 /* > T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T. */
50 /* > The input eigenvalues should have been computed by SLARRE. */
51 /* > \endverbatim */
52 /* Arguments: */
53 /* ========== */
54 /* > \param[in] N */
55 /* > \verbatim */
56 /* > N is INTEGER */
57 /* > The order of the matrix. N >= 0. */
58 /* > \endverbatim */
59 /* > */
60 /* > \param[in] VL */
61 /* > \verbatim */
62 /* > VL is REAL */
63 /* > \endverbatim */
64 /* > */
65 /* > \param[in] VU */
66 /* > \verbatim */
67 /* > VU is REAL */
68 /* > Lower and upper bounds of the interval that contains the desired */
69 /* > eigenvalues. VL < VU. Needed to compute gaps on the left or right */
70 /* > end of the extremal eigenvalues in the desired RANGE. */
71 /* > \endverbatim */
72 /* > */
73 /* > \param[in,out] D */
74 /* > \verbatim */
75 /* > D is REAL array, dimension (N) */
76 /* > On entry, the N diagonal elements of the diagonal matrix D. */
77 /* > On exit, D may be overwritten. */
78 /* > \endverbatim */
79 /* > */
80 /* > \param[in,out] L */
81 /* > \verbatim */
82 /* > L is REAL array, dimension (N) */
83 /* > On entry, the (N-1) subdiagonal elements of the unit */
84 /* > bidiagonal matrix L are in elements 1 to N-1 of L */
85 /* > (if the matrix is not splitted.) At the end of each block */
86 /* > is stored the corresponding shift as given by SLARRE. */
87 /* > On exit, L is overwritten. */
88 /* > \endverbatim */
89 /* > */
90 /* > \param[in] PIVMIN */
91 /* > \verbatim */
92 /* > PIVMIN is REAL */
93 /* > The minimum pivot allowed in the Sturm sequence. */
94 /* > \endverbatim */
95 /* > */
96 /* > \param[in] ISPLIT */
97 /* > \verbatim */
98 /* > ISPLIT is INTEGER array, dimension (N) */
99 /* > The splitting points, at which T breaks up into blocks. */
100 /* > The first block consists of rows/columns 1 to */
101 /* > ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 */
102 /* > through ISPLIT( 2 ), etc. */
103 /* > \endverbatim */
104 /* > */
105 /* > \param[in] M */
106 /* > \verbatim */
107 /* > M is INTEGER */
108 /* > The total number of input eigenvalues. 0 <= M <= N. */
109 /* > \endverbatim */
110 /* > */
111 /* > \param[in] DOL */
112 /* > \verbatim */
113 /* > DOL is INTEGER */
114 /* > \endverbatim */
115 /* > */
116 /* > \param[in] DOU */
117 /* > \verbatim */
118 /* > DOU is INTEGER */
119 /* > If the user wants to compute only selected eigenvectors from all */
120 /* > the eigenvalues supplied, he can specify an index range DOL:DOU. */
121 /* > Or else the setting DOL=1, DOU=M should be applied. */
122 /* > Note that DOL and DOU refer to the order in which the eigenvalues */
123 /* > are stored in W. */
124 /* > If the user wants to compute only selected eigenpairs, then */
125 /* > the columns DOL-1 to DOU+1 of the eigenvector space Z contain the */
126 /* > computed eigenvectors. All other columns of Z are set to zero. */
127 /* > \endverbatim */
128 /* > */
129 /* > \param[in] MINRGP */
130 /* > \verbatim */
131 /* > MINRGP is REAL */
132 /* > \endverbatim */
133 /* > */
134 /* > \param[in] RTOL1 */
135 /* > \verbatim */
136 /* > RTOL1 is REAL */
137 /* > \endverbatim */
138 /* > */
139 /* > \param[in] RTOL2 */
140 /* > \verbatim */
141 /* > RTOL2 is REAL */
142 /* > Parameters for bisection. */
143 /* > An interval [LEFT,RIGHT] has converged if */
144 /* > RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) */
145 /* > \endverbatim */
146 /* > */
147 /* > \param[in,out] W */
148 /* > \verbatim */
149 /* > W is REAL array, dimension (N) */
150 /* > The first M elements of W contain the APPROXIMATE eigenvalues for */
151 /* > which eigenvectors are to be computed. The eigenvalues */
152 /* > should be grouped by split-off block and ordered from */
153 /* > smallest to largest within the block ( The output array */
154 /* > W from SLARRE is expected here ). Furthermore, they are with */
155 /* > respect to the shift of the corresponding root representation */
156 /* > for their block. On exit, W holds the eigenvalues of the */
157 /* > UNshifted matrix. */
158 /* > \endverbatim */
159 /* > */
160 /* > \param[in,out] WERR */
161 /* > \verbatim */
162 /* > WERR is REAL array, dimension (N) */
163 /* > The first M elements contain the semiwidth of the uncertainty */
164 /* > interval of the corresponding eigenvalue in W */
165 /* > \endverbatim */
166 /* > */
167 /* > \param[in,out] WGAP */
168 /* > \verbatim */
169 /* > WGAP is REAL array, dimension (N) */
170 /* > The separation from the right neighbor eigenvalue in W. */
171 /* > \endverbatim */
172 /* > */
173 /* > \param[in] IBLOCK */
174 /* > \verbatim */
175 /* > IBLOCK is INTEGER array, dimension (N) */
176 /* > The indices of the blocks (submatrices) associated with the */
177 /* > corresponding eigenvalues in W;
178 IBLOCK(i)=1 if eigenvalue */
179 /* > W(i) belongs to the first block from the top, =2 if W(i) */
180 /* > belongs to the second block, etc. */
181 /* > \endverbatim */
182 /* > */
183 /* > \param[in] INDEXW */
184 /* > \verbatim */
185 /* > INDEXW is INTEGER array, dimension (N) */
186 /* > The indices of the eigenvalues within each block (submatrix);
187 */
188 /* > for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the */
189 /* > i-th eigenvalue W(i) is the 10-th eigenvalue in the second block. */
190 /* > \endverbatim */
191 /* > */
192 /* > \param[in] GERS */
193 /* > \verbatim */
194 /* > GERS is REAL array, dimension (2*N) */
195 /* > The N Gerschgorin intervals (the i-th Gerschgorin interval */
196 /* > is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should */
197 /* > be computed from the original UNshifted matrix. */
198 /* > \endverbatim */
199 /* > */
200 /* > \param[out] Z */
201 /* > \verbatim */
202 /* > Z is array, dimension (LDZ, max(1,M) ) */
203 /* > If INFO = 0, the first M columns of Z contain the */
204 /* > orthonormal eigenvectors of the matrix T */
205 /* > corresponding to the input eigenvalues, with the i-th */
206 /* > column of Z holding the eigenvector associated with W(i). */
207 /* > Note: the user must ensure that at least max(1,M) columns are */
208 /* > supplied in the array Z. */
209 /* > \endverbatim */
210 /* > */
211 /* > \param[in] LDZ */
212 /* > \verbatim */
213 /* > LDZ is INTEGER */
214 /* > The leading dimension of the array Z. LDZ >= 1, and if */
215 /* > JOBZ = 'V', LDZ >= max(1,N). */
216 /* > \endverbatim */
217 /* > */
218 /* > \param[out] ISUPPZ */
219 /* > \verbatim */
220 /* > ISUPPZ is INTEGER array, dimension ( 2*max(1,M) ) */
221 /* > The support of the eigenvectors in Z, i.e., the indices */
222 /* > indicating the nonzero elements in Z. The I-th eigenvector */
223 /* > is nonzero only in elements ISUPPZ( 2*I-1 ) through */
224 /* > ISUPPZ( 2*I ). */
225 /* > \endverbatim */
226 /* > */
227 /* > \param[out] WORK */
228 /* > \verbatim */
229 /* > WORK is REAL array, dimension (12*N) */
230 /* > \endverbatim */
231 /* > */
232 /* > \param[out] IWORK */
233 /* > \verbatim */
234 /* > IWORK is INTEGER array, dimension (7*N) */
235 /* > \endverbatim */
236 /* > */
237 /* > \param[out] INFO */
238 /* > \verbatim */
239 /* > INFO is INTEGER */
240 /* > = 0: successful exit */
241 /* > */
242 /* > > 0: A problem occured in CLARRV. */
243 /* > < 0: One of the called subroutines signaled an internal problem. */
244 /* > Needs inspection of the corresponding parameter IINFO */
245 /* > for further information. */
246 /* > */
247 /* > =-1: Problem in SLARRB when refining a child's eigenvalues. */
248 /* > =-2: Problem in SLARRF when computing the RRR of a child. */
249 /* > When a child is inside a tight cluster, it can be difficult */
250 /* > to find an RRR. A partial remedy from the user's point of */
251 /* > view is to make the parameter MINRGP smaller and recompile. */
252 /* > However, as the orthogonality of the computed vectors is */
253 /* > proportional to 1/MINRGP, the user should be aware that */
254 /* > he might be trading in precision when he decreases MINRGP. */
255 /* > =-3: Problem in SLARRB when refining a single eigenvalue */
256 /* > after the Rayleigh correction was rejected. */
257 /* > = 5: The Rayleigh Quotient Iteration failed to converge to */
258 /* > full accuracy in MAXITR steps. */
259 /* > \endverbatim */
260 /* Authors: */
261 /* ======== */
262 /* > \author Univ. of Tennessee */
263 /* > \author Univ. of California Berkeley */
264 /* > \author Univ. of Colorado Denver */
265 /* > \author NAG Ltd. */
266 /* > \date September 2012 */
267 /* > \ingroup complexOTHERauxiliary */
268 /* > \par Contributors: */
269 /* ================== */
270 /* > */
271 /* > Beresford Parlett, University of California, Berkeley, USA \n */
272 /* > Jim Demmel, University of California, Berkeley, USA \n */
273 /* > Inderjit Dhillon, University of Texas, Austin, USA \n */
274 /* > Osni Marques, LBNL/NERSC, USA \n */
275 /* > Christof Voemel, University of California, Berkeley, USA */
276 /* ===================================================================== */
277 /* Subroutine */
clarrv_(integer * n,real * vl,real * vu,real * d__,real * l,real * pivmin,integer * isplit,integer * m,integer * dol,integer * dou,real * minrgp,real * rtol1,real * rtol2,real * w,real * werr,real * wgap,integer * iblock,integer * indexw,real * gers,complex * z__,integer * ldz,integer * isuppz,real * work,integer * iwork,integer * info)278 int clarrv_(integer *n, real *vl, real *vu, real *d__, real * l, real *pivmin, integer *isplit, integer *m, integer *dol, integer * dou, real *minrgp, real *rtol1, real *rtol2, real *w, real *werr, real *wgap, integer *iblock, integer *indexw, real *gers, complex * z__, integer *ldz, integer *isuppz, real *work, integer *iwork, integer *info)
279 {
280 /* System generated locals */
281 integer z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5, i__6;
282 real r__1, r__2;
283 complex q__1;
284 logical L__1;
285 /* Builtin functions */
286 double log(doublereal);
287 /* Local variables */
288 integer minwsize, i__, j, k, p, q, miniwsize, ii;
289 real gl;
290 integer im, in;
291 real gu, gap, eps, tau, tol, tmp;
292 integer zto;
293 real ztz;
294 integer iend, jblk;
295 real lgap;
296 integer done;
297 real rgap, left;
298 integer wend, iter;
299 real bstw;
300 integer itmp1, indld;
301 real fudge;
302 integer idone;
303 real sigma;
304 integer iinfo, iindr;
305 real resid;
306 logical eskip;
307 real right;
308 integer nclus, zfrom;
309 extern /* Subroutine */
310 int scopy_(integer *, real *, integer *, real *, integer *);
311 real rqtol;
312 integer iindc1, iindc2, indin1, indin2;
313 extern /* Subroutine */
314 int clar1v_(integer *, integer *, integer *, real *, real *, real *, real *, real *, real *, real *, complex *, logical *, integer *, real *, real *, integer *, integer *, real * , real *, real *, real *);
315 logical stp2ii;
316 real lambda;
317 integer ibegin, indeig;
318 logical needbs;
319 integer indlld;
320 real sgndef, mingma;
321 extern real slamch_(char *);
322 integer oldien, oldncl, wbegin;
323 real spdiam;
324 integer negcnt;
325 extern /* Subroutine */
326 int claset_(char *, integer *, integer *, complex *, complex *, complex *, integer *);
327 integer oldcls;
328 real savgap;
329 integer ndepth;
330 real ssigma;
331 extern /* Subroutine */
332 int csscal_(integer *, real *, complex *, integer *);
333 logical usedbs;
334 integer iindwk, offset;
335 real gaptol;
336 extern /* Subroutine */
337 int slarrb_(integer *, real *, real *, integer *, integer *, real *, real *, integer *, real *, real *, real *, real *, integer *, real *, real *, integer *, integer *);
338 integer newcls, oldfst, indwrk, windex, oldlst;
339 logical usedrq;
340 integer newfst, newftt, parity, windmn, windpl, isupmn, newlst, zusedl;
341 real bstres;
342 integer newsiz, zusedu, zusedw;
343 real nrminv, rqcorr;
344 logical tryrqc;
345 integer isupmx;
346 extern /* Subroutine */
347 int slarrf_(integer *, real *, real *, real *, integer *, integer *, real *, real *, real *, real *, real *, real *, real *, real *, real *, real *, real *, integer *);
348 /* -- LAPACK auxiliary routine (version 3.4.2) -- */
349 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
350 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
351 /* September 2012 */
352 /* .. Scalar Arguments .. */
353 /* .. */
354 /* .. Array Arguments .. */
355 /* .. */
356 /* ===================================================================== */
357 /* .. Parameters .. */
358 /* .. */
359 /* .. Local Scalars .. */
360 /* .. */
361 /* .. External Functions .. */
362 /* .. */
363 /* .. External Subroutines .. */
364 /* .. */
365 /* .. Intrinsic Functions .. */
366 /* .. */
367 /* .. Executable Statements .. */
368 /* .. */
369 /* The first N entries of WORK are reserved for the eigenvalues */
370 /* Parameter adjustments */
371 --d__;
372 --l;
373 --isplit;
374 --w;
375 --werr;
376 --wgap;
377 --iblock;
378 --indexw;
379 --gers;
380 z_dim1 = *ldz;
381 z_offset = 1 + z_dim1;
382 z__ -= z_offset;
383 --isuppz;
384 --work;
385 --iwork;
386 /* Function Body */
387 indld = *n + 1;
388 indlld = (*n << 1) + 1;
389 indin1 = *n * 3 + 1;
390 indin2 = (*n << 2) + 1;
391 indwrk = *n * 5 + 1;
392 minwsize = *n * 12;
393 i__1 = minwsize;
394 for (i__ = 1;
395 i__ <= i__1;
396 ++i__)
397 {
398 work[i__] = 0.f;
399 /* L5: */
400 }
401 /* IWORK(IINDR+1:IINDR+N) hold the twist indices R for the */
402 /* factorization used to compute the FP vector */
403 iindr = 0;
404 /* IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current */
405 /* layer and the one above. */
406 iindc1 = *n;
407 iindc2 = *n << 1;
408 iindwk = *n * 3 + 1;
409 miniwsize = *n * 7;
410 i__1 = miniwsize;
411 for (i__ = 1;
412 i__ <= i__1;
413 ++i__)
414 {
415 iwork[i__] = 0;
416 /* L10: */
417 }
418 zusedl = 1;
419 if (*dol > 1)
420 {
421 /* Set lower bound for use of Z */
422 zusedl = *dol - 1;
423 }
424 zusedu = *m;
425 if (*dou < *m)
426 {
427 /* Set lower bound for use of Z */
428 zusedu = *dou + 1;
429 }
430 /* The width of the part of Z that is used */
431 zusedw = zusedu - zusedl + 1;
432 claset_("Full", n, &zusedw, &c_b1, &c_b1, &z__[zusedl * z_dim1 + 1], ldz);
433 eps = slamch_("Precision");
434 rqtol = eps * 2.f;
435 /* Set expert flags for standard code. */
436 tryrqc = TRUE_;
437 if (*dol == 1 && *dou == *m)
438 {
439 }
440 else
441 {
442 /* Only selected eigenpairs are computed. Since the other evalues */
443 /* are not refined by RQ iteration, bisection has to compute to full */
444 /* accuracy. */
445 *rtol1 = eps * 4.f;
446 *rtol2 = eps * 4.f;
447 }
448 /* The entries WBEGIN:WEND in W, WERR, WGAP correspond to the */
449 /* desired eigenvalues. The support of the nonzero eigenvector */
450 /* entries is contained in the interval IBEGIN:IEND. */
451 /* Remark that if k eigenpairs are desired, then the eigenvectors */
452 /* are stored in k contiguous columns of Z. */
453 /* DONE is the number of eigenvectors already computed */
454 done = 0;
455 ibegin = 1;
456 wbegin = 1;
457 i__1 = iblock[*m];
458 for (jblk = 1;
459 jblk <= i__1;
460 ++jblk)
461 {
462 iend = isplit[jblk];
463 sigma = l[iend];
464 /* Find the eigenvectors of the submatrix indexed IBEGIN */
465 /* through IEND. */
466 wend = wbegin - 1;
467 L15:
468 if (wend < *m)
469 {
470 if (iblock[wend + 1] == jblk)
471 {
472 ++wend;
473 goto L15;
474 }
475 }
476 if (wend < wbegin)
477 {
478 ibegin = iend + 1;
479 goto L170;
480 }
481 else if (wend < *dol || wbegin > *dou)
482 {
483 ibegin = iend + 1;
484 wbegin = wend + 1;
485 goto L170;
486 }
487 /* Find local spectral diameter of the block */
488 gl = gers[(ibegin << 1) - 1];
489 gu = gers[ibegin * 2];
490 i__2 = iend;
491 for (i__ = ibegin + 1;
492 i__ <= i__2;
493 ++i__)
494 {
495 /* Computing MIN */
496 r__1 = gers[(i__ << 1) - 1];
497 gl = min(r__1,gl);
498 /* Computing MAX */
499 r__1 = gers[i__ * 2];
500 gu = max(r__1,gu);
501 /* L20: */
502 }
503 spdiam = gu - gl;
504 /* OLDIEN is the last index of the previous block */
505 oldien = ibegin - 1;
506 /* Calculate the size of the current block */
507 in = iend - ibegin + 1;
508 /* The number of eigenvalues in the current block */
509 im = wend - wbegin + 1;
510 /* This is for a 1x1 block */
511 if (ibegin == iend)
512 {
513 ++done;
514 i__2 = ibegin + wbegin * z_dim1;
515 z__[i__2].r = 1.f;
516 z__[i__2].i = 0.f; // , expr subst
517 isuppz[(wbegin << 1) - 1] = ibegin;
518 isuppz[wbegin * 2] = ibegin;
519 w[wbegin] += sigma;
520 work[wbegin] = w[wbegin];
521 ibegin = iend + 1;
522 ++wbegin;
523 goto L170;
524 }
525 /* The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND) */
526 /* Note that these can be approximations, in this case, the corresp. */
527 /* entries of WERR give the size of the uncertainty interval. */
528 /* The eigenvalue approximations will be refined when necessary as */
529 /* high relative accuracy is required for the computation of the */
530 /* corresponding eigenvectors. */
531 scopy_(&im, &w[wbegin], &c__1, &work[wbegin], &c__1);
532 /* We store in W the eigenvalue approximations w.r.t. the original */
533 /* matrix T. */
534 i__2 = im;
535 for (i__ = 1;
536 i__ <= i__2;
537 ++i__)
538 {
539 w[wbegin + i__ - 1] += sigma;
540 /* L30: */
541 }
542 /* NDEPTH is the current depth of the representation tree */
543 ndepth = 0;
544 /* PARITY is either 1 or 0 */
545 parity = 1;
546 /* NCLUS is the number of clusters for the next level of the */
547 /* representation tree, we start with NCLUS = 1 for the root */
548 nclus = 1;
549 iwork[iindc1 + 1] = 1;
550 iwork[iindc1 + 2] = im;
551 /* IDONE is the number of eigenvectors already computed in the current */
552 /* block */
553 idone = 0;
554 /* loop while( IDONE.LT.IM ) */
555 /* generate the representation tree for the current block and */
556 /* compute the eigenvectors */
557 L40:
558 if (idone < im)
559 {
560 /* This is a crude protection against infinitely deep trees */
561 if (ndepth > *m)
562 {
563 *info = -2;
564 return 0;
565 }
566 /* breadth first processing of the current level of the representation */
567 /* tree: OLDNCL = number of clusters on current level */
568 oldncl = nclus;
569 /* reset NCLUS to count the number of child clusters */
570 nclus = 0;
571 parity = 1 - parity;
572 if (parity == 0)
573 {
574 oldcls = iindc1;
575 newcls = iindc2;
576 }
577 else
578 {
579 oldcls = iindc2;
580 newcls = iindc1;
581 }
582 /* Process the clusters on the current level */
583 i__2 = oldncl;
584 for (i__ = 1;
585 i__ <= i__2;
586 ++i__)
587 {
588 j = oldcls + (i__ << 1);
589 /* OLDFST, OLDLST = first, last index of current cluster. */
590 /* cluster indices start with 1 and are relative */
591 /* to WBEGIN when accessing W, WGAP, WERR, Z */
592 oldfst = iwork[j - 1];
593 oldlst = iwork[j];
594 if (ndepth > 0)
595 {
596 /* Retrieve relatively robust representation (RRR) of cluster */
597 /* that has been computed at the previous level */
598 /* The RRR is stored in Z and overwritten once the eigenvectors */
599 /* have been computed or when the cluster is refined */
600 if (*dol == 1 && *dou == *m)
601 {
602 /* Get representation from location of the leftmost evalue */
603 /* of the cluster */
604 j = wbegin + oldfst - 1;
605 }
606 else
607 {
608 if (wbegin + oldfst - 1 < *dol)
609 {
610 /* Get representation from the left end of Z array */
611 j = *dol - 1;
612 }
613 else if (wbegin + oldfst - 1 > *dou)
614 {
615 /* Get representation from the right end of Z array */
616 j = *dou;
617 }
618 else
619 {
620 j = wbegin + oldfst - 1;
621 }
622 }
623 i__3 = in - 1;
624 for (k = 1;
625 k <= i__3;
626 ++k)
627 {
628 i__4 = ibegin + k - 1 + j * z_dim1;
629 d__[ibegin + k - 1] = z__[i__4].r;
630 i__4 = ibegin + k - 1 + (j + 1) * z_dim1;
631 l[ibegin + k - 1] = z__[i__4].r;
632 /* L45: */
633 }
634 i__3 = iend + j * z_dim1;
635 d__[iend] = z__[i__3].r;
636 i__3 = iend + (j + 1) * z_dim1;
637 sigma = z__[i__3].r;
638 /* Set the corresponding entries in Z to zero */
639 claset_("Full", &in, &c__2, &c_b1, &c_b1, &z__[ibegin + j * z_dim1], ldz);
640 }
641 /* Compute DL and DLL of current RRR */
642 i__3 = iend - 1;
643 for (j = ibegin;
644 j <= i__3;
645 ++j)
646 {
647 tmp = d__[j] * l[j];
648 work[indld - 1 + j] = tmp;
649 work[indlld - 1 + j] = tmp * l[j];
650 /* L50: */
651 }
652 if (ndepth > 0)
653 {
654 /* P and Q are index of the first and last eigenvalue to compute */
655 /* within the current block */
656 p = indexw[wbegin - 1 + oldfst];
657 q = indexw[wbegin - 1 + oldlst];
658 /* Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET */
659 /* through the Q-OFFSET elements of these arrays are to be used. */
660 /* OFFSET = P-OLDFST */
661 offset = indexw[wbegin] - 1;
662 /* perform limited bisection (if necessary) to get approximate */
663 /* eigenvalues to the precision needed. */
664 slarrb_(&in, &d__[ibegin], &work[indlld + ibegin - 1], &p, &q, rtol1, rtol2, &offset, &work[wbegin], &wgap[ wbegin], &werr[wbegin], &work[indwrk], &iwork[ iindwk], pivmin, &spdiam, &in, &iinfo);
665 if (iinfo != 0)
666 {
667 *info = -1;
668 return 0;
669 }
670 /* We also recompute the extremal gaps. W holds all eigenvalues */
671 /* of the unshifted matrix and must be used for computation */
672 /* of WGAP, the entries of WORK might stem from RRRs with */
673 /* different shifts. The gaps from WBEGIN-1+OLDFST to */
674 /* WBEGIN-1+OLDLST are correctly computed in SLARRB. */
675 /* However, we only allow the gaps to become greater since */
676 /* this is what should happen when we decrease WERR */
677 if (oldfst > 1)
678 {
679 /* Computing MAX */
680 r__1 = wgap[wbegin + oldfst - 2];
681 r__2 = w[wbegin + oldfst - 1] - werr[wbegin + oldfst - 1] - w[ wbegin + oldfst - 2] - werr[wbegin + oldfst - 2]; // , expr subst
682 wgap[wbegin + oldfst - 2] = max(r__1,r__2);
683 }
684 if (wbegin + oldlst - 1 < wend)
685 {
686 /* Computing MAX */
687 r__1 = wgap[wbegin + oldlst - 1];
688 r__2 = w[wbegin + oldlst] - werr[wbegin + oldlst] - w[wbegin + oldlst - 1] - werr[wbegin + oldlst - 1]; // , expr subst
689 wgap[wbegin + oldlst - 1] = max(r__1,r__2);
690 }
691 /* Each time the eigenvalues in WORK get refined, we store */
692 /* the newly found approximation with all shifts applied in W */
693 i__3 = oldlst;
694 for (j = oldfst;
695 j <= i__3;
696 ++j)
697 {
698 w[wbegin + j - 1] = work[wbegin + j - 1] + sigma;
699 /* L53: */
700 }
701 }
702 /* Process the current node. */
703 newfst = oldfst;
704 i__3 = oldlst;
705 for (j = oldfst;
706 j <= i__3;
707 ++j)
708 {
709 if (j == oldlst)
710 {
711 /* we are at the right end of the cluster, this is also the */
712 /* boundary of the child cluster */
713 newlst = j;
714 }
715 else if (wgap[wbegin + j - 1] >= *minrgp * (r__1 = work[ wbegin + j - 1], f2c_abs(r__1)))
716 {
717 /* the right relative gap is big enough, the child cluster */
718 /* (NEWFST,..,NEWLST) is well separated from the following */
719 newlst = j;
720 }
721 else
722 {
723 /* inside a child cluster, the relative gap is not */
724 /* big enough. */
725 goto L140;
726 }
727 /* Compute size of child cluster found */
728 newsiz = newlst - newfst + 1;
729 /* NEWFTT is the place in Z where the new RRR or the computed */
730 /* eigenvector is to be stored */
731 if (*dol == 1 && *dou == *m)
732 {
733 /* Store representation at location of the leftmost evalue */
734 /* of the cluster */
735 newftt = wbegin + newfst - 1;
736 }
737 else
738 {
739 if (wbegin + newfst - 1 < *dol)
740 {
741 /* Store representation at the left end of Z array */
742 newftt = *dol - 1;
743 }
744 else if (wbegin + newfst - 1 > *dou)
745 {
746 /* Store representation at the right end of Z array */
747 newftt = *dou;
748 }
749 else
750 {
751 newftt = wbegin + newfst - 1;
752 }
753 }
754 if (newsiz > 1)
755 {
756 /* Current child is not a singleton but a cluster. */
757 /* Compute and store new representation of child. */
758 /* Compute left and right cluster gap. */
759 /* LGAP and RGAP are not computed from WORK because */
760 /* the eigenvalue approximations may stem from RRRs */
761 /* different shifts. However, W hold all eigenvalues */
762 /* of the unshifted matrix. Still, the entries in WGAP */
763 /* have to be computed from WORK since the entries */
764 /* in W might be of the same order so that gaps are not */
765 /* exhibited correctly for very close eigenvalues. */
766 if (newfst == 1)
767 {
768 /* Computing MAX */
769 r__1 = 0.f;
770 r__2 = w[wbegin] - werr[wbegin] - *vl; // , expr subst
771 lgap = max(r__1,r__2);
772 }
773 else
774 {
775 lgap = wgap[wbegin + newfst - 2];
776 }
777 rgap = wgap[wbegin + newlst - 1];
778 /* Compute left- and rightmost eigenvalue of child */
779 /* to high precision in order to shift as close */
780 /* as possible and obtain as large relative gaps */
781 /* as possible */
782 for (k = 1;
783 k <= 2;
784 ++k)
785 {
786 if (k == 1)
787 {
788 p = indexw[wbegin - 1 + newfst];
789 }
790 else
791 {
792 p = indexw[wbegin - 1 + newlst];
793 }
794 offset = indexw[wbegin] - 1;
795 slarrb_(&in, &d__[ibegin], &work[indlld + ibegin - 1], &p, &p, &rqtol, &rqtol, &offset, & work[wbegin], &wgap[wbegin], &werr[wbegin] , &work[indwrk], &iwork[iindwk], pivmin, & spdiam, &in, &iinfo);
796 /* L55: */
797 }
798 if (wbegin + newlst - 1 < *dol || wbegin + newfst - 1 > *dou)
799 {
800 /* if the cluster contains no desired eigenvalues */
801 /* skip the computation of that branch of the rep. tree */
802 /* We could skip before the refinement of the extremal */
803 /* eigenvalues of the child, but then the representation */
804 /* tree could be different from the one when nothing is */
805 /* skipped. For this reason we skip at this place. */
806 idone = idone + newlst - newfst + 1;
807 goto L139;
808 }
809 /* Compute RRR of child cluster. */
810 /* Note that the new RRR is stored in Z */
811 /* SLARRF needs LWORK = 2*N */
812 slarrf_(&in, &d__[ibegin], &l[ibegin], &work[indld + ibegin - 1], &newfst, &newlst, &work[wbegin], &wgap[wbegin], &werr[wbegin], &spdiam, &lgap, &rgap, pivmin, &tau, &work[indin1], &work[ indin2], &work[indwrk], &iinfo);
813 /* In the complex case, SLARRF cannot write */
814 /* the new RRR directly into Z and needs an intermediate */
815 /* workspace */
816 i__4 = in - 1;
817 for (k = 1;
818 k <= i__4;
819 ++k)
820 {
821 i__5 = ibegin + k - 1 + newftt * z_dim1;
822 i__6 = indin1 + k - 1;
823 q__1.r = work[i__6];
824 q__1.i = 0.f; // , expr subst
825 z__[i__5].r = q__1.r;
826 z__[i__5].i = q__1.i; // , expr subst
827 i__5 = ibegin + k - 1 + (newftt + 1) * z_dim1;
828 i__6 = indin2 + k - 1;
829 q__1.r = work[i__6];
830 q__1.i = 0.f; // , expr subst
831 z__[i__5].r = q__1.r;
832 z__[i__5].i = q__1.i; // , expr subst
833 /* L56: */
834 }
835 i__4 = iend + newftt * z_dim1;
836 i__5 = indin1 + in - 1;
837 q__1.r = work[i__5];
838 q__1.i = 0.f; // , expr subst
839 z__[i__4].r = q__1.r;
840 z__[i__4].i = q__1.i; // , expr subst
841 if (iinfo == 0)
842 {
843 /* a new RRR for the cluster was found by SLARRF */
844 /* update shift and store it */
845 ssigma = sigma + tau;
846 i__4 = iend + (newftt + 1) * z_dim1;
847 q__1.r = ssigma;
848 q__1.i = 0.f; // , expr subst
849 z__[i__4].r = q__1.r;
850 z__[i__4].i = q__1.i; // , expr subst
851 /* WORK() are the midpoints and WERR() the semi-width */
852 /* Note that the entries in W are unchanged. */
853 i__4 = newlst;
854 for (k = newfst;
855 k <= i__4;
856 ++k)
857 {
858 fudge = eps * 3.f * (r__1 = work[wbegin + k - 1], f2c_abs(r__1));
859 work[wbegin + k - 1] -= tau;
860 fudge += eps * 4.f * (r__1 = work[wbegin + k - 1], f2c_abs(r__1));
861 /* Fudge errors */
862 werr[wbegin + k - 1] += fudge;
863 /* Gaps are not fudged. Provided that WERR is small */
864 /* when eigenvalues are close, a zero gap indicates */
865 /* that a new representation is needed for resolving */
866 /* the cluster. A fudge could lead to a wrong decision */
867 /* of judging eigenvalues 'separated' which in */
868 /* reality are not. This could have a negative impact */
869 /* on the orthogonality of the computed eigenvectors. */
870 /* L116: */
871 }
872 ++nclus;
873 k = newcls + (nclus << 1);
874 iwork[k - 1] = newfst;
875 iwork[k] = newlst;
876 }
877 else
878 {
879 *info = -2;
880 return 0;
881 }
882 }
883 else
884 {
885 /* Compute eigenvector of singleton */
886 iter = 0;
887 tol = log((real) in) * 4.f * eps;
888 k = newfst;
889 windex = wbegin + k - 1;
890 /* Computing MAX */
891 i__4 = windex - 1;
892 windmn = max(i__4,1);
893 /* Computing MIN */
894 i__4 = windex + 1;
895 windpl = min(i__4,*m);
896 lambda = work[windex];
897 ++done;
898 /* Check if eigenvector computation is to be skipped */
899 if (windex < *dol || windex > *dou)
900 {
901 eskip = TRUE_;
902 goto L125;
903 }
904 else
905 {
906 eskip = FALSE_;
907 }
908 left = work[windex] - werr[windex];
909 right = work[windex] + werr[windex];
910 indeig = indexw[windex];
911 /* Note that since we compute the eigenpairs for a child, */
912 /* all eigenvalue approximations are w.r.t the same shift. */
913 /* In this case, the entries in WORK should be used for */
914 /* computing the gaps since they exhibit even very small */
915 /* differences in the eigenvalues, as opposed to the */
916 /* entries in W which might "look" the same. */
917 if (k == 1)
918 {
919 /* In the case RANGE='I' and with not much initial */
920 /* accuracy in LAMBDA and VL, the formula */
921 /* LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA ) */
922 /* can lead to an overestimation of the left gap and */
923 /* thus to inadequately early RQI 'convergence'. */
924 /* Prevent this by forcing a small left gap. */
925 /* Computing MAX */
926 r__1 = f2c_abs(left);
927 r__2 = f2c_abs(right); // , expr subst
928 lgap = eps * max(r__1,r__2);
929 }
930 else
931 {
932 lgap = wgap[windmn];
933 }
934 if (k == im)
935 {
936 /* In the case RANGE='I' and with not much initial */
937 /* accuracy in LAMBDA and VU, the formula */
938 /* can lead to an overestimation of the right gap and */
939 /* thus to inadequately early RQI 'convergence'. */
940 /* Prevent this by forcing a small right gap. */
941 /* Computing MAX */
942 r__1 = f2c_abs(left);
943 r__2 = f2c_abs(right); // , expr subst
944 rgap = eps * max(r__1,r__2);
945 }
946 else
947 {
948 rgap = wgap[windex];
949 }
950 gap = min(lgap,rgap);
951 if (k == 1 || k == im)
952 {
953 /* The eigenvector support can become wrong */
954 /* because significant entries could be cut off due to a */
955 /* large GAPTOL parameter in LAR1V. Prevent this. */
956 gaptol = 0.f;
957 }
958 else
959 {
960 gaptol = gap * eps;
961 }
962 isupmn = in;
963 isupmx = 1;
964 /* Update WGAP so that it holds the minimum gap */
965 /* to the left or the right. This is crucial in the */
966 /* case where bisection is used to ensure that the */
967 /* eigenvalue is refined up to the required precision. */
968 /* The correct value is restored afterwards. */
969 savgap = wgap[windex];
970 wgap[windex] = gap;
971 /* We want to use the Rayleigh Quotient Correction */
972 /* as often as possible since it converges quadratically */
973 /* when we are close enough to the desired eigenvalue. */
974 /* However, the Rayleigh Quotient can have the wrong sign */
975 /* and lead us away from the desired eigenvalue. In this */
976 /* case, the best we can do is to use bisection. */
977 usedbs = FALSE_;
978 usedrq = FALSE_;
979 /* Bisection is initially turned off unless it is forced */
980 needbs = ! tryrqc;
981 L120: /* Check if bisection should be used to refine eigenvalue */
982 if (needbs)
983 {
984 /* Take the bisection as new iterate */
985 usedbs = TRUE_;
986 itmp1 = iwork[iindr + windex];
987 offset = indexw[wbegin] - 1;
988 r__1 = eps * 2.f;
989 slarrb_(&in, &d__[ibegin], &work[indlld + ibegin - 1], &indeig, &indeig, &c_b28, &r__1, & offset, &work[wbegin], &wgap[wbegin], & werr[wbegin], &work[indwrk], &iwork[ iindwk], pivmin, &spdiam, &itmp1, &iinfo);
990 if (iinfo != 0)
991 {
992 *info = -3;
993 return 0;
994 }
995 lambda = work[windex];
996 /* Reset twist index from inaccurate LAMBDA to */
997 /* force computation of true MINGMA */
998 iwork[iindr + windex] = 0;
999 }
1000 /* Given LAMBDA, compute the eigenvector. */
1001 L__1 = ! usedbs;
1002 clar1v_(&in, &c__1, &in, &lambda, &d__[ibegin], &l[ ibegin], &work[indld + ibegin - 1], &work[ indlld + ibegin - 1], pivmin, &gaptol, &z__[ ibegin + windex * z_dim1], &L__1, &negcnt, & ztz, &mingma, &iwork[iindr + windex], &isuppz[ (windex << 1) - 1], &nrminv, &resid, &rqcorr, &work[indwrk]);
1003 if (iter == 0)
1004 {
1005 bstres = resid;
1006 bstw = lambda;
1007 }
1008 else if (resid < bstres)
1009 {
1010 bstres = resid;
1011 bstw = lambda;
1012 }
1013 /* Computing MIN */
1014 i__4 = isupmn;
1015 i__5 = isuppz[(windex << 1) - 1]; // , expr subst
1016 isupmn = min(i__4,i__5);
1017 /* Computing MAX */
1018 i__4 = isupmx;
1019 i__5 = isuppz[windex * 2]; // , expr subst
1020 isupmx = max(i__4,i__5);
1021 ++iter;
1022 /* sin alpha <= |resid|/gap */
1023 /* Note that both the residual and the gap are */
1024 /* proportional to the matrix, so ||T|| doesn't play */
1025 /* a role in the quotient */
1026 /* Convergence test for Rayleigh-Quotient iteration */
1027 /* (omitted when Bisection has been used) */
1028 if (resid > tol * gap && f2c_abs(rqcorr) > rqtol * f2c_abs( lambda) && ! usedbs)
1029 {
1030 /* We need to check that the RQCORR update doesn't */
1031 /* move the eigenvalue away from the desired one and */
1032 /* towards a neighbor. -> protection with bisection */
1033 if (indeig <= negcnt)
1034 {
1035 /* The wanted eigenvalue lies to the left */
1036 sgndef = -1.f;
1037 }
1038 else
1039 {
1040 /* The wanted eigenvalue lies to the right */
1041 sgndef = 1.f;
1042 }
1043 /* We only use the RQCORR if it improves the */
1044 /* the iterate reasonably. */
1045 if (rqcorr * sgndef >= 0.f && lambda + rqcorr <= right && lambda + rqcorr >= left)
1046 {
1047 usedrq = TRUE_;
1048 /* Store new midpoint of bisection interval in WORK */
1049 if (sgndef == 1.f)
1050 {
1051 /* The current LAMBDA is on the left of the true */
1052 /* eigenvalue */
1053 left = lambda;
1054 /* We prefer to assume that the error estimate */
1055 /* is correct. We could make the interval not */
1056 /* as a bracket but to be modified if the RQCORR */
1057 /* chooses to. In this case, the RIGHT side should */
1058 /* be modified as follows: */
1059 /* RIGHT = MAX(RIGHT, LAMBDA + RQCORR) */
1060 }
1061 else
1062 {
1063 /* The current LAMBDA is on the right of the true */
1064 /* eigenvalue */
1065 right = lambda;
1066 /* See comment about assuming the error estimate is */
1067 /* correct above. */
1068 /* LEFT = MIN(LEFT, LAMBDA + RQCORR) */
1069 }
1070 work[windex] = (right + left) * .5f;
1071 /* Take RQCORR since it has the correct sign and */
1072 /* improves the iterate reasonably */
1073 lambda += rqcorr;
1074 /* Update width of error interval */
1075 werr[windex] = (right - left) * .5f;
1076 }
1077 else
1078 {
1079 needbs = TRUE_;
1080 }
1081 if (right - left < rqtol * f2c_abs(lambda))
1082 {
1083 /* The eigenvalue is computed to bisection accuracy */
1084 /* compute eigenvector and stop */
1085 usedbs = TRUE_;
1086 goto L120;
1087 }
1088 else if (iter < 10)
1089 {
1090 goto L120;
1091 }
1092 else if (iter == 10)
1093 {
1094 needbs = TRUE_;
1095 goto L120;
1096 }
1097 else
1098 {
1099 *info = 5;
1100 return 0;
1101 }
1102 }
1103 else
1104 {
1105 stp2ii = FALSE_;
1106 if (usedrq && usedbs && bstres <= resid)
1107 {
1108 lambda = bstw;
1109 stp2ii = TRUE_;
1110 }
1111 if (stp2ii)
1112 {
1113 /* improve error angle by second step */
1114 L__1 = ! usedbs;
1115 clar1v_(&in, &c__1, &in, &lambda, &d__[ibegin] , &l[ibegin], &work[indld + ibegin - 1], &work[indlld + ibegin - 1], pivmin, &gaptol, &z__[ibegin + windex * z_dim1], &L__1, &negcnt, &ztz, & mingma, &iwork[iindr + windex], & isuppz[(windex << 1) - 1], &nrminv, & resid, &rqcorr, &work[indwrk]);
1116 }
1117 work[windex] = lambda;
1118 }
1119 /* Compute FP-vector support w.r.t. whole matrix */
1120 isuppz[(windex << 1) - 1] += oldien;
1121 isuppz[windex * 2] += oldien;
1122 zfrom = isuppz[(windex << 1) - 1];
1123 zto = isuppz[windex * 2];
1124 isupmn += oldien;
1125 isupmx += oldien;
1126 /* Ensure vector is ok if support in the RQI has changed */
1127 if (isupmn < zfrom)
1128 {
1129 i__4 = zfrom - 1;
1130 for (ii = isupmn;
1131 ii <= i__4;
1132 ++ii)
1133 {
1134 i__5 = ii + windex * z_dim1;
1135 z__[i__5].r = 0.f;
1136 z__[i__5].i = 0.f; // , expr subst
1137 /* L122: */
1138 }
1139 }
1140 if (isupmx > zto)
1141 {
1142 i__4 = isupmx;
1143 for (ii = zto + 1;
1144 ii <= i__4;
1145 ++ii)
1146 {
1147 i__5 = ii + windex * z_dim1;
1148 z__[i__5].r = 0.f;
1149 z__[i__5].i = 0.f; // , expr subst
1150 /* L123: */
1151 }
1152 }
1153 i__4 = zto - zfrom + 1;
1154 csscal_(&i__4, &nrminv, &z__[zfrom + windex * z_dim1], &c__1);
1155 L125: /* Update W */
1156 w[windex] = lambda + sigma;
1157 /* Recompute the gaps on the left and right */
1158 /* But only allow them to become larger and not */
1159 /* smaller (which can only happen through "bad" */
1160 /* cancellation and doesn't reflect the theory */
1161 /* where the initial gaps are underestimated due */
1162 /* to WERR being too crude.) */
1163 if (! eskip)
1164 {
1165 if (k > 1)
1166 {
1167 /* Computing MAX */
1168 r__1 = wgap[windmn];
1169 r__2 = w[windex] - werr[ windex] - w[windmn] - werr[windmn]; // , expr subst
1170 wgap[windmn] = max(r__1,r__2);
1171 }
1172 if (windex < wend)
1173 {
1174 /* Computing MAX */
1175 r__1 = savgap;
1176 r__2 = w[windpl] - werr[windpl] - w[windex] - werr[windex]; // , expr subst
1177 wgap[windex] = max(r__1,r__2);
1178 }
1179 }
1180 ++idone;
1181 }
1182 /* here ends the code for the current child */
1183 L139: /* Proceed to any remaining child nodes */
1184 newfst = j + 1;
1185 L140:
1186 ;
1187 }
1188 /* L150: */
1189 }
1190 ++ndepth;
1191 goto L40;
1192 }
1193 ibegin = iend + 1;
1194 wbegin = wend + 1;
1195 L170:
1196 ;
1197 }
1198 return 0;
1199 /* End of CLARRV */
1200 }
1201 /* clarrv_ */
1202