1 /* ./src_f77/dstevd.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
dstevd_(char * jobz,integer * n,doublereal * d__,doublereal * e,doublereal * z__,integer * ldz,doublereal * work,integer * lwork,integer * iwork,integer * liwork,integer * info,ftnlen jobz_len)12 /* Subroutine */ int dstevd_(char *jobz, integer *n, doublereal *d__,
13 doublereal *e, doublereal *z__, integer *ldz, doublereal *work,
14 integer *lwork, integer *iwork, integer *liwork, integer *info,
15 ftnlen jobz_len)
16 {
17 /* System generated locals */
18 integer z_dim1, z_offset, i__1;
19 doublereal d__1;
20
21 /* Builtin functions */
22 double sqrt(doublereal);
23
24 /* Local variables */
25 static doublereal eps, rmin, rmax, tnrm;
26 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
27 integer *);
28 static doublereal sigma;
29 extern logical lsame_(char *, char *, ftnlen, ftnlen);
30 static integer lwmin;
31 static logical wantz;
32 extern doublereal dlamch_(char *, ftnlen);
33 static integer iscale;
34 extern /* Subroutine */ int dstedc_(char *, integer *, doublereal *,
35 doublereal *, doublereal *, integer *, doublereal *, integer *,
36 integer *, integer *, integer *, ftnlen);
37 static doublereal safmin;
38 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
39 static doublereal bignum;
40 extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *,
41 ftnlen);
42 extern /* Subroutine */ int dsterf_(integer *, doublereal *, doublereal *,
43 integer *);
44 static integer liwmin;
45 static doublereal smlnum;
46 static logical lquery;
47
48
49 /* -- LAPACK driver routine (version 3.0) -- */
50 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
51 /* Courant Institute, Argonne National Lab, and Rice University */
52 /* June 30, 1999 */
53
54 /* .. Scalar Arguments .. */
55 /* .. */
56 /* .. Array Arguments .. */
57 /* .. */
58
59 /* Purpose */
60 /* ======= */
61
62 /* DSTEVD computes all eigenvalues and, optionally, eigenvectors of a */
63 /* real symmetric tridiagonal matrix. If eigenvectors are desired, it */
64 /* uses a divide and conquer algorithm. */
65
66 /* The divide and conquer algorithm makes very mild assumptions about */
67 /* floating point arithmetic. It will work on machines with a guard */
68 /* digit in add/subtract, or on those binary machines without guard */
69 /* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or */
70 /* Cray-2. It could conceivably fail on hexadecimal or decimal machines */
71 /* without guard digits, but we know of none. */
72
73 /* Arguments */
74 /* ========= */
75
76 /* JOBZ (input) CHARACTER*1 */
77 /* = 'N': Compute eigenvalues only; */
78 /* = 'V': Compute eigenvalues and eigenvectors. */
79
80 /* N (input) INTEGER */
81 /* The order of the matrix. N >= 0. */
82
83 /* D (input/output) DOUBLE PRECISION array, dimension (N) */
84 /* On entry, the n diagonal elements of the tridiagonal matrix */
85 /* A. */
86 /* On exit, if INFO = 0, the eigenvalues in ascending order. */
87
88 /* E (input/output) DOUBLE PRECISION array, dimension (N) */
89 /* On entry, the (n-1) subdiagonal elements of the tridiagonal */
90 /* matrix A, stored in elements 1 to N-1 of E; E(N) need not */
91 /* be set, but is used by the routine. */
92 /* On exit, the contents of E are destroyed. */
93
94 /* Z (output) DOUBLE PRECISION array, dimension (LDZ, N) */
95 /* If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal */
96 /* eigenvectors of the matrix A, with the i-th column of Z */
97 /* holding the eigenvector associated with D(i). */
98 /* If JOBZ = 'N', then Z is not referenced. */
99
100 /* LDZ (input) INTEGER */
101 /* The leading dimension of the array Z. LDZ >= 1, and if */
102 /* JOBZ = 'V', LDZ >= max(1,N). */
103
104 /* WORK (workspace/output) DOUBLE PRECISION array, */
105 /* dimension (LWORK) */
106 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
107
108 /* LWORK (input) INTEGER */
109 /* The dimension of the array WORK. */
110 /* If JOBZ = 'N' or N <= 1 then LWORK must be at least 1. */
111 /* If JOBZ = 'V' and N > 1 then LWORK must be at least */
112 /* ( 1 + 4*N + N**2 ). */
113
114 /* If LWORK = -1, then a workspace query is assumed; the routine */
115 /* only calculates the optimal size of the WORK array, returns */
116 /* this value as the first entry of the WORK array, and no error */
117 /* message related to LWORK is issued by XERBLA. */
118
119 /* IWORK (workspace/output) INTEGER array, dimension (LIWORK) */
120 /* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. */
121
122 /* LIWORK (input) INTEGER */
123 /* The dimension of the array IWORK. */
124 /* If JOBZ = 'N' or N <= 1 then LIWORK must be at least 1. */
125 /* If JOBZ = 'V' and N > 1 then LIWORK must be at least 3+5*N. */
126
127 /* If LIWORK = -1, then a workspace query is assumed; the */
128 /* routine only calculates the optimal size of the IWORK array, */
129 /* returns this value as the first entry of the IWORK array, and */
130 /* no error message related to LIWORK is issued by XERBLA. */
131
132 /* INFO (output) INTEGER */
133 /* = 0: successful exit */
134 /* < 0: if INFO = -i, the i-th argument had an illegal value */
135 /* > 0: if INFO = i, the algorithm failed to converge; i */
136 /* off-diagonal elements of E did not converge to zero. */
137
138 /* ===================================================================== */
139
140 /* .. Parameters .. */
141 /* .. */
142 /* .. Local Scalars .. */
143 /* .. */
144 /* .. External Functions .. */
145 /* .. */
146 /* .. External Subroutines .. */
147 /* .. */
148 /* .. Intrinsic Functions .. */
149 /* .. */
150 /* .. Executable Statements .. */
151
152 /* Test the input parameters. */
153
154 /* Parameter adjustments */
155 --d__;
156 --e;
157 z_dim1 = *ldz;
158 z_offset = 1 + z_dim1;
159 z__ -= z_offset;
160 --work;
161 --iwork;
162
163 /* Function Body */
164 wantz = lsame_(jobz, "V", (ftnlen)1, (ftnlen)1);
165 lquery = *lwork == -1 || *liwork == -1;
166
167 *info = 0;
168 liwmin = 1;
169 lwmin = 1;
170 if (*n > 1 && wantz) {
171 /* Computing 2nd power */
172 i__1 = *n;
173 lwmin = (*n << 2) + 1 + i__1 * i__1;
174 liwmin = *n * 5 + 3;
175 }
176
177 if (! (wantz || lsame_(jobz, "N", (ftnlen)1, (ftnlen)1))) {
178 *info = -1;
179 } else if (*n < 0) {
180 *info = -2;
181 } else if (*ldz < 1 || wantz && *ldz < *n) {
182 *info = -6;
183 } else if (*lwork < lwmin && ! lquery) {
184 *info = -8;
185 } else if (*liwork < liwmin && ! lquery) {
186 *info = -10;
187 }
188
189 if (*info == 0) {
190 work[1] = (doublereal) lwmin;
191 iwork[1] = liwmin;
192 }
193
194 if (*info != 0) {
195 i__1 = -(*info);
196 xerbla_("DSTEVD", &i__1, (ftnlen)6);
197 return 0;
198 } else if (lquery) {
199 return 0;
200 }
201
202 /* Quick return if possible */
203
204 if (*n == 0) {
205 return 0;
206 }
207
208 if (*n == 1) {
209 if (wantz) {
210 z__[z_dim1 + 1] = 1.;
211 }
212 return 0;
213 }
214
215 /* Get machine constants. */
216
217 safmin = dlamch_("Safe minimum", (ftnlen)12);
218 eps = dlamch_("Precision", (ftnlen)9);
219 smlnum = safmin / eps;
220 bignum = 1. / smlnum;
221 rmin = sqrt(smlnum);
222 rmax = sqrt(bignum);
223
224 /* Scale matrix to allowable range, if necessary. */
225
226 iscale = 0;
227 tnrm = dlanst_("M", n, &d__[1], &e[1], (ftnlen)1);
228 if (tnrm > 0. && tnrm < rmin) {
229 iscale = 1;
230 sigma = rmin / tnrm;
231 } else if (tnrm > rmax) {
232 iscale = 1;
233 sigma = rmax / tnrm;
234 }
235 if (iscale == 1) {
236 dscal_(n, &sigma, &d__[1], &c__1);
237 i__1 = *n - 1;
238 dscal_(&i__1, &sigma, &e[1], &c__1);
239 }
240
241 /* For eigenvalues only, call DSTERF. For eigenvalues and */
242 /* eigenvectors, call DSTEDC. */
243
244 if (! wantz) {
245 dsterf_(n, &d__[1], &e[1], info);
246 } else {
247 dstedc_("I", n, &d__[1], &e[1], &z__[z_offset], ldz, &work[1], lwork,
248 &iwork[1], liwork, info, (ftnlen)1);
249 }
250
251 /* If matrix was scaled, then rescale eigenvalues appropriately. */
252
253 if (iscale == 1) {
254 d__1 = 1. / sigma;
255 dscal_(n, &d__1, &d__[1], &c__1);
256 }
257
258 work[1] = (doublereal) lwmin;
259 iwork[1] = liwmin;
260
261 return 0;
262
263 /* End of DSTEVD */
264
265 } /* dstevd_ */
266
267