1 /* ./src_f77/chbevx.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 complex c_b1 = {0.f,0.f};
11 static complex c_b2 = {1.f,0.f};
12 static real c_b16 = 1.f;
13 static integer c__1 = 1;
14
chbevx_(char * jobz,char * range,char * uplo,integer * n,integer * kd,complex * ab,integer * ldab,complex * q,integer * ldq,real * vl,real * vu,integer * il,integer * iu,real * abstol,integer * m,real * w,complex * z__,integer * ldz,complex * work,real * rwork,integer * iwork,integer * ifail,integer * info,ftnlen jobz_len,ftnlen range_len,ftnlen uplo_len)15 /* Subroutine */ int chbevx_(char *jobz, char *range, char *uplo, integer *n,
16 integer *kd, complex *ab, integer *ldab, complex *q, integer *ldq,
17 real *vl, real *vu, integer *il, integer *iu, real *abstol, integer *
18 m, real *w, complex *z__, integer *ldz, complex *work, real *rwork,
19 integer *iwork, integer *ifail, integer *info, ftnlen jobz_len,
20 ftnlen range_len, ftnlen uplo_len)
21 {
22 /* System generated locals */
23 integer ab_dim1, ab_offset, q_dim1, q_offset, z_dim1, z_offset, i__1,
24 i__2;
25 real r__1, r__2;
26
27 /* Builtin functions */
28 double sqrt(doublereal);
29
30 /* Local variables */
31 static integer i__, j, jj;
32 static real eps, vll, vuu, tmp1;
33 static integer indd, inde;
34 static real anrm;
35 static integer imax;
36 static real rmin, rmax;
37 static complex ctmp1;
38 static integer itmp1, indee;
39 static real sigma;
40 extern logical lsame_(char *, char *, ftnlen, ftnlen);
41 extern /* Subroutine */ int cgemv_(char *, integer *, integer *, complex *
42 , complex *, integer *, complex *, integer *, complex *, complex *
43 , integer *, ftnlen);
44 static integer iinfo;
45 extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
46 static char order[1];
47 extern /* Subroutine */ int ccopy_(integer *, complex *, integer *,
48 complex *, integer *), cswap_(integer *, complex *, integer *,
49 complex *, integer *);
50 static logical lower;
51 extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
52 integer *);
53 static logical wantz;
54 extern doublereal clanhb_(char *, char *, integer *, integer *, complex *,
55 integer *, real *, ftnlen, ftnlen);
56 static logical alleig, indeig;
57 static integer iscale, indibl;
58 extern /* Subroutine */ int clascl_(char *, integer *, integer *, real *,
59 real *, integer *, integer *, complex *, integer *, integer *,
60 ftnlen), chbtrd_(char *, char *, integer *, integer *, complex *,
61 integer *, real *, real *, complex *, integer *, complex *,
62 integer *, ftnlen, ftnlen);
63 static logical valeig;
64 extern doublereal slamch_(char *, ftnlen);
65 extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex
66 *, integer *, complex *, integer *, ftnlen);
67 static real safmin;
68 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
69 static real abstll, bignum;
70 static integer indiwk, indisp;
71 extern /* Subroutine */ int cstein_(integer *, real *, real *, integer *,
72 real *, integer *, integer *, complex *, integer *, real *,
73 integer *, integer *, integer *);
74 static integer indrwk, indwrk;
75 extern /* Subroutine */ int csteqr_(char *, integer *, real *, real *,
76 complex *, integer *, real *, integer *, ftnlen), ssterf_(integer
77 *, real *, real *, integer *);
78 static integer nsplit;
79 extern /* Subroutine */ int sstebz_(char *, char *, integer *, real *,
80 real *, integer *, integer *, real *, real *, real *, integer *,
81 integer *, real *, integer *, integer *, real *, integer *,
82 integer *, ftnlen, ftnlen);
83 static real smlnum;
84
85
86 /* -- LAPACK driver routine (version 3.0) -- */
87 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
88 /* Courant Institute, Argonne National Lab, and Rice University */
89 /* June 30, 1999 */
90
91 /* .. Scalar Arguments .. */
92 /* .. */
93 /* .. Array Arguments .. */
94 /* .. */
95
96 /* Purpose */
97 /* ======= */
98
99 /* CHBEVX computes selected eigenvalues and, optionally, eigenvectors */
100 /* of a complex Hermitian band matrix A. Eigenvalues and eigenvectors */
101 /* can be selected by specifying either a range of values or a range of */
102 /* indices for the desired eigenvalues. */
103
104 /* Arguments */
105 /* ========= */
106
107 /* JOBZ (input) CHARACTER*1 */
108 /* = 'N': Compute eigenvalues only; */
109 /* = 'V': Compute eigenvalues and eigenvectors. */
110
111 /* RANGE (input) CHARACTER*1 */
112 /* = 'A': all eigenvalues will be found; */
113 /* = 'V': all eigenvalues in the half-open interval (VL,VU] */
114 /* will be found; */
115 /* = 'I': the IL-th through IU-th eigenvalues will be found. */
116
117 /* UPLO (input) CHARACTER*1 */
118 /* = 'U': Upper triangle of A is stored; */
119 /* = 'L': Lower triangle of A is stored. */
120
121 /* N (input) INTEGER */
122 /* The order of the matrix A. N >= 0. */
123
124 /* KD (input) INTEGER */
125 /* The number of superdiagonals of the matrix A if UPLO = 'U', */
126 /* or the number of subdiagonals if UPLO = 'L'. KD >= 0. */
127
128 /* AB (input/output) COMPLEX array, dimension (LDAB, N) */
129 /* On entry, the upper or lower triangle of the Hermitian band */
130 /* matrix A, stored in the first KD+1 rows of the array. The */
131 /* j-th column of A is stored in the j-th column of the array AB */
132 /* as follows: */
133 /* if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; */
134 /* if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). */
135
136 /* On exit, AB is overwritten by values generated during the */
137 /* reduction to tridiagonal form. */
138
139 /* LDAB (input) INTEGER */
140 /* The leading dimension of the array AB. LDAB >= KD + 1. */
141
142 /* Q (output) COMPLEX array, dimension (LDQ, N) */
143 /* If JOBZ = 'V', the N-by-N unitary matrix used in the */
144 /* reduction to tridiagonal form. */
145 /* If JOBZ = 'N', the array Q is not referenced. */
146
147 /* LDQ (input) INTEGER */
148 /* The leading dimension of the array Q. If JOBZ = 'V', then */
149 /* LDQ >= max(1,N). */
150
151 /* VL (input) REAL */
152 /* VU (input) REAL */
153 /* If RANGE='V', the lower and upper bounds of the interval to */
154 /* be searched for eigenvalues. VL < VU. */
155 /* Not referenced if RANGE = 'A' or 'I'. */
156
157 /* IL (input) INTEGER */
158 /* IU (input) INTEGER */
159 /* If RANGE='I', the indices (in ascending order) of the */
160 /* smallest and largest eigenvalues to be returned. */
161 /* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */
162 /* Not referenced if RANGE = 'A' or 'V'. */
163
164 /* ABSTOL (input) REAL */
165 /* The absolute error tolerance for the eigenvalues. */
166 /* An approximate eigenvalue is accepted as converged */
167 /* when it is determined to lie in an interval [a,b] */
168 /* of width less than or equal to */
169
170 /* ABSTOL + EPS * max( |a|,|b| ) , */
171
172 /* where EPS is the machine precision. If ABSTOL is less than */
173 /* or equal to zero, then EPS*|T| will be used in its place, */
174 /* where |T| is the 1-norm of the tridiagonal matrix obtained */
175 /* by reducing AB to tridiagonal form. */
176
177 /* Eigenvalues will be computed most accurately when ABSTOL is */
178 /* set to twice the underflow threshold 2*SLAMCH('S'), not zero. */
179 /* If this routine returns with INFO>0, indicating that some */
180 /* eigenvectors did not converge, try setting ABSTOL to */
181 /* 2*SLAMCH('S'). */
182
183 /* See "Computing Small Singular Values of Bidiagonal Matrices */
184 /* with Guaranteed High Relative Accuracy," by Demmel and */
185 /* Kahan, LAPACK Working Note #3. */
186
187 /* M (output) INTEGER */
188 /* The total number of eigenvalues found. 0 <= M <= N. */
189 /* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. */
190
191 /* W (output) REAL array, dimension (N) */
192 /* The first M elements contain the selected eigenvalues in */
193 /* ascending order. */
194
195 /* Z (output) COMPLEX array, dimension (LDZ, max(1,M)) */
196 /* If JOBZ = 'V', then if INFO = 0, the first M columns of Z */
197 /* contain the orthonormal eigenvectors of the matrix A */
198 /* corresponding to the selected eigenvalues, with the i-th */
199 /* column of Z holding the eigenvector associated with W(i). */
200 /* If an eigenvector fails to converge, then that column of Z */
201 /* contains the latest approximation to the eigenvector, and the */
202 /* index of the eigenvector is returned in IFAIL. */
203 /* If JOBZ = 'N', then Z is not referenced. */
204 /* Note: the user must ensure that at least max(1,M) columns are */
205 /* supplied in the array Z; if RANGE = 'V', the exact value of M */
206 /* is not known in advance and an upper bound must be used. */
207
208 /* LDZ (input) INTEGER */
209 /* The leading dimension of the array Z. LDZ >= 1, and if */
210 /* JOBZ = 'V', LDZ >= max(1,N). */
211
212 /* WORK (workspace) COMPLEX array, dimension (N) */
213
214 /* RWORK (workspace) REAL array, dimension (7*N) */
215
216 /* IWORK (workspace) INTEGER array, dimension (5*N) */
217
218 /* IFAIL (output) INTEGER array, dimension (N) */
219 /* If JOBZ = 'V', then if INFO = 0, the first M elements of */
220 /* IFAIL are zero. If INFO > 0, then IFAIL contains the */
221 /* indices of the eigenvectors that failed to converge. */
222 /* If JOBZ = 'N', then IFAIL is not referenced. */
223
224 /* INFO (output) INTEGER */
225 /* = 0: successful exit */
226 /* < 0: if INFO = -i, the i-th argument had an illegal value */
227 /* > 0: if INFO = i, then i eigenvectors failed to converge. */
228 /* Their indices are stored in array IFAIL. */
229
230 /* ===================================================================== */
231
232 /* .. Parameters .. */
233 /* .. */
234 /* .. Local Scalars .. */
235 /* .. */
236 /* .. External Functions .. */
237 /* .. */
238 /* .. External Subroutines .. */
239 /* .. */
240 /* .. Intrinsic Functions .. */
241 /* .. */
242 /* .. Executable Statements .. */
243
244 /* Test the input parameters. */
245
246 /* Parameter adjustments */
247 ab_dim1 = *ldab;
248 ab_offset = 1 + ab_dim1;
249 ab -= ab_offset;
250 q_dim1 = *ldq;
251 q_offset = 1 + q_dim1;
252 q -= q_offset;
253 --w;
254 z_dim1 = *ldz;
255 z_offset = 1 + z_dim1;
256 z__ -= z_offset;
257 --work;
258 --rwork;
259 --iwork;
260 --ifail;
261
262 /* Function Body */
263 wantz = lsame_(jobz, "V", (ftnlen)1, (ftnlen)1);
264 alleig = lsame_(range, "A", (ftnlen)1, (ftnlen)1);
265 valeig = lsame_(range, "V", (ftnlen)1, (ftnlen)1);
266 indeig = lsame_(range, "I", (ftnlen)1, (ftnlen)1);
267 lower = lsame_(uplo, "L", (ftnlen)1, (ftnlen)1);
268
269 *info = 0;
270 if (! (wantz || lsame_(jobz, "N", (ftnlen)1, (ftnlen)1))) {
271 *info = -1;
272 } else if (! (alleig || valeig || indeig)) {
273 *info = -2;
274 } else if (! (lower || lsame_(uplo, "U", (ftnlen)1, (ftnlen)1))) {
275 *info = -3;
276 } else if (*n < 0) {
277 *info = -4;
278 } else if (*kd < 0) {
279 *info = -5;
280 } else if (*ldab < *kd + 1) {
281 *info = -7;
282 } else if (wantz && *ldq < max(1,*n)) {
283 *info = -9;
284 } else {
285 if (valeig) {
286 if (*n > 0 && *vu <= *vl) {
287 *info = -11;
288 }
289 } else if (indeig) {
290 if (*il < 1 || *il > max(1,*n)) {
291 *info = -12;
292 } else if (*iu < min(*n,*il) || *iu > *n) {
293 *info = -13;
294 }
295 }
296 }
297 if (*info == 0) {
298 if (*ldz < 1 || wantz && *ldz < *n) {
299 *info = -18;
300 }
301 }
302
303 if (*info != 0) {
304 i__1 = -(*info);
305 xerbla_("CHBEVX", &i__1, (ftnlen)6);
306 return 0;
307 }
308
309 /* Quick return if possible */
310
311 *m = 0;
312 if (*n == 0) {
313 return 0;
314 }
315
316 if (*n == 1) {
317 *m = 1;
318 if (lower) {
319 i__1 = ab_dim1 + 1;
320 ctmp1.r = ab[i__1].r, ctmp1.i = ab[i__1].i;
321 } else {
322 i__1 = *kd + 1 + ab_dim1;
323 ctmp1.r = ab[i__1].r, ctmp1.i = ab[i__1].i;
324 }
325 tmp1 = ctmp1.r;
326 if (valeig) {
327 if (! (*vl < tmp1 && *vu >= tmp1)) {
328 *m = 0;
329 }
330 }
331 if (*m == 1) {
332 w[1] = ctmp1.r;
333 if (wantz) {
334 i__1 = z_dim1 + 1;
335 z__[i__1].r = 1.f, z__[i__1].i = 0.f;
336 }
337 }
338 return 0;
339 }
340
341 /* Get machine constants. */
342
343 safmin = slamch_("Safe minimum", (ftnlen)12);
344 eps = slamch_("Precision", (ftnlen)9);
345 smlnum = safmin / eps;
346 bignum = 1.f / smlnum;
347 rmin = sqrt(smlnum);
348 /* Computing MIN */
349 r__1 = sqrt(bignum), r__2 = 1.f / sqrt(sqrt(safmin));
350 rmax = dmin(r__1,r__2);
351
352 /* Scale matrix to allowable range, if necessary. */
353
354 iscale = 0;
355 abstll = *abstol;
356 if (valeig) {
357 vll = *vl;
358 vuu = *vu;
359 } else {
360 vll = 0.f;
361 vuu = 0.f;
362 }
363 anrm = clanhb_("M", uplo, n, kd, &ab[ab_offset], ldab, &rwork[1], (ftnlen)
364 1, (ftnlen)1);
365 if (anrm > 0.f && anrm < rmin) {
366 iscale = 1;
367 sigma = rmin / anrm;
368 } else if (anrm > rmax) {
369 iscale = 1;
370 sigma = rmax / anrm;
371 }
372 if (iscale == 1) {
373 if (lower) {
374 clascl_("B", kd, kd, &c_b16, &sigma, n, n, &ab[ab_offset], ldab,
375 info, (ftnlen)1);
376 } else {
377 clascl_("Q", kd, kd, &c_b16, &sigma, n, n, &ab[ab_offset], ldab,
378 info, (ftnlen)1);
379 }
380 if (*abstol > 0.f) {
381 abstll = *abstol * sigma;
382 }
383 if (valeig) {
384 vll = *vl * sigma;
385 vuu = *vu * sigma;
386 }
387 }
388
389 /* Call CHBTRD to reduce Hermitian band matrix to tridiagonal form. */
390
391 indd = 1;
392 inde = indd + *n;
393 indrwk = inde + *n;
394 indwrk = 1;
395 chbtrd_(jobz, uplo, n, kd, &ab[ab_offset], ldab, &rwork[indd], &rwork[
396 inde], &q[q_offset], ldq, &work[indwrk], &iinfo, (ftnlen)1, (
397 ftnlen)1);
398
399 /* If all eigenvalues are desired and ABSTOL is less than or equal */
400 /* to zero, then call SSTERF or CSTEQR. If this fails for some */
401 /* eigenvalue, then try SSTEBZ. */
402
403 if ((alleig || indeig && *il == 1 && *iu == *n) && *abstol <= 0.f) {
404 scopy_(n, &rwork[indd], &c__1, &w[1], &c__1);
405 indee = indrwk + (*n << 1);
406 if (! wantz) {
407 i__1 = *n - 1;
408 scopy_(&i__1, &rwork[inde], &c__1, &rwork[indee], &c__1);
409 ssterf_(n, &w[1], &rwork[indee], info);
410 } else {
411 clacpy_("A", n, n, &q[q_offset], ldq, &z__[z_offset], ldz, (
412 ftnlen)1);
413 i__1 = *n - 1;
414 scopy_(&i__1, &rwork[inde], &c__1, &rwork[indee], &c__1);
415 csteqr_(jobz, n, &w[1], &rwork[indee], &z__[z_offset], ldz, &
416 rwork[indrwk], info, (ftnlen)1);
417 if (*info == 0) {
418 i__1 = *n;
419 for (i__ = 1; i__ <= i__1; ++i__) {
420 ifail[i__] = 0;
421 /* L10: */
422 }
423 }
424 }
425 if (*info == 0) {
426 *m = *n;
427 goto L30;
428 }
429 *info = 0;
430 }
431
432 /* Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN. */
433
434 if (wantz) {
435 *(unsigned char *)order = 'B';
436 } else {
437 *(unsigned char *)order = 'E';
438 }
439 indibl = 1;
440 indisp = indibl + *n;
441 indiwk = indisp + *n;
442 sstebz_(range, order, n, &vll, &vuu, il, iu, &abstll, &rwork[indd], &
443 rwork[inde], m, &nsplit, &w[1], &iwork[indibl], &iwork[indisp], &
444 rwork[indrwk], &iwork[indiwk], info, (ftnlen)1, (ftnlen)1);
445
446 if (wantz) {
447 cstein_(n, &rwork[indd], &rwork[inde], m, &w[1], &iwork[indibl], &
448 iwork[indisp], &z__[z_offset], ldz, &rwork[indrwk], &iwork[
449 indiwk], &ifail[1], info);
450
451 /* Apply unitary matrix used in reduction to tridiagonal */
452 /* form to eigenvectors returned by CSTEIN. */
453
454 i__1 = *m;
455 for (j = 1; j <= i__1; ++j) {
456 ccopy_(n, &z__[j * z_dim1 + 1], &c__1, &work[1], &c__1);
457 cgemv_("N", n, n, &c_b2, &q[q_offset], ldq, &work[1], &c__1, &
458 c_b1, &z__[j * z_dim1 + 1], &c__1, (ftnlen)1);
459 /* L20: */
460 }
461 }
462
463 /* If matrix was scaled, then rescale eigenvalues appropriately. */
464
465 L30:
466 if (iscale == 1) {
467 if (*info == 0) {
468 imax = *m;
469 } else {
470 imax = *info - 1;
471 }
472 r__1 = 1.f / sigma;
473 sscal_(&imax, &r__1, &w[1], &c__1);
474 }
475
476 /* If eigenvalues are not in order, then sort them, along with */
477 /* eigenvectors. */
478
479 if (wantz) {
480 i__1 = *m - 1;
481 for (j = 1; j <= i__1; ++j) {
482 i__ = 0;
483 tmp1 = w[j];
484 i__2 = *m;
485 for (jj = j + 1; jj <= i__2; ++jj) {
486 if (w[jj] < tmp1) {
487 i__ = jj;
488 tmp1 = w[jj];
489 }
490 /* L40: */
491 }
492
493 if (i__ != 0) {
494 itmp1 = iwork[indibl + i__ - 1];
495 w[i__] = w[j];
496 iwork[indibl + i__ - 1] = iwork[indibl + j - 1];
497 w[j] = tmp1;
498 iwork[indibl + j - 1] = itmp1;
499 cswap_(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j * z_dim1 + 1],
500 &c__1);
501 if (*info != 0) {
502 itmp1 = ifail[i__];
503 ifail[i__] = ifail[j];
504 ifail[j] = itmp1;
505 }
506 }
507 /* L50: */
508 }
509 }
510
511 return 0;
512
513 /* End of CHBEVX */
514
515 } /* chbevx_ */
516
517