1 /* ./src_f77/cheevd.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__0 = 0;
11 static real c_b13 = 1.f;
12 static integer c__1 = 1;
13
cheevd_(char * jobz,char * uplo,integer * n,complex * a,integer * lda,real * w,complex * work,integer * lwork,real * rwork,integer * lrwork,integer * iwork,integer * liwork,integer * info,ftnlen jobz_len,ftnlen uplo_len)14 /* Subroutine */ int cheevd_(char *jobz, char *uplo, integer *n, complex *a,
15 integer *lda, real *w, complex *work, integer *lwork, real *rwork,
16 integer *lrwork, integer *iwork, integer *liwork, integer *info,
17 ftnlen jobz_len, ftnlen uplo_len)
18 {
19 /* System generated locals */
20 integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
21 real r__1, r__2;
22
23 /* Builtin functions */
24 double sqrt(doublereal);
25
26 /* Local variables */
27 static real eps;
28 static integer inde;
29 static real anrm;
30 static integer imax;
31 static real rmin, rmax;
32 static integer lopt;
33 static real sigma;
34 extern logical lsame_(char *, char *, ftnlen, ftnlen);
35 static integer iinfo;
36 extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
37 static integer lwmin, liopt;
38 static logical lower;
39 static integer llrwk, lropt;
40 static logical wantz;
41 static integer indwk2, llwrk2;
42 extern doublereal clanhe_(char *, char *, integer *, complex *, integer *,
43 real *, ftnlen, ftnlen);
44 static integer iscale;
45 extern /* Subroutine */ int clascl_(char *, integer *, integer *, real *,
46 real *, integer *, integer *, complex *, integer *, integer *,
47 ftnlen), cstedc_(char *, integer *, real *, real *, complex *,
48 integer *, complex *, integer *, real *, integer *, integer *,
49 integer *, integer *, ftnlen);
50 extern doublereal slamch_(char *, ftnlen);
51 extern /* Subroutine */ int chetrd_(char *, integer *, complex *, integer
52 *, real *, real *, complex *, complex *, integer *, integer *,
53 ftnlen), clacpy_(char *, integer *, integer *, complex *, integer
54 *, complex *, integer *, ftnlen);
55 static real safmin;
56 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
57 static real bignum;
58 static integer indtau, indrwk, indwrk, liwmin;
59 extern /* Subroutine */ int ssterf_(integer *, real *, real *, integer *);
60 static integer lrwmin;
61 extern /* Subroutine */ int cunmtr_(char *, char *, char *, integer *,
62 integer *, complex *, integer *, complex *, complex *, integer *,
63 complex *, integer *, integer *, ftnlen, ftnlen, ftnlen);
64 static integer llwork;
65 static real smlnum;
66 static logical lquery;
67
68
69 /* -- LAPACK driver routine (version 3.0) -- */
70 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
71 /* Courant Institute, Argonne National Lab, and Rice University */
72 /* June 30, 1999 */
73
74 /* .. Scalar Arguments .. */
75 /* .. */
76 /* .. Array Arguments .. */
77 /* .. */
78
79 /* Purpose */
80 /* ======= */
81
82 /* CHEEVD computes all eigenvalues and, optionally, eigenvectors of a */
83 /* complex Hermitian matrix A. If eigenvectors are desired, it uses a */
84 /* divide and conquer algorithm. */
85
86 /* The divide and conquer algorithm makes very mild assumptions about */
87 /* floating point arithmetic. It will work on machines with a guard */
88 /* digit in add/subtract, or on those binary machines without guard */
89 /* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or */
90 /* Cray-2. It could conceivably fail on hexadecimal or decimal machines */
91 /* without guard digits, but we know of none. */
92
93 /* Arguments */
94 /* ========= */
95
96 /* JOBZ (input) CHARACTER*1 */
97 /* = 'N': Compute eigenvalues only; */
98 /* = 'V': Compute eigenvalues and eigenvectors. */
99
100 /* UPLO (input) CHARACTER*1 */
101 /* = 'U': Upper triangle of A is stored; */
102 /* = 'L': Lower triangle of A is stored. */
103
104 /* N (input) INTEGER */
105 /* The order of the matrix A. N >= 0. */
106
107 /* A (input/output) COMPLEX array, dimension (LDA, N) */
108 /* On entry, the Hermitian matrix A. If UPLO = 'U', the */
109 /* leading N-by-N upper triangular part of A contains the */
110 /* upper triangular part of the matrix A. If UPLO = 'L', */
111 /* the leading N-by-N lower triangular part of A contains */
112 /* the lower triangular part of the matrix A. */
113 /* On exit, if JOBZ = 'V', then if INFO = 0, A contains the */
114 /* orthonormal eigenvectors of the matrix A. */
115 /* If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') */
116 /* or the upper triangle (if UPLO='U') of A, including the */
117 /* diagonal, is destroyed. */
118
119 /* LDA (input) INTEGER */
120 /* The leading dimension of the array A. LDA >= max(1,N). */
121
122 /* W (output) REAL array, dimension (N) */
123 /* If INFO = 0, the eigenvalues in ascending order. */
124
125 /* WORK (workspace/output) COMPLEX array, dimension (LWORK) */
126 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
127
128 /* LWORK (input) INTEGER */
129 /* The length of the array WORK. */
130 /* If N <= 1, LWORK must be at least 1. */
131 /* If JOBZ = 'N' and N > 1, LWORK must be at least N + 1. */
132 /* If JOBZ = 'V' and N > 1, LWORK must be at least 2*N + N**2. */
133
134 /* If LWORK = -1, then a workspace query is assumed; the routine */
135 /* only calculates the optimal size of the WORK array, returns */
136 /* this value as the first entry of the WORK array, and no error */
137 /* message related to LWORK is issued by XERBLA. */
138
139 /* RWORK (workspace/output) REAL array, */
140 /* dimension (LRWORK) */
141 /* On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. */
142
143 /* LRWORK (input) INTEGER */
144 /* The dimension of the array RWORK. */
145 /* If N <= 1, LRWORK must be at least 1. */
146 /* If JOBZ = 'N' and N > 1, LRWORK must be at least N. */
147 /* If JOBZ = 'V' and N > 1, LRWORK must be at least */
148 /* 1 + 5*N + 2*N**2. */
149
150 /* If LRWORK = -1, then a workspace query is assumed; the */
151 /* routine only calculates the optimal size of the RWORK array, */
152 /* returns this value as the first entry of the RWORK array, and */
153 /* no error message related to LRWORK is issued by XERBLA. */
154
155 /* IWORK (workspace/output) INTEGER array, dimension (LIWORK) */
156 /* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. */
157
158 /* LIWORK (input) INTEGER */
159 /* The dimension of the array IWORK. */
160 /* If N <= 1, LIWORK must be at least 1. */
161 /* If JOBZ = 'N' and N > 1, LIWORK must be at least 1. */
162 /* If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N. */
163
164 /* If LIWORK = -1, then a workspace query is assumed; the */
165 /* routine only calculates the optimal size of the IWORK array, */
166 /* returns this value as the first entry of the IWORK array, and */
167 /* no error message related to LIWORK is issued by XERBLA. */
168
169 /* INFO (output) INTEGER */
170 /* = 0: successful exit */
171 /* < 0: if INFO = -i, the i-th argument had an illegal value */
172 /* > 0: if INFO = i, the algorithm failed to converge; i */
173 /* off-diagonal elements of an intermediate tridiagonal */
174 /* form did not converge to zero. */
175
176 /* Further Details */
177 /* =============== */
178
179 /* Based on contributions by */
180 /* Jeff Rutter, Computer Science Division, University of California */
181 /* at Berkeley, USA */
182
183 /* ===================================================================== */
184
185 /* .. Parameters .. */
186 /* .. */
187 /* .. Local Scalars .. */
188 /* .. */
189 /* .. External Functions .. */
190 /* .. */
191 /* .. External Subroutines .. */
192 /* .. */
193 /* .. Intrinsic Functions .. */
194 /* .. */
195 /* .. Executable Statements .. */
196
197 /* Test the input parameters. */
198
199 /* Parameter adjustments */
200 a_dim1 = *lda;
201 a_offset = 1 + a_dim1;
202 a -= a_offset;
203 --w;
204 --work;
205 --rwork;
206 --iwork;
207
208 /* Function Body */
209 wantz = lsame_(jobz, "V", (ftnlen)1, (ftnlen)1);
210 lower = lsame_(uplo, "L", (ftnlen)1, (ftnlen)1);
211 lquery = *lwork == -1 || *lrwork == -1 || *liwork == -1;
212
213 *info = 0;
214 if (*n <= 1) {
215 lwmin = 1;
216 lrwmin = 1;
217 liwmin = 1;
218 lopt = lwmin;
219 lropt = lrwmin;
220 liopt = liwmin;
221 } else {
222 if (wantz) {
223 lwmin = (*n << 1) + *n * *n;
224 /* Computing 2nd power */
225 i__1 = *n;
226 lrwmin = *n * 5 + 1 + (i__1 * i__1 << 1);
227 liwmin = *n * 5 + 3;
228 } else {
229 lwmin = *n + 1;
230 lrwmin = *n;
231 liwmin = 1;
232 }
233 lopt = lwmin;
234 lropt = lrwmin;
235 liopt = liwmin;
236 }
237 if (! (wantz || lsame_(jobz, "N", (ftnlen)1, (ftnlen)1))) {
238 *info = -1;
239 } else if (! (lower || lsame_(uplo, "U", (ftnlen)1, (ftnlen)1))) {
240 *info = -2;
241 } else if (*n < 0) {
242 *info = -3;
243 } else if (*lda < max(1,*n)) {
244 *info = -5;
245 } else if (*lwork < lwmin && ! lquery) {
246 *info = -8;
247 } else if (*lrwork < lrwmin && ! lquery) {
248 *info = -10;
249 } else if (*liwork < liwmin && ! lquery) {
250 *info = -12;
251 }
252
253 if (*info == 0) {
254 work[1].r = (real) lopt, work[1].i = 0.f;
255 rwork[1] = (real) lropt;
256 iwork[1] = liopt;
257 }
258
259 if (*info != 0) {
260 i__1 = -(*info);
261 xerbla_("CHEEVD", &i__1, (ftnlen)6);
262 return 0;
263 } else if (lquery) {
264 return 0;
265 }
266
267 /* Quick return if possible */
268
269 if (*n == 0) {
270 return 0;
271 }
272
273 if (*n == 1) {
274 i__1 = a_dim1 + 1;
275 w[1] = a[i__1].r;
276 if (wantz) {
277 i__1 = a_dim1 + 1;
278 a[i__1].r = 1.f, a[i__1].i = 0.f;
279 }
280 return 0;
281 }
282
283 /* Get machine constants. */
284
285 safmin = slamch_("Safe minimum", (ftnlen)12);
286 eps = slamch_("Precision", (ftnlen)9);
287 smlnum = safmin / eps;
288 bignum = 1.f / smlnum;
289 rmin = sqrt(smlnum);
290 rmax = sqrt(bignum);
291
292 /* Scale matrix to allowable range, if necessary. */
293
294 anrm = clanhe_("M", uplo, n, &a[a_offset], lda, &rwork[1], (ftnlen)1, (
295 ftnlen)1);
296 iscale = 0;
297 if (anrm > 0.f && anrm < rmin) {
298 iscale = 1;
299 sigma = rmin / anrm;
300 } else if (anrm > rmax) {
301 iscale = 1;
302 sigma = rmax / anrm;
303 }
304 if (iscale == 1) {
305 clascl_(uplo, &c__0, &c__0, &c_b13, &sigma, n, n, &a[a_offset], lda,
306 info, (ftnlen)1);
307 }
308
309 /* Call CHETRD to reduce Hermitian matrix to tridiagonal form. */
310
311 inde = 1;
312 indtau = 1;
313 indwrk = indtau + *n;
314 indrwk = inde + *n;
315 indwk2 = indwrk + *n * *n;
316 llwork = *lwork - indwrk + 1;
317 llwrk2 = *lwork - indwk2 + 1;
318 llrwk = *lrwork - indrwk + 1;
319 chetrd_(uplo, n, &a[a_offset], lda, &w[1], &rwork[inde], &work[indtau], &
320 work[indwrk], &llwork, &iinfo, (ftnlen)1);
321 /* Computing MAX */
322 i__1 = indwrk;
323 r__1 = (real) lopt, r__2 = (real) (*n) + work[i__1].r;
324 lopt = dmax(r__1,r__2);
325
326 /* For eigenvalues only, call SSTERF. For eigenvectors, first call */
327 /* CSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the */
328 /* tridiagonal matrix, then call CUNMTR to multiply it to the */
329 /* Householder transformations represented as Householder vectors in */
330 /* A. */
331
332 if (! wantz) {
333 ssterf_(n, &w[1], &rwork[inde], info);
334 } else {
335 cstedc_("I", n, &w[1], &rwork[inde], &work[indwrk], n, &work[indwk2],
336 &llwrk2, &rwork[indrwk], &llrwk, &iwork[1], liwork, info, (
337 ftnlen)1);
338 cunmtr_("L", uplo, "N", n, n, &a[a_offset], lda, &work[indtau], &work[
339 indwrk], n, &work[indwk2], &llwrk2, &iinfo, (ftnlen)1, (
340 ftnlen)1, (ftnlen)1);
341 clacpy_("A", n, n, &work[indwrk], n, &a[a_offset], lda, (ftnlen)1);
342 /* Computing MAX */
343 /* Computing 2nd power */
344 i__3 = *n;
345 i__4 = indwk2;
346 i__1 = lopt, i__2 = *n + i__3 * i__3 + (integer) work[i__4].r;
347 lopt = max(i__1,i__2);
348 }
349
350 /* If matrix was scaled, then rescale eigenvalues appropriately. */
351
352 if (iscale == 1) {
353 if (*info == 0) {
354 imax = *n;
355 } else {
356 imax = *info - 1;
357 }
358 r__1 = 1.f / sigma;
359 sscal_(&imax, &r__1, &w[1], &c__1);
360 }
361
362 work[1].r = (real) lopt, work[1].i = 0.f;
363 rwork[1] = (real) lropt;
364 iwork[1] = liopt;
365
366 return 0;
367
368 /* End of CHEEVD */
369
370 } /* cheevd_ */
371
372