3 /* Table of constant values */
5 static integer c__2 = 2;
6 static integer c__1 = 1;
7 static integer c_n1 = -1;
9 /* Subroutine */ int dstein_(integer *n, doublereal *d__, doublereal *e,
10 integer *m, doublereal *w, integer *iblock, integer *isplit,
11 doublereal *z__, integer *ldz, doublereal *work, integer *iwork,
12 integer *ifail, integer *info)
14 /* System generated locals */
15 integer z_dim1, z_offset, i__1, i__2, i__3;
16 doublereal d__1, d__2, d__3, d__4, d__5;
18 /* Builtin functions */
19 double sqrt(doublereal);
22 integer i__, j, b1, j1, bn;
23 doublereal xj, scl, eps, sep, nrm, tol;
25 doublereal xjm, ztr, eps1;
27 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
30 extern doublereal dnrm2_(integer *, doublereal *, integer *);
31 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
33 integer iseed[4], gpind, iinfo;
34 extern doublereal dasum_(integer *, doublereal *, integer *);
35 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
36 doublereal *, integer *), daxpy_(integer *, doublereal *,
37 doublereal *, integer *, doublereal *, integer *);
39 integer indrv1, indrv2, indrv3, indrv4, indrv5;
40 extern doublereal dlamch_(char *);
41 extern /* Subroutine */ int dlagtf_(integer *, doublereal *, doublereal *,
42 doublereal *, doublereal *, doublereal *, doublereal *, integer *
44 extern integer idamax_(integer *, doublereal *, integer *);
45 extern /* Subroutine */ int xerbla_(char *, integer *), dlagts_(
46 integer *, integer *, doublereal *, doublereal *, doublereal *,
47 doublereal *, integer *, doublereal *, doublereal *, integer *);
49 extern /* Subroutine */ int dlarnv_(integer *, integer *, integer *,
52 doublereal onenrm, dtpcrt, pertol;
55 /* -- LAPACK routine (version 3.1) -- */
56 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
59 /* .. Scalar Arguments .. */
61 /* .. Array Arguments .. */
67 /* DSTEIN computes the eigenvectors of a real symmetric tridiagonal */
68 /* matrix T corresponding to specified eigenvalues, using inverse */
71 /* The maximum number of iterations allowed for each eigenvector is */
72 /* specified by an internal parameter MAXITS (currently set to 5). */
77 /* N (input) INTEGER */
78 /* The order of the matrix. N >= 0. */
80 /* D (input) DOUBLE PRECISION array, dimension (N) */
81 /* The n diagonal elements of the tridiagonal matrix T. */
83 /* E (input) DOUBLE PRECISION array, dimension (N-1) */
84 /* The (n-1) subdiagonal elements of the tridiagonal matrix */
85 /* T, in elements 1 to N-1. */
87 /* M (input) INTEGER */
88 /* The number of eigenvectors to be found. 0 <= M <= N. */
90 /* W (input) DOUBLE PRECISION array, dimension (N) */
91 /* The first M elements of W contain the eigenvalues for */
92 /* which eigenvectors are to be computed. The eigenvalues */
93 /* should be grouped by split-off block and ordered from */
94 /* smallest to largest within the block. ( The output array */
95 /* W from DSTEBZ with ORDER = 'B' is expected here. ) */
97 /* IBLOCK (input) INTEGER array, dimension (N) */
98 /* The submatrix indices associated with the corresponding */
99 /* eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to */
100 /* the first submatrix from the top, =2 if W(i) belongs to */
101 /* the second submatrix, etc. ( The output array IBLOCK */
102 /* from DSTEBZ is expected here. ) */
104 /* ISPLIT (input) INTEGER array, dimension (N) */
105 /* The splitting points, at which T breaks up into submatrices. */
106 /* The first submatrix consists of rows/columns 1 to */
107 /* ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 */
108 /* through ISPLIT( 2 ), etc. */
109 /* ( The output array ISPLIT from DSTEBZ is expected here. ) */
111 /* Z (output) DOUBLE PRECISION 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. */
117 /* LDZ (input) INTEGER */
118 /* The leading dimension of the array Z. LDZ >= max(1,N). */
120 /* WORK (workspace) DOUBLE PRECISION array, dimension (5*N) */
122 /* IWORK (workspace) INTEGER array, dimension (N) */
124 /* IFAIL (output) INTEGER array, dimension (M) */
125 /* On normal exit, all elements of IFAIL are zero. */
126 /* If one or more eigenvectors fail to converge after */
127 /* MAXITS iterations, then their indices are stored in */
130 /* INFO (output) INTEGER */
131 /* = 0: successful exit. */
132 /* < 0: if INFO = -i, the i-th argument had an illegal value */
133 /* > 0: if INFO = i, then i eigenvectors failed to converge */
134 /* in MAXITS iterations. Their indices are stored in */
137 /* Internal Parameters */
138 /* =================== */
140 /* MAXITS INTEGER, default = 5 */
141 /* The maximum number of iterations performed. */
143 /* EXTRA INTEGER, default = 2 */
144 /* The number of iterations performed after norm growth */
145 /* criterion is satisfied, should be at least 1. */
147 /* ===================================================================== */
149 /* .. Parameters .. */
151 /* .. Local Scalars .. */
153 /* .. Local Arrays .. */
155 /* .. External Functions .. */
157 /* .. External Subroutines .. */
159 /* .. Intrinsic Functions .. */
161 /* .. Executable Statements .. */
163 /* Test the input parameters. */
165 /* Parameter adjustments */
172 z_offset = 1 + z_dim1;
181 for (i__ = 1; i__ <= i__1; ++i__) {
188 } else if (*m < 0 || *m > *n) {
190 } else if (*ldz < max(1,*n)) {
194 for (j = 2; j <= i__1; ++j) {
195 if (iblock[j] < iblock[j - 1]) {
199 if (iblock[j] == iblock[j - 1] && w[j] < w[j - 1]) {
211 xerbla_("DSTEIN", &i__1);
215 /* Quick return if possible */
217 if (*n == 0 || *m == 0) {
219 } else if (*n == 1) {
220 z__[z_dim1 + 1] = 1.;
224 /* Get machine constants. */
226 eps = dlamch_("Precision");
228 /* Initialize seed for random number generator DLARNV. */
230 for (i__ = 1; i__ <= 4; ++i__) {
235 /* Initialize pointers. */
238 indrv2 = indrv1 + *n;
239 indrv3 = indrv2 + *n;
240 indrv4 = indrv3 + *n;
241 indrv5 = indrv4 + *n;
243 /* Compute eigenvectors of matrix blocks. */
247 for (nblk = 1; nblk <= i__1; ++nblk) {
249 /* Find starting and ending indices of block nblk. */
254 b1 = isplit[nblk - 1] + 1;
257 blksiz = bn - b1 + 1;
263 /* Compute reorthogonalization criterion and stopping criterion. */
265 onenrm = (d__1 = d__[b1], abs(d__1)) + (d__2 = e[b1], abs(d__2));
267 d__3 = onenrm, d__4 = (d__1 = d__[bn], abs(d__1)) + (d__2 = e[bn - 1],
269 onenrm = max(d__3,d__4);
271 for (i__ = b1 + 1; i__ <= i__2; ++i__) {
273 d__4 = onenrm, d__5 = (d__1 = d__[i__], abs(d__1)) + (d__2 = e[
274 i__ - 1], abs(d__2)) + (d__3 = e[i__], abs(d__3));
275 onenrm = max(d__4,d__5);
278 ortol = onenrm * .001;
280 dtpcrt = sqrt(.1 / blksiz);
282 /* Loop through eigenvalues of block nblk. */
287 for (j = j1; j <= i__2; ++j) {
288 if (iblock[j] != nblk) {
295 /* Skip all the work if the block size is one. */
298 work[indrv1 + 1] = 1.;
302 /* If eigenvalues j and j-1 are too close, add a relatively */
303 /* small perturbation. */
306 eps1 = (d__1 = eps * xj, abs(d__1));
317 /* Get random starting vector. */
319 dlarnv_(&c__2, iseed, &blksiz, &work[indrv1 + 1]);
321 /* Copy the matrix T so it won't be destroyed in factorization. */
323 dcopy_(&blksiz, &d__[b1], &c__1, &work[indrv4 + 1], &c__1);
325 dcopy_(&i__3, &e[b1], &c__1, &work[indrv2 + 2], &c__1);
327 dcopy_(&i__3, &e[b1], &c__1, &work[indrv3 + 1], &c__1);
329 /* Compute LU factors with partial pivoting ( PT = LU ) */
332 dlagtf_(&blksiz, &work[indrv4 + 1], &xj, &work[indrv2 + 2], &work[
333 indrv3 + 1], &tol, &work[indrv5 + 1], &iwork[1], &iinfo);
335 /* Update iteration count. */
343 /* Normalize and scale the righthand side vector Pb. */
346 d__2 = eps, d__3 = (d__1 = work[indrv4 + blksiz], abs(d__1));
347 scl = blksiz * onenrm * max(d__2,d__3) / dasum_(&blksiz, &work[
349 dscal_(&blksiz, &scl, &work[indrv1 + 1], &c__1);
351 /* Solve the system LU = Pb. */
353 dlagts_(&c_n1, &blksiz, &work[indrv4 + 1], &work[indrv2 + 2], &
354 work[indrv3 + 1], &work[indrv5 + 1], &iwork[1], &work[
355 indrv1 + 1], &tol, &iinfo);
357 /* Reorthogonalize by modified Gram-Schmidt if eigenvalues are */
363 if ((d__1 = xj - xjm, abs(d__1)) > ortol) {
368 for (i__ = gpind; i__ <= i__3; ++i__) {
369 ztr = -ddot_(&blksiz, &work[indrv1 + 1], &c__1, &z__[b1 +
370 i__ * z_dim1], &c__1);
371 daxpy_(&blksiz, &ztr, &z__[b1 + i__ * z_dim1], &c__1, &
372 work[indrv1 + 1], &c__1);
377 /* Check the infinity norm of the iterate. */
380 jmax = idamax_(&blksiz, &work[indrv1 + 1], &c__1);
381 nrm = (d__1 = work[indrv1 + jmax], abs(d__1));
383 /* Continue for additional iterations after norm reaches */
384 /* stopping criterion. */
396 /* If stopping criterion was not satisfied, update info and */
397 /* store eigenvector number in array ifail. */
403 /* Accept iterate as jth eigenvector. */
406 scl = 1. / dnrm2_(&blksiz, &work[indrv1 + 1], &c__1);
407 jmax = idamax_(&blksiz, &work[indrv1 + 1], &c__1);
408 if (work[indrv1 + jmax] < 0.) {
411 dscal_(&blksiz, &scl, &work[indrv1 + 1], &c__1);
414 for (i__ = 1; i__ <= i__3; ++i__) {
415 z__[i__ + j * z_dim1] = 0.;
419 for (i__ = 1; i__ <= i__3; ++i__) {
420 z__[b1 + i__ - 1 + j * z_dim1] = work[indrv1 + i__];
424 /* Save the shift to check eigenvalue spacing at next */