1 /* ./src_f77/sstein.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 integer c__2 = 2;
11 static integer c__1 = 1;
12 static integer c_n1 = -1;
13
sstein_(integer * n,real * d__,real * e,integer * m,real * w,integer * iblock,integer * isplit,real * z__,integer * ldz,real * work,integer * iwork,integer * ifail,integer * info)14 /* Subroutine */ int sstein_(integer *n, real *d__, real *e, integer *m, real
15 *w, integer *iblock, integer *isplit, real *z__, integer *ldz, real *
16 work, integer *iwork, integer *ifail, integer *info)
17 {
18 /* System generated locals */
19 integer z_dim1, z_offset, i__1, i__2, i__3;
20 real r__1, r__2, r__3, r__4, r__5;
21
22 /* Builtin functions */
23 double sqrt(doublereal);
24
25 /* Local variables */
26 static integer i__, j, b1, j1, bn;
27 static real xj, scl, eps, ctr, sep, nrm, tol;
28 static integer its;
29 static real xjm, eps1;
30 static integer jblk, nblk, jmax;
31 extern doublereal sdot_(integer *, real *, integer *, real *, integer *),
32 snrm2_(integer *, real *, integer *);
33 static integer iseed[4], gpind, iinfo;
34 extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
35 extern doublereal sasum_(integer *, real *, integer *);
36 extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
37 integer *);
38 static real ortol;
39 extern /* Subroutine */ int saxpy_(integer *, real *, real *, integer *,
40 real *, integer *);
41 static integer indrv1, indrv2, indrv3, indrv4, indrv5;
42 extern doublereal slamch_(char *, ftnlen);
43 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen), slagtf_(
44 integer *, real *, real *, real *, real *, real *, real *,
45 integer *, integer *);
46 static integer nrmchk;
47 extern integer isamax_(integer *, real *, integer *);
48 extern /* Subroutine */ int slagts_(integer *, integer *, real *, real *,
49 real *, real *, integer *, real *, real *, integer *);
50 static integer blksiz;
51 static real onenrm, pertol;
52 extern /* Subroutine */ int slarnv_(integer *, integer *, integer *, real
53 *);
54 static real stpcrt;
55
56
57 /* -- LAPACK routine (version 3.0) -- */
58 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
59 /* Courant Institute, Argonne National Lab, and Rice University */
60 /* September 30, 1994 */
61
62 /* .. Scalar Arguments .. */
63 /* .. */
64 /* .. Array Arguments .. */
65 /* .. */
66
67 /* Purpose */
68 /* ======= */
69
70 /* SSTEIN computes the eigenvectors of a real symmetric tridiagonal */
71 /* matrix T corresponding to specified eigenvalues, using inverse */
72 /* iteration. */
73
74 /* The maximum number of iterations allowed for each eigenvector is */
75 /* specified by an internal parameter MAXITS (currently set to 5). */
76
77 /* Arguments */
78 /* ========= */
79
80 /* N (input) INTEGER */
81 /* The order of the matrix. N >= 0. */
82
83 /* D (input) REAL array, dimension (N) */
84 /* The n diagonal elements of the tridiagonal matrix T. */
85
86 /* E (input) REAL array, dimension (N) */
87 /* The (n-1) subdiagonal elements of the tridiagonal matrix */
88 /* T, in elements 1 to N-1. E(N) need not be set. */
89
90 /* M (input) INTEGER */
91 /* The number of eigenvectors to be found. 0 <= M <= N. */
92
93 /* W (input) REAL array, dimension (N) */
94 /* The first M elements of W contain the eigenvalues for */
95 /* which eigenvectors are to be computed. The eigenvalues */
96 /* should be grouped by split-off block and ordered from */
97 /* smallest to largest within the block. ( The output array */
98 /* W from SSTEBZ with ORDER = 'B' is expected here. ) */
99
100 /* IBLOCK (input) INTEGER array, dimension (N) */
101 /* The submatrix indices associated with the corresponding */
102 /* eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to */
103 /* the first submatrix from the top, =2 if W(i) belongs to */
104 /* the second submatrix, etc. ( The output array IBLOCK */
105 /* from SSTEBZ is expected here. ) */
106
107 /* ISPLIT (input) INTEGER array, dimension (N) */
108 /* The splitting points, at which T breaks up into submatrices. */
109 /* The first submatrix consists of rows/columns 1 to */
110 /* ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 */
111 /* through ISPLIT( 2 ), etc. */
112 /* ( The output array ISPLIT from SSTEBZ is expected here. ) */
113
114 /* Z (output) REAL array, dimension (LDZ, M) */
115 /* The computed eigenvectors. The eigenvector associated */
116 /* with the eigenvalue W(i) is stored in the i-th column of */
117 /* Z. Any vector which fails to converge is set to its current */
118 /* iterate after MAXITS iterations. */
119
120 /* LDZ (input) INTEGER */
121 /* The leading dimension of the array Z. LDZ >= max(1,N). */
122
123 /* WORK (workspace) REAL array, dimension (5*N) */
124
125 /* IWORK (workspace) INTEGER array, dimension (N) */
126
127 /* IFAIL (output) INTEGER array, dimension (M) */
128 /* On normal exit, all elements of IFAIL are zero. */
129 /* If one or more eigenvectors fail to converge after */
130 /* MAXITS iterations, then their indices are stored in */
131 /* array IFAIL. */
132
133 /* INFO (output) INTEGER */
134 /* = 0: successful exit. */
135 /* < 0: if INFO = -i, the i-th argument had an illegal value */
136 /* > 0: if INFO = i, then i eigenvectors failed to converge */
137 /* in MAXITS iterations. Their indices are stored in */
138 /* array IFAIL. */
139
140 /* Internal Parameters */
141 /* =================== */
142
143 /* MAXITS INTEGER, default = 5 */
144 /* The maximum number of iterations performed. */
145
146 /* EXTRA INTEGER, default = 2 */
147 /* The number of iterations performed after norm growth */
148 /* criterion is satisfied, should be at least 1. */
149
150 /* ===================================================================== */
151
152 /* .. Parameters .. */
153 /* .. */
154 /* .. Local Scalars .. */
155 /* .. */
156 /* .. Local Arrays .. */
157 /* .. */
158 /* .. External Functions .. */
159 /* .. */
160 /* .. External Subroutines .. */
161 /* .. */
162 /* .. Intrinsic Functions .. */
163 /* .. */
164 /* .. Executable Statements .. */
165
166 /* Test the input parameters. */
167
168 /* Parameter adjustments */
169 --d__;
170 --e;
171 --w;
172 --iblock;
173 --isplit;
174 z_dim1 = *ldz;
175 z_offset = 1 + z_dim1;
176 z__ -= z_offset;
177 --work;
178 --iwork;
179 --ifail;
180
181 /* Function Body */
182 *info = 0;
183 i__1 = *m;
184 for (i__ = 1; i__ <= i__1; ++i__) {
185 ifail[i__] = 0;
186 /* L10: */
187 }
188
189 if (*n < 0) {
190 *info = -1;
191 } else if (*m < 0 || *m > *n) {
192 *info = -4;
193 } else if (*ldz < max(1,*n)) {
194 *info = -9;
195 } else {
196 i__1 = *m;
197 for (j = 2; j <= i__1; ++j) {
198 if (iblock[j] < iblock[j - 1]) {
199 *info = -6;
200 goto L30;
201 }
202 if (iblock[j] == iblock[j - 1] && w[j] < w[j - 1]) {
203 *info = -5;
204 goto L30;
205 }
206 /* L20: */
207 }
208 L30:
209 ;
210 }
211
212 if (*info != 0) {
213 i__1 = -(*info);
214 xerbla_("SSTEIN", &i__1, (ftnlen)6);
215 return 0;
216 }
217
218 /* Quick return if possible */
219
220 if (*n == 0 || *m == 0) {
221 return 0;
222 } else if (*n == 1) {
223 z__[z_dim1 + 1] = 1.f;
224 return 0;
225 }
226
227 /* Get machine constants. */
228
229 eps = slamch_("Precision", (ftnlen)9);
230
231 /* Initialize seed for random number generator SLARNV. */
232
233 for (i__ = 1; i__ <= 4; ++i__) {
234 iseed[i__ - 1] = 1;
235 /* L40: */
236 }
237
238 /* Initialize pointers. */
239
240 indrv1 = 0;
241 indrv2 = indrv1 + *n;
242 indrv3 = indrv2 + *n;
243 indrv4 = indrv3 + *n;
244 indrv5 = indrv4 + *n;
245
246 /* Compute eigenvectors of matrix blocks. */
247
248 j1 = 1;
249 i__1 = iblock[*m];
250 for (nblk = 1; nblk <= i__1; ++nblk) {
251
252 /* Find starting and ending indices of block nblk. */
253
254 if (nblk == 1) {
255 b1 = 1;
256 } else {
257 b1 = isplit[nblk - 1] + 1;
258 }
259 bn = isplit[nblk];
260 blksiz = bn - b1 + 1;
261 if (blksiz == 1) {
262 goto L60;
263 }
264 gpind = b1;
265
266 /* Compute reorthogonalization criterion and stopping criterion. */
267
268 onenrm = (r__1 = d__[b1], dabs(r__1)) + (r__2 = e[b1], dabs(r__2));
269 /* Computing MAX */
270 r__3 = onenrm, r__4 = (r__1 = d__[bn], dabs(r__1)) + (r__2 = e[bn - 1]
271 , dabs(r__2));
272 onenrm = dmax(r__3,r__4);
273 i__2 = bn - 1;
274 for (i__ = b1 + 1; i__ <= i__2; ++i__) {
275 /* Computing MAX */
276 r__4 = onenrm, r__5 = (r__1 = d__[i__], dabs(r__1)) + (r__2 = e[
277 i__ - 1], dabs(r__2)) + (r__3 = e[i__], dabs(r__3));
278 onenrm = dmax(r__4,r__5);
279 /* L50: */
280 }
281 ortol = onenrm * .001f;
282
283 stpcrt = sqrt(.1f / blksiz);
284
285 /* Loop through eigenvalues of block nblk. */
286
287 L60:
288 jblk = 0;
289 i__2 = *m;
290 for (j = j1; j <= i__2; ++j) {
291 if (iblock[j] != nblk) {
292 j1 = j;
293 goto L160;
294 }
295 ++jblk;
296 xj = w[j];
297
298 /* Skip all the work if the block size is one. */
299
300 if (blksiz == 1) {
301 work[indrv1 + 1] = 1.f;
302 goto L120;
303 }
304
305 /* If eigenvalues j and j-1 are too close, add a relatively */
306 /* small perturbation. */
307
308 if (jblk > 1) {
309 eps1 = (r__1 = eps * xj, dabs(r__1));
310 pertol = eps1 * 10.f;
311 sep = xj - xjm;
312 if (sep < pertol) {
313 xj = xjm + pertol;
314 }
315 }
316
317 its = 0;
318 nrmchk = 0;
319
320 /* Get random starting vector. */
321
322 slarnv_(&c__2, iseed, &blksiz, &work[indrv1 + 1]);
323
324 /* Copy the matrix T so it won't be destroyed in factorization. */
325
326 scopy_(&blksiz, &d__[b1], &c__1, &work[indrv4 + 1], &c__1);
327 i__3 = blksiz - 1;
328 scopy_(&i__3, &e[b1], &c__1, &work[indrv2 + 2], &c__1);
329 i__3 = blksiz - 1;
330 scopy_(&i__3, &e[b1], &c__1, &work[indrv3 + 1], &c__1);
331
332 /* Compute LU factors with partial pivoting ( PT = LU ) */
333
334 tol = 0.f;
335 slagtf_(&blksiz, &work[indrv4 + 1], &xj, &work[indrv2 + 2], &work[
336 indrv3 + 1], &tol, &work[indrv5 + 1], &iwork[1], &iinfo);
337
338 /* Update iteration count. */
339
340 L70:
341 ++its;
342 if (its > 5) {
343 goto L100;
344 }
345
346 /* Normalize and scale the righthand side vector Pb. */
347
348 /* Computing MAX */
349 r__2 = eps, r__3 = (r__1 = work[indrv4 + blksiz], dabs(r__1));
350 scl = blksiz * onenrm * dmax(r__2,r__3) / sasum_(&blksiz, &work[
351 indrv1 + 1], &c__1);
352 sscal_(&blksiz, &scl, &work[indrv1 + 1], &c__1);
353
354 /* Solve the system LU = Pb. */
355
356 slagts_(&c_n1, &blksiz, &work[indrv4 + 1], &work[indrv2 + 2], &
357 work[indrv3 + 1], &work[indrv5 + 1], &iwork[1], &work[
358 indrv1 + 1], &tol, &iinfo);
359
360 /* Reorthogonalize by modified Gram-Schmidt if eigenvalues are */
361 /* close enough. */
362
363 if (jblk == 1) {
364 goto L90;
365 }
366 if ((r__1 = xj - xjm, dabs(r__1)) > ortol) {
367 gpind = j;
368 }
369 if (gpind != j) {
370 i__3 = j - 1;
371 for (i__ = gpind; i__ <= i__3; ++i__) {
372 ctr = -sdot_(&blksiz, &work[indrv1 + 1], &c__1, &z__[b1 +
373 i__ * z_dim1], &c__1);
374 saxpy_(&blksiz, &ctr, &z__[b1 + i__ * z_dim1], &c__1, &
375 work[indrv1 + 1], &c__1);
376 /* L80: */
377 }
378 }
379
380 /* Check the infinity norm of the iterate. */
381
382 L90:
383 jmax = isamax_(&blksiz, &work[indrv1 + 1], &c__1);
384 nrm = (r__1 = work[indrv1 + jmax], dabs(r__1));
385
386 /* Continue for additional iterations after norm reaches */
387 /* stopping criterion. */
388
389 if (nrm < stpcrt) {
390 goto L70;
391 }
392 ++nrmchk;
393 if (nrmchk < 3) {
394 goto L70;
395 }
396
397 goto L110;
398
399 /* If stopping criterion was not satisfied, update info and */
400 /* store eigenvector number in array ifail. */
401
402 L100:
403 ++(*info);
404 ifail[*info] = j;
405
406 /* Accept iterate as jth eigenvector. */
407
408 L110:
409 scl = 1.f / snrm2_(&blksiz, &work[indrv1 + 1], &c__1);
410 jmax = isamax_(&blksiz, &work[indrv1 + 1], &c__1);
411 if (work[indrv1 + jmax] < 0.f) {
412 scl = -scl;
413 }
414 sscal_(&blksiz, &scl, &work[indrv1 + 1], &c__1);
415 L120:
416 i__3 = *n;
417 for (i__ = 1; i__ <= i__3; ++i__) {
418 z__[i__ + j * z_dim1] = 0.f;
419 /* L130: */
420 }
421 i__3 = blksiz;
422 for (i__ = 1; i__ <= i__3; ++i__) {
423 z__[b1 + i__ - 1 + j * z_dim1] = work[indrv1 + i__];
424 /* L140: */
425 }
426
427 /* Save the shift to check eigenvalue spacing at next */
428 /* iteration. */
429
430 xjm = xj;
431
432 /* L150: */
433 }
434 L160:
435 ;
436 }
437
438 return 0;
439
440 /* End of SSTEIN */
441
442 } /* sstein_ */
443
444