1 /* ../netlib/zstein.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 integer c__2 = 2;
5 static integer c__1 = 1;
6 static integer c_n1 = -1;
7 /* > \brief \b ZSTEIN */
8 /* =========== DOCUMENTATION =========== */
9 /* Online html documentation available at */
10 /* http://www.netlib.org/lapack/explore-html/ */
11 /* > \htmlonly */
12 /* > Download ZSTEIN + dependencies */
13 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstein. f"> */
14 /* > [TGZ]</a> */
15 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstein. f"> */
16 /* > [ZIP]</a> */
17 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstein. f"> */
18 /* > [TXT]</a> */
19 /* > \endhtmlonly */
20 /* Definition: */
21 /* =========== */
22 /* SUBROUTINE ZSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, */
23 /* IWORK, IFAIL, INFO ) */
24 /* .. Scalar Arguments .. */
25 /* INTEGER INFO, LDZ, M, N */
26 /* .. */
27 /* .. Array Arguments .. */
28 /* INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ), */
29 /* $ IWORK( * ) */
30 /* DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) */
31 /* COMPLEX*16 Z( LDZ, * ) */
32 /* .. */
33 /* > \par Purpose: */
34 /* ============= */
35 /* > */
36 /* > \verbatim */
37 /* > */
38 /* > ZSTEIN computes the eigenvectors of a real symmetric tridiagonal */
39 /* > matrix T corresponding to specified eigenvalues, using inverse */
40 /* > iteration. */
41 /* > */
42 /* > The maximum number of iterations allowed for each eigenvector is */
43 /* > specified by an internal parameter MAXITS (currently set to 5). */
44 /* > */
45 /* > Although the eigenvectors are real, they are stored in a complex */
46 /* > array, which may be passed to ZUNMTR or ZUPMTR for back */
47 /* > transformation to the eigenvectors of a complex Hermitian matrix */
48 /* > which was reduced to tridiagonal form. */
49 /* > */
50 /* > \endverbatim */
51 /* Arguments: */
52 /* ========== */
53 /* > \param[in] N */
54 /* > \verbatim */
55 /* > N is INTEGER */
56 /* > The order of the matrix. N >= 0. */
57 /* > \endverbatim */
58 /* > */
59 /* > \param[in] D */
60 /* > \verbatim */
61 /* > D is DOUBLE PRECISION array, dimension (N) */
62 /* > The n diagonal elements of the tridiagonal matrix T. */
63 /* > \endverbatim */
64 /* > */
65 /* > \param[in] E */
66 /* > \verbatim */
67 /* > E is DOUBLE PRECISION array, dimension (N-1) */
68 /* > The (n-1) subdiagonal elements of the tridiagonal matrix */
69 /* > T, stored in elements 1 to N-1. */
70 /* > \endverbatim */
71 /* > */
72 /* > \param[in] M */
73 /* > \verbatim */
74 /* > M is INTEGER */
75 /* > The number of eigenvectors to be found. 0 <= M <= N. */
76 /* > \endverbatim */
77 /* > */
78 /* > \param[in] W */
79 /* > \verbatim */
80 /* > W is DOUBLE PRECISION array, dimension (N) */
81 /* > The first M elements of W contain the eigenvalues for */
82 /* > which eigenvectors are to be computed. The eigenvalues */
83 /* > should be grouped by split-off block and ordered from */
84 /* > smallest to largest within the block. ( The output array */
85 /* > W from DSTEBZ with ORDER = 'B' is expected here. ) */
86 /* > \endverbatim */
87 /* > */
88 /* > \param[in] IBLOCK */
89 /* > \verbatim */
90 /* > IBLOCK is INTEGER array, dimension (N) */
91 /* > The submatrix indices associated with the corresponding */
92 /* > eigenvalues in W;
93 IBLOCK(i)=1 if eigenvalue W(i) belongs to */
94 /* > the first submatrix from the top, =2 if W(i) belongs to */
95 /* > the second submatrix, etc. ( The output array IBLOCK */
96 /* > from DSTEBZ is expected here. ) */
97 /* > \endverbatim */
98 /* > */
99 /* > \param[in] ISPLIT */
100 /* > \verbatim */
101 /* > ISPLIT is INTEGER array, dimension (N) */
102 /* > The splitting points, at which T breaks up into submatrices. */
103 /* > The first submatrix consists of rows/columns 1 to */
104 /* > ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 */
105 /* > through ISPLIT( 2 ), etc. */
106 /* > ( The output array ISPLIT from DSTEBZ is expected here. ) */
107 /* > \endverbatim */
108 /* > */
109 /* > \param[out] Z */
110 /* > \verbatim */
111 /* > Z is COMPLEX*16 array, dimension (LDZ, M) */
112 /* > The computed eigenvectors. The eigenvector associated */
113 /* > with the eigenvalue W(i) is stored in the i-th column of */
114 /* > Z. Any vector which fails to converge is set to its current */
115 /* > iterate after MAXITS iterations. */
116 /* > The imaginary parts of the eigenvectors are set to zero. */
117 /* > \endverbatim */
118 /* > */
119 /* > \param[in] LDZ */
120 /* > \verbatim */
121 /* > LDZ is INTEGER */
122 /* > The leading dimension of the array Z. LDZ >= max(1,N). */
123 /* > \endverbatim */
124 /* > */
125 /* > \param[out] WORK */
126 /* > \verbatim */
127 /* > WORK is DOUBLE PRECISION array, dimension (5*N) */
128 /* > \endverbatim */
129 /* > */
130 /* > \param[out] IWORK */
131 /* > \verbatim */
132 /* > IWORK is INTEGER array, dimension (N) */
133 /* > \endverbatim */
134 /* > */
135 /* > \param[out] IFAIL */
136 /* > \verbatim */
137 /* > IFAIL is INTEGER array, dimension (M) */
138 /* > On normal exit, all elements of IFAIL are zero. */
139 /* > If one or more eigenvectors fail to converge after */
140 /* > MAXITS iterations, then their indices are stored in */
141 /* > array IFAIL. */
142 /* > \endverbatim */
143 /* > */
144 /* > \param[out] INFO */
145 /* > \verbatim */
146 /* > INFO is INTEGER */
147 /* > = 0: successful exit */
148 /* > < 0: if INFO = -i, the i-th argument had an illegal value */
149 /* > > 0: if INFO = i, then i eigenvectors failed to converge */
150 /* > in MAXITS iterations. Their indices are stored in */
151 /* > array IFAIL. */
152 /* > \endverbatim */
153 /* > \par Internal Parameters: */
154 /* ========================= */
155 /* > */
156 /* > \verbatim */
157 /* > MAXITS INTEGER, default = 5 */
158 /* > The maximum number of iterations performed. */
159 /* > */
160 /* > EXTRA INTEGER, default = 2 */
161 /* > The number of iterations performed after norm growth */
162 /* > criterion is satisfied, should be at least 1. */
163 /* > \endverbatim */
164 /* Authors: */
165 /* ======== */
166 /* > \author Univ. of Tennessee */
167 /* > \author Univ. of California Berkeley */
168 /* > \author Univ. of Colorado Denver */
169 /* > \author NAG Ltd. */
170 /* > \date November 2011 */
171 /* > \ingroup complex16OTHERcomputational */
172 /* ===================================================================== */
173 /* Subroutine */
zstein_(integer * n,doublereal * d__,doublereal * e,integer * m,doublereal * w,integer * iblock,integer * isplit,doublecomplex * z__,integer * ldz,doublereal * work,integer * iwork,integer * ifail,integer * info)174 int zstein_(integer *n, doublereal *d__, doublereal *e, integer *m, doublereal *w, integer *iblock, integer *isplit, doublecomplex *z__, integer *ldz, doublereal *work, integer *iwork, integer *ifail, integer *info)
175 {
176 /* System generated locals */
177 integer z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5;
178 doublereal d__1, d__2, d__3, d__4, d__5;
179 doublecomplex z__1;
180 /* Builtin functions */
181 double sqrt(doublereal);
182 /* Local variables */
183 integer i__, j, b1, j1, bn, jr;
184 doublereal xj, scl, eps, sep, nrm, tol;
185 integer its;
186 doublereal xjm, ztr, eps1;
187 integer jblk, nblk, jmax;
188 extern doublereal dnrm2_(integer *, doublereal *, integer *);
189 extern /* Subroutine */
190 int dscal_(integer *, doublereal *, doublereal *, integer *);
191 integer iseed[4], gpind, iinfo;
192 extern doublereal dasum_(integer *, doublereal *, integer *);
193 extern /* Subroutine */
194 int dcopy_(integer *, doublereal *, integer *, doublereal *, integer *);
195 doublereal ortol;
196 integer indrv1, indrv2, indrv3, indrv4, indrv5;
197 extern doublereal dlamch_(char *);
198 extern /* Subroutine */
199 int dlagtf_(integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, integer * , integer *);
200 extern integer idamax_(integer *, doublereal *, integer *);
201 extern /* Subroutine */
202 int xerbla_(char *, integer *), dlagts_( integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, integer *, doublereal *, doublereal *, integer *);
203 integer nrmchk;
204 extern /* Subroutine */
205 int dlarnv_(integer *, integer *, integer *, doublereal *);
206 integer blksiz;
207 doublereal onenrm, dtpcrt, pertol;
208 /* -- LAPACK computational routine (version 3.4.0) -- */
209 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
210 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
211 /* November 2011 */
212 /* .. Scalar Arguments .. */
213 /* .. */
214 /* .. Array Arguments .. */
215 /* .. */
216 /* ===================================================================== */
217 /* .. Parameters .. */
218 /* .. */
219 /* .. Local Scalars .. */
220 /* .. */
221 /* .. Local Arrays .. */
222 /* .. */
223 /* .. External Functions .. */
224 /* .. */
225 /* .. External Subroutines .. */
226 /* .. */
227 /* .. Intrinsic Functions .. */
228 /* .. */
229 /* .. Executable Statements .. */
230 /* Test the input parameters. */
231 /* Parameter adjustments */
232 --d__;
233 --e;
234 --w;
235 --iblock;
236 --isplit;
237 z_dim1 = *ldz;
238 z_offset = 1 + z_dim1;
239 z__ -= z_offset;
240 --work;
241 --iwork;
242 --ifail;
243 /* Function Body */
244 *info = 0;
245 i__1 = *m;
246 for (i__ = 1;
247 i__ <= i__1;
248 ++i__)
249 {
250 ifail[i__] = 0;
251 /* L10: */
252 }
253 if (*n < 0)
254 {
255 *info = -1;
256 }
257 else if (*m < 0 || *m > *n)
258 {
259 *info = -4;
260 }
261 else if (*ldz < max(1,*n))
262 {
263 *info = -9;
264 }
265 else
266 {
267 i__1 = *m;
268 for (j = 2;
269 j <= i__1;
270 ++j)
271 {
272 if (iblock[j] < iblock[j - 1])
273 {
274 *info = -6;
275 goto L30;
276 }
277 if (iblock[j] == iblock[j - 1] && w[j] < w[j - 1])
278 {
279 *info = -5;
280 goto L30;
281 }
282 /* L20: */
283 }
284 L30:
285 ;
286 }
287 if (*info != 0)
288 {
289 i__1 = -(*info);
290 xerbla_("ZSTEIN", &i__1);
291 return 0;
292 }
293 /* Quick return if possible */
294 if (*n == 0 || *m == 0)
295 {
296 return 0;
297 }
298 else if (*n == 1)
299 {
300 i__1 = z_dim1 + 1;
301 z__[i__1].r = 1.;
302 z__[i__1].i = 0.; // , expr subst
303 return 0;
304 }
305 /* Get machine constants. */
306 eps = dlamch_("Precision");
307 /* Initialize seed for random number generator DLARNV. */
308 for (i__ = 1;
309 i__ <= 4;
310 ++i__)
311 {
312 iseed[i__ - 1] = 1;
313 /* L40: */
314 }
315 /* Initialize pointers. */
316 indrv1 = 0;
317 indrv2 = indrv1 + *n;
318 indrv3 = indrv2 + *n;
319 indrv4 = indrv3 + *n;
320 indrv5 = indrv4 + *n;
321 /* Compute eigenvectors of matrix blocks. */
322 j1 = 1;
323 i__1 = iblock[*m];
324 for (nblk = 1;
325 nblk <= i__1;
326 ++nblk)
327 {
328 /* Find starting and ending indices of block nblk. */
329 if (nblk == 1)
330 {
331 b1 = 1;
332 }
333 else
334 {
335 b1 = isplit[nblk - 1] + 1;
336 }
337 bn = isplit[nblk];
338 blksiz = bn - b1 + 1;
339 if (blksiz == 1)
340 {
341 goto L60;
342 }
343 gpind = b1;
344 /* Compute reorthogonalization criterion and stopping criterion. */
345 onenrm = (d__1 = d__[b1], f2c_abs(d__1)) + (d__2 = e[b1], f2c_abs(d__2));
346 /* Computing MAX */
347 d__3 = onenrm;
348 d__4 = (d__1 = d__[bn], f2c_abs(d__1)) + (d__2 = e[bn - 1], f2c_abs(d__2)); // , expr subst
349 onenrm = max(d__3,d__4);
350 i__2 = bn - 1;
351 for (i__ = b1 + 1;
352 i__ <= i__2;
353 ++i__)
354 {
355 /* Computing MAX */
356 d__4 = onenrm;
357 d__5 = (d__1 = d__[i__], f2c_abs(d__1)) + (d__2 = e[ i__ - 1], f2c_abs(d__2)) + (d__3 = e[i__], f2c_abs(d__3)); // , expr subst
358 onenrm = max(d__4,d__5);
359 /* L50: */
360 }
361 ortol = onenrm * .001;
362 dtpcrt = sqrt(.1 / blksiz);
363 /* Loop through eigenvalues of block nblk. */
364 L60:
365 jblk = 0;
366 i__2 = *m;
367 for (j = j1;
368 j <= i__2;
369 ++j)
370 {
371 if (iblock[j] != nblk)
372 {
373 j1 = j;
374 goto L180;
375 }
376 ++jblk;
377 xj = w[j];
378 /* Skip all the work if the block size is one. */
379 if (blksiz == 1)
380 {
381 work[indrv1 + 1] = 1.;
382 goto L140;
383 }
384 /* If eigenvalues j and j-1 are too close, add a relatively */
385 /* small perturbation. */
386 if (jblk > 1)
387 {
388 eps1 = (d__1 = eps * xj, f2c_abs(d__1));
389 pertol = eps1 * 10.;
390 sep = xj - xjm;
391 if (sep < pertol)
392 {
393 xj = xjm + pertol;
394 }
395 }
396 its = 0;
397 nrmchk = 0;
398 /* Get random starting vector. */
399 dlarnv_(&c__2, iseed, &blksiz, &work[indrv1 + 1]);
400 /* Copy the matrix T so it won't be destroyed in factorization. */
401 dcopy_(&blksiz, &d__[b1], &c__1, &work[indrv4 + 1], &c__1);
402 i__3 = blksiz - 1;
403 dcopy_(&i__3, &e[b1], &c__1, &work[indrv2 + 2], &c__1);
404 i__3 = blksiz - 1;
405 dcopy_(&i__3, &e[b1], &c__1, &work[indrv3 + 1], &c__1);
406 /* Compute LU factors with partial pivoting ( PT = LU ) */
407 tol = 0.;
408 dlagtf_(&blksiz, &work[indrv4 + 1], &xj, &work[indrv2 + 2], &work[ indrv3 + 1], &tol, &work[indrv5 + 1], &iwork[1], &iinfo);
409 /* Update iteration count. */
410 L70:
411 ++its;
412 if (its > 5)
413 {
414 goto L120;
415 }
416 /* Normalize and scale the righthand side vector Pb. */
417 /* Computing MAX */
418 d__2 = eps;
419 d__3 = (d__1 = work[indrv4 + blksiz], f2c_abs(d__1)); // , expr subst
420 scl = blksiz * onenrm * max(d__2,d__3) / dasum_(&blksiz, &work[ indrv1 + 1], &c__1);
421 dscal_(&blksiz, &scl, &work[indrv1 + 1], &c__1);
422 /* Solve the system LU = Pb. */
423 dlagts_(&c_n1, &blksiz, &work[indrv4 + 1], &work[indrv2 + 2], & work[indrv3 + 1], &work[indrv5 + 1], &iwork[1], &work[ indrv1 + 1], &tol, &iinfo);
424 /* Reorthogonalize by modified Gram-Schmidt if eigenvalues are */
425 /* close enough. */
426 if (jblk == 1)
427 {
428 goto L110;
429 }
430 if ((d__1 = xj - xjm, f2c_abs(d__1)) > ortol)
431 {
432 gpind = j;
433 }
434 if (gpind != j)
435 {
436 i__3 = j - 1;
437 for (i__ = gpind;
438 i__ <= i__3;
439 ++i__)
440 {
441 ztr = 0.;
442 i__4 = blksiz;
443 for (jr = 1;
444 jr <= i__4;
445 ++jr)
446 {
447 i__5 = b1 - 1 + jr + i__ * z_dim1;
448 ztr += work[indrv1 + jr] * z__[i__5].r;
449 /* L80: */
450 }
451 i__4 = blksiz;
452 for (jr = 1;
453 jr <= i__4;
454 ++jr)
455 {
456 i__5 = b1 - 1 + jr + i__ * z_dim1;
457 work[indrv1 + jr] -= ztr * z__[i__5].r;
458 /* L90: */
459 }
460 /* L100: */
461 }
462 }
463 /* Check the infinity norm of the iterate. */
464 L110:
465 jmax = idamax_(&blksiz, &work[indrv1 + 1], &c__1);
466 nrm = (d__1 = work[indrv1 + jmax], f2c_abs(d__1));
467 /* Continue for additional iterations after norm reaches */
468 /* stopping criterion. */
469 if (nrm < dtpcrt)
470 {
471 goto L70;
472 }
473 ++nrmchk;
474 if (nrmchk < 3)
475 {
476 goto L70;
477 }
478 goto L130;
479 /* If stopping criterion was not satisfied, update info and */
480 /* store eigenvector number in array ifail. */
481 L120:
482 ++(*info);
483 ifail[*info] = j;
484 /* Accept iterate as jth eigenvector. */
485 L130:
486 scl = 1. / dnrm2_(&blksiz, &work[indrv1 + 1], &c__1);
487 jmax = idamax_(&blksiz, &work[indrv1 + 1], &c__1);
488 if (work[indrv1 + jmax] < 0.)
489 {
490 scl = -scl;
491 }
492 dscal_(&blksiz, &scl, &work[indrv1 + 1], &c__1);
493 L140:
494 i__3 = *n;
495 for (i__ = 1;
496 i__ <= i__3;
497 ++i__)
498 {
499 i__4 = i__ + j * z_dim1;
500 z__[i__4].r = 0.;
501 z__[i__4].i = 0.; // , expr subst
502 /* L150: */
503 }
504 i__3 = blksiz;
505 for (i__ = 1;
506 i__ <= i__3;
507 ++i__)
508 {
509 i__4 = b1 + i__ - 1 + j * z_dim1;
510 i__5 = indrv1 + i__;
511 z__1.r = work[i__5];
512 z__1.i = 0.; // , expr subst
513 z__[i__4].r = z__1.r;
514 z__[i__4].i = z__1.i; // , expr subst
515 /* L160: */
516 }
517 /* Save the shift to check eigenvalue spacing at next */
518 /* iteration. */
519 xjm = xj;
520 /* L170: */
521 }
522 L180:
523 ;
524 }
525 return 0;
526 /* End of ZSTEIN */
527 }
528 /* zstein_ */
529