1 /* ./src_f77/dstegr.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__1 = 1;
11 static doublereal c_b14 = 0.;
12
dstegr_(char * jobz,char * range,integer * n,doublereal * d__,doublereal * e,doublereal * vl,doublereal * vu,integer * il,integer * iu,doublereal * abstol,integer * m,doublereal * w,doublereal * z__,integer * ldz,integer * isuppz,doublereal * work,integer * lwork,integer * iwork,integer * liwork,integer * info,ftnlen jobz_len,ftnlen range_len)13 /* Subroutine */ int dstegr_(char *jobz, char *range, integer *n, doublereal *
14 d__, doublereal *e, doublereal *vl, doublereal *vu, integer *il,
15 integer *iu, doublereal *abstol, integer *m, doublereal *w,
16 doublereal *z__, integer *ldz, integer *isuppz, doublereal *work,
17 integer *lwork, integer *iwork, integer *liwork, integer *info,
18 ftnlen jobz_len, ftnlen range_len)
19 {
20 /* System generated locals */
21 integer z_dim1, z_offset, i__1, i__2;
22 doublereal d__1, d__2;
23
24 /* Builtin functions */
25 double sqrt(doublereal);
26
27 /* Local variables */
28 static integer i__, j, jj;
29 static doublereal eps, tol, tmp;
30 static integer iend;
31 static doublereal rmin, rmax;
32 static integer itmp;
33 static doublereal tnrm;
34 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
35 integer *);
36 static doublereal scale;
37 extern logical lsame_(char *, char *, ftnlen, ftnlen);
38 static integer iinfo;
39 extern /* Subroutine */ int dswap_(integer *, doublereal *, integer *,
40 doublereal *, integer *);
41 static integer lwmin;
42 static logical wantz;
43 extern doublereal dlamch_(char *, ftnlen);
44 static logical alleig;
45 static integer ibegin;
46 static logical indeig;
47 static integer iindbl;
48 static logical valeig;
49 extern /* Subroutine */ int dlarre_(integer *, doublereal *, doublereal *,
50 doublereal *, integer *, integer *, integer *, doublereal *,
51 doublereal *, doublereal *, doublereal *, integer *), dlaset_(
52 char *, integer *, integer *, doublereal *, doublereal *,
53 doublereal *, integer *, ftnlen);
54 static doublereal safmin;
55 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
56 static doublereal bignum;
57 static integer iindwk, indgrs, indwof;
58 extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *,
59 ftnlen);
60 extern /* Subroutine */ int dlarrv_(integer *, doublereal *, doublereal *,
61 integer *, integer *, doublereal *, integer *, doublereal *,
62 doublereal *, doublereal *, integer *, integer *, doublereal *,
63 integer *, integer *);
64 static doublereal thresh;
65 static integer iinspl, indwrk, liwmin, nsplit;
66 static doublereal smlnum;
67 static logical lquery;
68
69
70 /* -- LAPACK computational routine (version 3.0) -- */
71 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
72 /* Courant Institute, Argonne National Lab, and Rice University */
73 /* June 30, 1999 */
74
75 /* .. Scalar Arguments .. */
76 /* .. */
77 /* .. Array Arguments .. */
78 /* .. */
79
80 /* Purpose */
81 /* ======= */
82
83 /* DSTEGR computes selected eigenvalues and, optionally, eigenvectors */
84 /* of a real symmetric tridiagonal matrix T. Eigenvalues and */
85 /* eigenvectors can be selected by specifying either a range of values */
86 /* or a range of indices for the desired eigenvalues. The eigenvalues */
87 /* are computed by the dqds algorithm, while orthogonal eigenvectors are */
88 /* computed from various ``good'' L D L^T representations (also known as */
89 /* Relatively Robust Representations). Gram-Schmidt orthogonalization is */
90 /* avoided as far as possible. More specifically, the various steps of */
91 /* the algorithm are as follows. For the i-th unreduced block of T, */
92 /* (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T */
93 /* is a relatively robust representation, */
94 /* (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high */
95 /* relative accuracy by the dqds algorithm, */
96 /* (c) If there is a cluster of close eigenvalues, "choose" sigma_i */
97 /* close to the cluster, and go to step (a), */
98 /* (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T, */
99 /* compute the corresponding eigenvector by forming a */
100 /* rank-revealing twisted factorization. */
101 /* The desired accuracy of the output can be specified by the input */
102 /* parameter ABSTOL. */
103
104 /* For more details, see "A new O(n^2) algorithm for the symmetric */
105 /* tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon, */
106 /* Computer Science Division Technical Report No. UCB/CSD-97-971, */
107 /* UC Berkeley, May 1997. */
108
109 /* Note 1 : Currently DSTEGR is only set up to find ALL the n */
110 /* eigenvalues and eigenvectors of T in O(n^2) time */
111 /* Note 2 : Currently the routine DSTEIN is called when an appropriate */
112 /* sigma_i cannot be chosen in step (c) above. DSTEIN invokes modified */
113 /* Gram-Schmidt when eigenvalues are close. */
114 /* Note 3 : DSTEGR works only on machines which follow ieee-754 */
115 /* floating-point standard in their handling of infinities and NaNs. */
116 /* Normal execution of DSTEGR may create NaNs and infinities and hence */
117 /* may abort due to a floating point exception in environments which */
118 /* do not conform to the ieee standard. */
119
120 /* Arguments */
121 /* ========= */
122
123 /* JOBZ (input) CHARACTER*1 */
124 /* = 'N': Compute eigenvalues only; */
125 /* = 'V': Compute eigenvalues and eigenvectors. */
126
127 /* RANGE (input) CHARACTER*1 */
128 /* = 'A': all eigenvalues will be found. */
129 /* = 'V': all eigenvalues in the half-open interval (VL,VU] */
130 /* will be found. */
131 /* = 'I': the IL-th through IU-th eigenvalues will be found. */
132 /* ********* Only RANGE = 'A' is currently supported ********************* */
133
134 /* N (input) INTEGER */
135 /* The order of the matrix. N >= 0. */
136
137 /* D (input/output) DOUBLE PRECISION array, dimension (N) */
138 /* On entry, the n diagonal elements of the tridiagonal matrix */
139 /* T. On exit, D is overwritten. */
140
141 /* E (input/output) DOUBLE PRECISION array, dimension (N) */
142 /* On entry, the (n-1) subdiagonal elements of the tridiagonal */
143 /* matrix T in elements 1 to N-1 of E; E(N) need not be set. */
144 /* On exit, E is overwritten. */
145
146 /* VL (input) DOUBLE PRECISION */
147 /* VU (input) DOUBLE PRECISION */
148 /* If RANGE='V', the lower and upper bounds of the interval to */
149 /* be searched for eigenvalues. VL < VU. */
150 /* Not referenced if RANGE = 'A' or 'I'. */
151
152 /* IL (input) INTEGER */
153 /* IU (input) INTEGER */
154 /* If RANGE='I', the indices (in ascending order) of the */
155 /* smallest and largest eigenvalues to be returned. */
156 /* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */
157 /* Not referenced if RANGE = 'A' or 'V'. */
158
159 /* ABSTOL (input) DOUBLE PRECISION */
160 /* The absolute error tolerance for the */
161 /* eigenvalues/eigenvectors. IF JOBZ = 'V', the eigenvalues and */
162 /* eigenvectors output have residual norms bounded by ABSTOL, */
163 /* and the dot products between different eigenvectors are */
164 /* bounded by ABSTOL. If ABSTOL is less than N*EPS*|T|, then */
165 /* N*EPS*|T| will be used in its place, where EPS is the */
166 /* machine precision and |T| is the 1-norm of the tridiagonal */
167 /* matrix. The eigenvalues are computed to an accuracy of */
168 /* EPS*|T| irrespective of ABSTOL. If high relative accuracy */
169 /* is important, set ABSTOL to DLAMCH( 'Safe minimum' ). */
170 /* See Barlow and Demmel "Computing Accurate Eigensystems of */
171 /* Scaled Diagonally Dominant Matrices", LAPACK Working Note #7 */
172 /* for a discussion of which matrices define their eigenvalues */
173 /* to high relative accuracy. */
174
175 /* M (output) INTEGER */
176 /* The total number of eigenvalues found. 0 <= M <= N. */
177 /* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. */
178
179 /* W (output) DOUBLE PRECISION array, dimension (N) */
180 /* The first M elements contain the selected eigenvalues in */
181 /* ascending order. */
182
183 /* Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) ) */
184 /* If JOBZ = 'V', then if INFO = 0, the first M columns of Z */
185 /* contain the orthonormal eigenvectors of the matrix T */
186 /* corresponding to the selected eigenvalues, with the i-th */
187 /* column of Z holding the eigenvector associated with W(i). */
188 /* If JOBZ = 'N', then Z is not referenced. */
189 /* Note: the user must ensure that at least max(1,M) columns are */
190 /* supplied in the array Z; if RANGE = 'V', the exact value of M */
191 /* is not known in advance and an upper bound must be used. */
192
193 /* LDZ (input) INTEGER */
194 /* The leading dimension of the array Z. LDZ >= 1, and if */
195 /* JOBZ = 'V', LDZ >= max(1,N). */
196
197 /* ISUPPZ (output) INTEGER ARRAY, dimension ( 2*max(1,M) ) */
198 /* The support of the eigenvectors in Z, i.e., the indices */
199 /* indicating the nonzero elements in Z. The i-th eigenvector */
200 /* is nonzero only in elements ISUPPZ( 2*i-1 ) through */
201 /* ISUPPZ( 2*i ). */
202
203 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) */
204 /* On exit, if INFO = 0, WORK(1) returns the optimal */
205 /* (and minimal) LWORK. */
206
207 /* LWORK (input) INTEGER */
208 /* The dimension of the array WORK. LWORK >= max(1,18*N) */
209
210 /* If LWORK = -1, then a workspace query is assumed; the routine */
211 /* only calculates the optimal size of the WORK array, returns */
212 /* this value as the first entry of the WORK array, and no error */
213 /* message related to LWORK is issued by XERBLA. */
214
215 /* IWORK (workspace/output) INTEGER array, dimension (LIWORK) */
216 /* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. */
217
218 /* LIWORK (input) INTEGER */
219 /* The dimension of the array IWORK. LIWORK >= max(1,10*N) */
220
221 /* If LIWORK = -1, then a workspace query is assumed; the */
222 /* routine only calculates the optimal size of the IWORK array, */
223 /* returns this value as the first entry of the IWORK array, and */
224 /* no error message related to LIWORK is issued by XERBLA. */
225
226 /* INFO (output) INTEGER */
227 /* = 0: successful exit */
228 /* < 0: if INFO = -i, the i-th argument had an illegal value */
229 /* > 0: if INFO = 1, internal error in DLARRE, */
230 /* if INFO = 2, internal error in DLARRV. */
231
232 /* Further Details */
233 /* =============== */
234
235 /* Based on contributions by */
236 /* Inderjit Dhillon, IBM Almaden, USA */
237 /* Osni Marques, LBNL/NERSC, USA */
238
239 /* ===================================================================== */
240
241 /* .. Parameters .. */
242 /* .. */
243 /* .. Local Scalars .. */
244 /* .. */
245 /* .. External Functions .. */
246 /* .. */
247 /* .. External Subroutines .. */
248 /* .. */
249 /* .. Intrinsic Functions .. */
250 /* .. */
251 /* .. Executable Statements .. */
252
253 /* Test the input parameters. */
254
255 /* Parameter adjustments */
256 --d__;
257 --e;
258 --w;
259 z_dim1 = *ldz;
260 z_offset = 1 + z_dim1;
261 z__ -= z_offset;
262 --isuppz;
263 --work;
264 --iwork;
265
266 /* Function Body */
267 wantz = lsame_(jobz, "V", (ftnlen)1, (ftnlen)1);
268 alleig = lsame_(range, "A", (ftnlen)1, (ftnlen)1);
269 valeig = lsame_(range, "V", (ftnlen)1, (ftnlen)1);
270 indeig = lsame_(range, "I", (ftnlen)1, (ftnlen)1);
271
272 lquery = *lwork == -1 || *liwork == -1;
273 lwmin = *n * 18;
274 liwmin = *n * 10;
275
276 *info = 0;
277 if (! (wantz || lsame_(jobz, "N", (ftnlen)1, (ftnlen)1))) {
278 *info = -1;
279 } else if (! (alleig || valeig || indeig)) {
280 *info = -2;
281
282 /* The following two lines need to be removed once the */
283 /* RANGE = 'V' and RANGE = 'I' options are provided. */
284
285 } else if (valeig || indeig) {
286 *info = -2;
287 } else if (*n < 0) {
288 *info = -3;
289 } else if (valeig && *n > 0 && *vu <= *vl) {
290 *info = -7;
291 } else if (indeig && *il < 1) {
292 *info = -8;
293 /* The following change should be made in DSTEVX also, otherwise */
294 /* IL can be specified as N+1 and IU as N. */
295 /* ELSE IF( INDEIG .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN */
296 } else if (indeig && (*iu < *il || *iu > *n)) {
297 *info = -9;
298 } else if (*ldz < 1 || wantz && *ldz < *n) {
299 *info = -14;
300 } else if (*lwork < lwmin && ! lquery) {
301 *info = -17;
302 } else if (*liwork < liwmin && ! lquery) {
303 *info = -19;
304 }
305 if (*info == 0) {
306 work[1] = (doublereal) lwmin;
307 iwork[1] = liwmin;
308 }
309
310 if (*info != 0) {
311 i__1 = -(*info);
312 xerbla_("DSTEGR", &i__1, (ftnlen)6);
313 return 0;
314 } else if (lquery) {
315 return 0;
316 }
317
318 /* Quick return if possible */
319
320 *m = 0;
321 if (*n == 0) {
322 return 0;
323 }
324
325 if (*n == 1) {
326 if (alleig || indeig) {
327 *m = 1;
328 w[1] = d__[1];
329 } else {
330 if (*vl < d__[1] && *vu >= d__[1]) {
331 *m = 1;
332 w[1] = d__[1];
333 }
334 }
335 if (wantz) {
336 z__[z_dim1 + 1] = 1.;
337 }
338 return 0;
339 }
340
341 /* Get machine constants. */
342
343 safmin = dlamch_("Safe minimum", (ftnlen)12);
344 eps = dlamch_("Precision", (ftnlen)9);
345 smlnum = safmin / eps;
346 bignum = 1. / smlnum;
347 rmin = sqrt(smlnum);
348 /* Computing MIN */
349 d__1 = sqrt(bignum), d__2 = 1. / sqrt(sqrt(safmin));
350 rmax = min(d__1,d__2);
351
352 /* Scale matrix to allowable range, if necessary. */
353
354 scale = 1.;
355 tnrm = dlanst_("M", n, &d__[1], &e[1], (ftnlen)1);
356 if (tnrm > 0. && tnrm < rmin) {
357 scale = rmin / tnrm;
358 } else if (tnrm > rmax) {
359 scale = rmax / tnrm;
360 }
361 if (scale != 1.) {
362 dscal_(n, &scale, &d__[1], &c__1);
363 i__1 = *n - 1;
364 dscal_(&i__1, &scale, &e[1], &c__1);
365 tnrm *= scale;
366 }
367 indgrs = 1;
368 indwof = (*n << 1) + 1;
369 indwrk = *n * 3 + 1;
370
371 iinspl = 1;
372 iindbl = *n + 1;
373 iindwk = (*n << 1) + 1;
374
375 dlaset_("Full", n, n, &c_b14, &c_b14, &z__[z_offset], ldz, (ftnlen)4);
376
377 /* Compute the desired eigenvalues of the tridiagonal after splitting */
378 /* into smaller subblocks if the corresponding of-diagonal elements */
379 /* are small */
380
381 thresh = eps * tnrm;
382 dlarre_(n, &d__[1], &e[1], &thresh, &nsplit, &iwork[iinspl], m, &w[1], &
383 work[indwof], &work[indgrs], &work[indwrk], &iinfo);
384 if (iinfo != 0) {
385 *info = 1;
386 return 0;
387 }
388
389 if (wantz) {
390
391 /* Compute the desired eigenvectors corresponding to the computed */
392 /* eigenvalues */
393
394 /* Computing MAX */
395 d__1 = *abstol, d__2 = (doublereal) (*n) * thresh;
396 tol = max(d__1,d__2);
397 ibegin = 1;
398 i__1 = nsplit;
399 for (i__ = 1; i__ <= i__1; ++i__) {
400 iend = iwork[iinspl + i__ - 1];
401 i__2 = iend;
402 for (j = ibegin; j <= i__2; ++j) {
403 iwork[iindbl + j - 1] = i__;
404 /* L10: */
405 }
406 ibegin = iend + 1;
407 /* L20: */
408 }
409
410 dlarrv_(n, &d__[1], &e[1], &iwork[iinspl], m, &w[1], &iwork[iindbl], &
411 work[indgrs], &tol, &z__[z_offset], ldz, &isuppz[1], &work[
412 indwrk], &iwork[iindwk], &iinfo);
413 if (iinfo != 0) {
414 *info = 2;
415 return 0;
416 }
417
418 }
419
420 ibegin = 1;
421 i__1 = nsplit;
422 for (i__ = 1; i__ <= i__1; ++i__) {
423 iend = iwork[iinspl + i__ - 1];
424 i__2 = iend;
425 for (j = ibegin; j <= i__2; ++j) {
426 w[j] += work[indwof + i__ - 1];
427 /* L30: */
428 }
429 ibegin = iend + 1;
430 /* L40: */
431 }
432
433 /* If matrix was scaled, then rescale eigenvalues appropriately. */
434
435 if (scale != 1.) {
436 d__1 = 1. / scale;
437 dscal_(m, &d__1, &w[1], &c__1);
438 }
439
440 /* If eigenvalues are not in order, then sort them, along with */
441 /* eigenvectors. */
442
443 if (nsplit > 1) {
444 i__1 = *m - 1;
445 for (j = 1; j <= i__1; ++j) {
446 i__ = 0;
447 tmp = w[j];
448 i__2 = *m;
449 for (jj = j + 1; jj <= i__2; ++jj) {
450 if (w[jj] < tmp) {
451 i__ = jj;
452 tmp = w[jj];
453 }
454 /* L50: */
455 }
456 if (i__ != 0) {
457 w[i__] = w[j];
458 w[j] = tmp;
459 if (wantz) {
460 dswap_(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j * z_dim1
461 + 1], &c__1);
462 itmp = isuppz[(i__ << 1) - 1];
463 isuppz[(i__ << 1) - 1] = isuppz[(j << 1) - 1];
464 isuppz[(j << 1) - 1] = itmp;
465 itmp = isuppz[i__ * 2];
466 isuppz[i__ * 2] = isuppz[j * 2];
467 isuppz[j * 2] = itmp;
468 }
469 }
470 /* L60: */
471 }
472 }
473
474 work[1] = (doublereal) lwmin;
475 iwork[1] = liwmin;
476 return 0;
477
478 /* End of DSTEGR */
479
480 } /* dstegr_ */
481
482