1 /* ./src_f77/dspev.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
dspev_(char * jobz,char * uplo,integer * n,doublereal * ap,doublereal * w,doublereal * z__,integer * ldz,doublereal * work,integer * info,ftnlen jobz_len,ftnlen uplo_len)12 /* Subroutine */ int dspev_(char *jobz, char *uplo, integer *n, doublereal *
13 ap, doublereal *w, doublereal *z__, integer *ldz, doublereal *work,
14 integer *info, ftnlen jobz_len, ftnlen uplo_len)
15 {
16 /* System generated locals */
17 integer z_dim1, z_offset, i__1;
18 doublereal d__1;
19
20 /* Builtin functions */
21 double sqrt(doublereal);
22
23 /* Local variables */
24 static doublereal eps;
25 static integer inde;
26 static doublereal anrm;
27 static integer imax;
28 static doublereal rmin, rmax;
29 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
30 integer *);
31 static doublereal sigma;
32 extern logical lsame_(char *, char *, ftnlen, ftnlen);
33 static integer iinfo;
34 static logical wantz;
35 extern doublereal dlamch_(char *, ftnlen);
36 static integer iscale;
37 static doublereal safmin;
38 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
39 static doublereal bignum;
40 extern doublereal dlansp_(char *, char *, integer *, doublereal *,
41 doublereal *, ftnlen, ftnlen);
42 static integer indtau;
43 extern /* Subroutine */ int dsterf_(integer *, doublereal *, doublereal *,
44 integer *);
45 static integer indwrk;
46 extern /* Subroutine */ int dopgtr_(char *, integer *, doublereal *,
47 doublereal *, doublereal *, integer *, doublereal *, integer *,
48 ftnlen), dsptrd_(char *, integer *, doublereal *, doublereal *,
49 doublereal *, doublereal *, integer *, ftnlen), dsteqr_(char *,
50 integer *, doublereal *, doublereal *, doublereal *, integer *,
51 doublereal *, integer *, ftnlen);
52 static doublereal smlnum;
53
54
55 /* -- LAPACK driver routine (version 3.0) -- */
56 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
57 /* Courant Institute, Argonne National Lab, and Rice University */
58 /* March 31, 1993 */
59
60 /* .. Scalar Arguments .. */
61 /* .. */
62 /* .. Array Arguments .. */
63 /* .. */
64
65 /* Purpose */
66 /* ======= */
67
68 /* DSPEV computes all the eigenvalues and, optionally, eigenvectors of a */
69 /* real symmetric matrix A in packed storage. */
70
71 /* Arguments */
72 /* ========= */
73
74 /* JOBZ (input) CHARACTER*1 */
75 /* = 'N': Compute eigenvalues only; */
76 /* = 'V': Compute eigenvalues and eigenvectors. */
77
78 /* UPLO (input) CHARACTER*1 */
79 /* = 'U': Upper triangle of A is stored; */
80 /* = 'L': Lower triangle of A is stored. */
81
82 /* N (input) INTEGER */
83 /* The order of the matrix A. N >= 0. */
84
85 /* AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2) */
86 /* On entry, the upper or lower triangle of the symmetric matrix */
87 /* A, packed columnwise in a linear array. The j-th column of A */
88 /* is stored in the array AP as follows: */
89 /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
90 /* if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. */
91
92 /* On exit, AP is overwritten by values generated during the */
93 /* reduction to tridiagonal form. If UPLO = 'U', the diagonal */
94 /* and first superdiagonal of the tridiagonal matrix T overwrite */
95 /* the corresponding elements of A, and if UPLO = 'L', the */
96 /* diagonal and first subdiagonal of T overwrite the */
97 /* corresponding elements of A. */
98
99 /* W (output) DOUBLE PRECISION array, dimension (N) */
100 /* If INFO = 0, the eigenvalues in ascending order. */
101
102 /* Z (output) DOUBLE PRECISION array, dimension (LDZ, N) */
103 /* If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal */
104 /* eigenvectors of the matrix A, with the i-th column of Z */
105 /* holding the eigenvector associated with W(i). */
106 /* If JOBZ = 'N', then Z is not referenced. */
107
108 /* LDZ (input) INTEGER */
109 /* The leading dimension of the array Z. LDZ >= 1, and if */
110 /* JOBZ = 'V', LDZ >= max(1,N). */
111
112 /* WORK (workspace) DOUBLE PRECISION array, dimension (3*N) */
113
114 /* INFO (output) INTEGER */
115 /* = 0: successful exit. */
116 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
117 /* > 0: if INFO = i, the algorithm failed to converge; i */
118 /* off-diagonal elements of an intermediate tridiagonal */
119 /* form did not converge to zero. */
120
121 /* ===================================================================== */
122
123 /* .. Parameters .. */
124 /* .. */
125 /* .. Local Scalars .. */
126 /* .. */
127 /* .. External Functions .. */
128 /* .. */
129 /* .. External Subroutines .. */
130 /* .. */
131 /* .. Intrinsic Functions .. */
132 /* .. */
133 /* .. Executable Statements .. */
134
135 /* Test the input parameters. */
136
137 /* Parameter adjustments */
138 --ap;
139 --w;
140 z_dim1 = *ldz;
141 z_offset = 1 + z_dim1;
142 z__ -= z_offset;
143 --work;
144
145 /* Function Body */
146 wantz = lsame_(jobz, "V", (ftnlen)1, (ftnlen)1);
147
148 *info = 0;
149 if (! (wantz || lsame_(jobz, "N", (ftnlen)1, (ftnlen)1))) {
150 *info = -1;
151 } else if (! (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) || lsame_(uplo,
152 "L", (ftnlen)1, (ftnlen)1))) {
153 *info = -2;
154 } else if (*n < 0) {
155 *info = -3;
156 } else if (*ldz < 1 || wantz && *ldz < *n) {
157 *info = -7;
158 }
159
160 if (*info != 0) {
161 i__1 = -(*info);
162 xerbla_("DSPEV ", &i__1, (ftnlen)6);
163 return 0;
164 }
165
166 /* Quick return if possible */
167
168 if (*n == 0) {
169 return 0;
170 }
171
172 if (*n == 1) {
173 w[1] = ap[1];
174 if (wantz) {
175 z__[z_dim1 + 1] = 1.;
176 }
177 return 0;
178 }
179
180 /* Get machine constants. */
181
182 safmin = dlamch_("Safe minimum", (ftnlen)12);
183 eps = dlamch_("Precision", (ftnlen)9);
184 smlnum = safmin / eps;
185 bignum = 1. / smlnum;
186 rmin = sqrt(smlnum);
187 rmax = sqrt(bignum);
188
189 /* Scale matrix to allowable range, if necessary. */
190
191 anrm = dlansp_("M", uplo, n, &ap[1], &work[1], (ftnlen)1, (ftnlen)1);
192 iscale = 0;
193 if (anrm > 0. && anrm < rmin) {
194 iscale = 1;
195 sigma = rmin / anrm;
196 } else if (anrm > rmax) {
197 iscale = 1;
198 sigma = rmax / anrm;
199 }
200 if (iscale == 1) {
201 i__1 = *n * (*n + 1) / 2;
202 dscal_(&i__1, &sigma, &ap[1], &c__1);
203 }
204
205 /* Call DSPTRD to reduce symmetric packed matrix to tridiagonal form. */
206
207 inde = 1;
208 indtau = inde + *n;
209 dsptrd_(uplo, n, &ap[1], &w[1], &work[inde], &work[indtau], &iinfo, (
210 ftnlen)1);
211
212 /* For eigenvalues only, call DSTERF. For eigenvectors, first call */
213 /* DOPGTR to generate the orthogonal matrix, then call DSTEQR. */
214
215 if (! wantz) {
216 dsterf_(n, &w[1], &work[inde], info);
217 } else {
218 indwrk = indtau + *n;
219 dopgtr_(uplo, n, &ap[1], &work[indtau], &z__[z_offset], ldz, &work[
220 indwrk], &iinfo, (ftnlen)1);
221 dsteqr_(jobz, n, &w[1], &work[inde], &z__[z_offset], ldz, &work[
222 indtau], info, (ftnlen)1);
223 }
224
225 /* If matrix was scaled, then rescale eigenvalues appropriately. */
226
227 if (iscale == 1) {
228 if (*info == 0) {
229 imax = *n;
230 } else {
231 imax = *info - 1;
232 }
233 d__1 = 1. / sigma;
234 dscal_(&imax, &d__1, &w[1], &c__1);
235 }
236
237 return 0;
238
239 /* End of DSPEV */
240
241 } /* dspev_ */
242
243