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