1 /* ./src_f77/dsbevd.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 doublereal c_b11 = 1.;
11 static doublereal c_b18 = 0.;
12 static integer c__1 = 1;
13
dsbevd_(char * jobz,char * uplo,integer * n,integer * kd,doublereal * ab,integer * ldab,doublereal * w,doublereal * z__,integer * ldz,doublereal * work,integer * lwork,integer * iwork,integer * liwork,integer * info,ftnlen jobz_len,ftnlen uplo_len)14 /* Subroutine */ int dsbevd_(char *jobz, char *uplo, integer *n, integer *kd,
15 doublereal *ab, integer *ldab, doublereal *w, doublereal *z__,
16 integer *ldz, doublereal *work, integer *lwork, integer *iwork,
17 integer *liwork, integer *info, ftnlen jobz_len, ftnlen uplo_len)
18 {
19 /* System generated locals */
20 integer ab_dim1, ab_offset, z_dim1, z_offset, i__1;
21 doublereal d__1;
22
23 /* Builtin functions */
24 double sqrt(doublereal);
25
26 /* Local variables */
27 static doublereal eps;
28 static integer inde;
29 static doublereal anrm, rmin, rmax;
30 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
31 integer *), dgemm_(char *, char *, integer *, integer *, integer *
32 , doublereal *, doublereal *, integer *, doublereal *, integer *,
33 doublereal *, doublereal *, integer *, ftnlen, ftnlen);
34 static doublereal sigma;
35 extern logical lsame_(char *, char *, ftnlen, ftnlen);
36 static integer iinfo, lwmin;
37 static logical lower, wantz;
38 static integer indwk2, llwrk2;
39 extern doublereal dlamch_(char *, ftnlen);
40 static integer iscale;
41 extern /* Subroutine */ int dlascl_(char *, integer *, integer *,
42 doublereal *, doublereal *, integer *, integer *, doublereal *,
43 integer *, integer *, ftnlen);
44 extern doublereal dlansb_(char *, char *, integer *, integer *,
45 doublereal *, integer *, doublereal *, ftnlen, ftnlen);
46 extern /* Subroutine */ int dstedc_(char *, integer *, doublereal *,
47 doublereal *, doublereal *, integer *, doublereal *, integer *,
48 integer *, integer *, integer *, ftnlen), dlacpy_(char *, integer
49 *, integer *, doublereal *, integer *, doublereal *, integer *,
50 ftnlen);
51 static doublereal safmin;
52 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
53 static doublereal bignum;
54 extern /* Subroutine */ int dsbtrd_(char *, char *, integer *, integer *,
55 doublereal *, integer *, doublereal *, doublereal *, doublereal *,
56 integer *, doublereal *, integer *, ftnlen, ftnlen), dsterf_(
57 integer *, doublereal *, doublereal *, integer *);
58 static integer indwrk, liwmin;
59 static doublereal smlnum;
60 static logical lquery;
61
62
63 /* -- LAPACK driver routine (version 3.0) -- */
64 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
65 /* Courant Institute, Argonne National Lab, and Rice University */
66 /* June 30, 1999 */
67
68 /* .. Scalar Arguments .. */
69 /* .. */
70 /* .. Array Arguments .. */
71 /* .. */
72
73 /* Purpose */
74 /* ======= */
75
76 /* DSBEVD computes all the eigenvalues and, optionally, eigenvectors of */
77 /* a real symmetric band matrix A. If eigenvectors are desired, it uses */
78 /* a divide and conquer algorithm. */
79
80 /* The divide and conquer algorithm makes very mild assumptions about */
81 /* floating point arithmetic. It will work on machines with a guard */
82 /* digit in add/subtract, or on those binary machines without guard */
83 /* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or */
84 /* Cray-2. It could conceivably fail on hexadecimal or decimal machines */
85 /* without guard digits, but we know of none. */
86
87 /* Arguments */
88 /* ========= */
89
90 /* JOBZ (input) CHARACTER*1 */
91 /* = 'N': Compute eigenvalues only; */
92 /* = 'V': Compute eigenvalues and eigenvectors. */
93
94 /* UPLO (input) CHARACTER*1 */
95 /* = 'U': Upper triangle of A is stored; */
96 /* = 'L': Lower triangle of A is stored. */
97
98 /* N (input) INTEGER */
99 /* The order of the matrix A. N >= 0. */
100
101 /* KD (input) INTEGER */
102 /* The number of superdiagonals of the matrix A if UPLO = 'U', */
103 /* or the number of subdiagonals if UPLO = 'L'. KD >= 0. */
104
105 /* AB (input/output) DOUBLE PRECISION array, dimension (LDAB, N) */
106 /* On entry, the upper or lower triangle of the symmetric band */
107 /* matrix A, stored in the first KD+1 rows of the array. The */
108 /* j-th column of A is stored in the j-th column of the array AB */
109 /* as follows: */
110 /* if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; */
111 /* if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). */
112
113 /* On exit, AB is overwritten by values generated during the */
114 /* reduction to tridiagonal form. If UPLO = 'U', the first */
115 /* superdiagonal and the diagonal of the tridiagonal matrix T */
116 /* are returned in rows KD and KD+1 of AB, and if UPLO = 'L', */
117 /* the diagonal and first subdiagonal of T are returned in the */
118 /* first two rows of AB. */
119
120 /* LDAB (input) INTEGER */
121 /* The leading dimension of the array AB. LDAB >= KD + 1. */
122
123 /* W (output) DOUBLE PRECISION array, dimension (N) */
124 /* If INFO = 0, the eigenvalues in ascending order. */
125
126 /* Z (output) DOUBLE PRECISION array, dimension (LDZ, N) */
127 /* If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal */
128 /* eigenvectors of the matrix A, with the i-th column of Z */
129 /* holding the eigenvector associated with W(i). */
130 /* If JOBZ = 'N', then Z is not referenced. */
131
132 /* LDZ (input) INTEGER */
133 /* The leading dimension of the array Z. LDZ >= 1, and if */
134 /* JOBZ = 'V', LDZ >= max(1,N). */
135
136 /* WORK (workspace/output) DOUBLE PRECISION array, */
137 /* dimension (LWORK) */
138 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
139
140 /* LWORK (input) INTEGER */
141 /* The dimension of the array WORK. */
142 /* IF N <= 1, LWORK must be at least 1. */
143 /* If JOBZ = 'N' and N > 2, LWORK must be at least 2*N. */
144 /* If JOBZ = 'V' and N > 2, LWORK must be at least */
145 /* ( 1 + 5*N + 2*N**2 ). */
146
147 /* If LWORK = -1, then a workspace query is assumed; the routine */
148 /* only calculates the optimal size of the WORK array, returns */
149 /* this value as the first entry of the WORK array, and no error */
150 /* message related to LWORK is issued by XERBLA. */
151
152 /* IWORK (workspace/output) INTEGER array, dimension (LIWORK) */
153 /* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. */
154
155 /* LIWORK (input) INTEGER */
156 /* The dimension of the array LIWORK. */
157 /* If JOBZ = 'N' or N <= 1, LIWORK must be at least 1. */
158 /* If JOBZ = 'V' and N > 2, LIWORK must be at least 3 + 5*N. */
159
160 /* If LIWORK = -1, then a workspace query is assumed; the */
161 /* routine only calculates the optimal size of the IWORK array, */
162 /* returns this value as the first entry of the IWORK array, and */
163 /* no error message related to LIWORK is issued by XERBLA. */
164
165 /* INFO (output) INTEGER */
166 /* = 0: successful exit */
167 /* < 0: if INFO = -i, the i-th argument had an illegal value */
168 /* > 0: if INFO = i, the algorithm failed to converge; i */
169 /* off-diagonal elements of an intermediate tridiagonal */
170 /* form did not converge to zero. */
171
172 /* ===================================================================== */
173
174 /* .. Parameters .. */
175 /* .. */
176 /* .. Local Scalars .. */
177 /* .. */
178 /* .. External Functions .. */
179 /* .. */
180 /* .. External Subroutines .. */
181 /* .. */
182 /* .. Intrinsic Functions .. */
183 /* .. */
184 /* .. Executable Statements .. */
185
186 /* Test the input parameters. */
187
188 /* Parameter adjustments */
189 ab_dim1 = *ldab;
190 ab_offset = 1 + ab_dim1;
191 ab -= ab_offset;
192 --w;
193 z_dim1 = *ldz;
194 z_offset = 1 + z_dim1;
195 z__ -= z_offset;
196 --work;
197 --iwork;
198
199 /* Function Body */
200 wantz = lsame_(jobz, "V", (ftnlen)1, (ftnlen)1);
201 lower = lsame_(uplo, "L", (ftnlen)1, (ftnlen)1);
202 lquery = *lwork == -1 || *liwork == -1;
203
204 *info = 0;
205 if (*n <= 1) {
206 liwmin = 1;
207 lwmin = 1;
208 } else {
209 if (wantz) {
210 liwmin = *n * 5 + 3;
211 /* Computing 2nd power */
212 i__1 = *n;
213 lwmin = *n * 5 + 1 + (i__1 * i__1 << 1);
214 } else {
215 liwmin = 1;
216 lwmin = *n << 1;
217 }
218 }
219 if (! (wantz || lsame_(jobz, "N", (ftnlen)1, (ftnlen)1))) {
220 *info = -1;
221 } else if (! (lower || lsame_(uplo, "U", (ftnlen)1, (ftnlen)1))) {
222 *info = -2;
223 } else if (*n < 0) {
224 *info = -3;
225 } else if (*kd < 0) {
226 *info = -4;
227 } else if (*ldab < *kd + 1) {
228 *info = -6;
229 } else if (*ldz < 1 || wantz && *ldz < *n) {
230 *info = -9;
231 } else if (*lwork < lwmin && ! lquery) {
232 *info = -11;
233 } else if (*liwork < liwmin && ! lquery) {
234 *info = -13;
235 }
236
237 if (*info == 0) {
238 work[1] = (doublereal) lwmin;
239 iwork[1] = liwmin;
240 }
241
242 if (*info != 0) {
243 i__1 = -(*info);
244 xerbla_("DSBEVD", &i__1, (ftnlen)6);
245 return 0;
246 } else if (lquery) {
247 return 0;
248 }
249
250 /* Quick return if possible */
251
252 if (*n == 0) {
253 return 0;
254 }
255
256 if (*n == 1) {
257 w[1] = ab[ab_dim1 + 1];
258 if (wantz) {
259 z__[z_dim1 + 1] = 1.;
260 }
261 return 0;
262 }
263
264 /* Get machine constants. */
265
266 safmin = dlamch_("Safe minimum", (ftnlen)12);
267 eps = dlamch_("Precision", (ftnlen)9);
268 smlnum = safmin / eps;
269 bignum = 1. / smlnum;
270 rmin = sqrt(smlnum);
271 rmax = sqrt(bignum);
272
273 /* Scale matrix to allowable range, if necessary. */
274
275 anrm = dlansb_("M", uplo, n, kd, &ab[ab_offset], ldab, &work[1], (ftnlen)
276 1, (ftnlen)1);
277 iscale = 0;
278 if (anrm > 0. && anrm < rmin) {
279 iscale = 1;
280 sigma = rmin / anrm;
281 } else if (anrm > rmax) {
282 iscale = 1;
283 sigma = rmax / anrm;
284 }
285 if (iscale == 1) {
286 if (lower) {
287 dlascl_("B", kd, kd, &c_b11, &sigma, n, n, &ab[ab_offset], ldab,
288 info, (ftnlen)1);
289 } else {
290 dlascl_("Q", kd, kd, &c_b11, &sigma, n, n, &ab[ab_offset], ldab,
291 info, (ftnlen)1);
292 }
293 }
294
295 /* Call DSBTRD to reduce symmetric band matrix to tridiagonal form. */
296
297 inde = 1;
298 indwrk = inde + *n;
299 indwk2 = indwrk + *n * *n;
300 llwrk2 = *lwork - indwk2 + 1;
301 dsbtrd_(jobz, uplo, n, kd, &ab[ab_offset], ldab, &w[1], &work[inde], &z__[
302 z_offset], ldz, &work[indwrk], &iinfo, (ftnlen)1, (ftnlen)1);
303
304 /* For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC. */
305
306 if (! wantz) {
307 dsterf_(n, &w[1], &work[inde], info);
308 } else {
309 dstedc_("I", n, &w[1], &work[inde], &work[indwrk], n, &work[indwk2], &
310 llwrk2, &iwork[1], liwork, info, (ftnlen)1);
311 dgemm_("N", "N", n, n, n, &c_b11, &z__[z_offset], ldz, &work[indwrk],
312 n, &c_b18, &work[indwk2], n, (ftnlen)1, (ftnlen)1);
313 dlacpy_("A", n, n, &work[indwk2], n, &z__[z_offset], ldz, (ftnlen)1);
314 }
315
316 /* If matrix was scaled, then rescale eigenvalues appropriately. */
317
318 if (iscale == 1) {
319 d__1 = 1. / sigma;
320 dscal_(n, &d__1, &w[1], &c__1);
321 }
322
323 work[1] = (doublereal) lwmin;
324 iwork[1] = liwmin;
325 return 0;
326
327 /* End of DSBEVD */
328
329 } /* dsbevd_ */
330
331