1 /* ./src_f77/zhegs2.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 = {1.,0.};
11 static integer c__1 = 1;
12
zhegs2_(integer * itype,char * uplo,integer * n,doublecomplex * a,integer * lda,doublecomplex * b,integer * ldb,integer * info,ftnlen uplo_len)13 /* Subroutine */ int zhegs2_(integer *itype, char *uplo, integer *n,
14 doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
15 integer *info, ftnlen uplo_len)
16 {
17 /* System generated locals */
18 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;
19 doublereal d__1, d__2;
20 doublecomplex z__1;
21
22 /* Local variables */
23 static integer k;
24 static doublecomplex ct;
25 static doublereal akk, bkk;
26 extern /* Subroutine */ int zher2_(char *, integer *, doublecomplex *,
27 doublecomplex *, integer *, doublecomplex *, integer *,
28 doublecomplex *, integer *, ftnlen);
29 extern logical lsame_(char *, char *, ftnlen, ftnlen);
30 static logical upper;
31 extern /* Subroutine */ int zaxpy_(integer *, doublecomplex *,
32 doublecomplex *, integer *, doublecomplex *, integer *), ztrmv_(
33 char *, char *, char *, integer *, doublecomplex *, integer *,
34 doublecomplex *, integer *, ftnlen, ftnlen, ftnlen), ztrsv_(char *
35 , char *, char *, integer *, doublecomplex *, integer *,
36 doublecomplex *, integer *, ftnlen, ftnlen, ftnlen), xerbla_(char
37 *, integer *, ftnlen), zdscal_(integer *, doublereal *,
38 doublecomplex *, integer *), zlacgv_(integer *, doublecomplex *,
39 integer *);
40
41
42 /* -- LAPACK routine (version 3.0) -- */
43 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
44 /* Courant Institute, Argonne National Lab, and Rice University */
45 /* September 30, 1994 */
46
47 /* .. Scalar Arguments .. */
48 /* .. */
49 /* .. Array Arguments .. */
50 /* .. */
51
52 /* Purpose */
53 /* ======= */
54
55 /* ZHEGS2 reduces a complex Hermitian-definite generalized */
56 /* eigenproblem to standard form. */
57
58 /* If ITYPE = 1, the problem is A*x = lambda*B*x, */
59 /* and A is overwritten by inv(U')*A*inv(U) or inv(L)*A*inv(L') */
60
61 /* If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or */
62 /* B*A*x = lambda*x, and A is overwritten by U*A*U` or L'*A*L. */
63
64 /* B must have been previously factorized as U'*U or L*L' by ZPOTRF. */
65
66 /* Arguments */
67 /* ========= */
68
69 /* ITYPE (input) INTEGER */
70 /* = 1: compute inv(U')*A*inv(U) or inv(L)*A*inv(L'); */
71 /* = 2 or 3: compute U*A*U' or L'*A*L. */
72
73 /* UPLO (input) CHARACTER */
74 /* Specifies whether the upper or lower triangular part of the */
75 /* Hermitian matrix A is stored, and how B has been factorized. */
76 /* = 'U': Upper triangular */
77 /* = 'L': Lower triangular */
78
79 /* N (input) INTEGER */
80 /* The order of the matrices A and B. N >= 0. */
81
82 /* A (input/output) COMPLEX*16 array, dimension (LDA,N) */
83 /* On entry, the Hermitian matrix A. If UPLO = 'U', the leading */
84 /* n by n upper triangular part of A contains the upper */
85 /* triangular part of the matrix A, and the strictly lower */
86 /* triangular part of A is not referenced. If UPLO = 'L', the */
87 /* leading n by n lower triangular part of A contains the lower */
88 /* triangular part of the matrix A, and the strictly upper */
89 /* triangular part of A is not referenced. */
90
91 /* On exit, if INFO = 0, the transformed matrix, stored in the */
92 /* same format as A. */
93
94 /* LDA (input) INTEGER */
95 /* The leading dimension of the array A. LDA >= max(1,N). */
96
97 /* B (input) COMPLEX*16 array, dimension (LDB,N) */
98 /* The triangular factor from the Cholesky factorization of B, */
99 /* as returned by ZPOTRF. */
100
101 /* LDB (input) INTEGER */
102 /* The leading dimension of the array B. LDB >= max(1,N). */
103
104 /* INFO (output) INTEGER */
105 /* = 0: successful exit. */
106 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
107
108 /* ===================================================================== */
109
110 /* .. Parameters .. */
111 /* .. */
112 /* .. Local Scalars .. */
113 /* .. */
114 /* .. External Subroutines .. */
115 /* .. */
116 /* .. Intrinsic Functions .. */
117 /* .. */
118 /* .. External Functions .. */
119 /* .. */
120 /* .. Executable Statements .. */
121
122 /* Test the input parameters. */
123
124 /* Parameter adjustments */
125 a_dim1 = *lda;
126 a_offset = 1 + a_dim1;
127 a -= a_offset;
128 b_dim1 = *ldb;
129 b_offset = 1 + b_dim1;
130 b -= b_offset;
131
132 /* Function Body */
133 *info = 0;
134 upper = lsame_(uplo, "U", (ftnlen)1, (ftnlen)1);
135 if (*itype < 1 || *itype > 3) {
136 *info = -1;
137 } else if (! upper && ! lsame_(uplo, "L", (ftnlen)1, (ftnlen)1)) {
138 *info = -2;
139 } else if (*n < 0) {
140 *info = -3;
141 } else if (*lda < max(1,*n)) {
142 *info = -5;
143 } else if (*ldb < max(1,*n)) {
144 *info = -7;
145 }
146 if (*info != 0) {
147 i__1 = -(*info);
148 xerbla_("ZHEGS2", &i__1, (ftnlen)6);
149 return 0;
150 }
151
152 if (*itype == 1) {
153 if (upper) {
154
155 /* Compute inv(U')*A*inv(U) */
156
157 i__1 = *n;
158 for (k = 1; k <= i__1; ++k) {
159
160 /* Update the upper triangle of A(k:n,k:n) */
161
162 i__2 = k + k * a_dim1;
163 akk = a[i__2].r;
164 i__2 = k + k * b_dim1;
165 bkk = b[i__2].r;
166 /* Computing 2nd power */
167 d__1 = bkk;
168 akk /= d__1 * d__1;
169 i__2 = k + k * a_dim1;
170 a[i__2].r = akk, a[i__2].i = 0.;
171 if (k < *n) {
172 i__2 = *n - k;
173 d__1 = 1. / bkk;
174 zdscal_(&i__2, &d__1, &a[k + (k + 1) * a_dim1], lda);
175 d__1 = akk * -.5;
176 ct.r = d__1, ct.i = 0.;
177 i__2 = *n - k;
178 zlacgv_(&i__2, &a[k + (k + 1) * a_dim1], lda);
179 i__2 = *n - k;
180 zlacgv_(&i__2, &b[k + (k + 1) * b_dim1], ldb);
181 i__2 = *n - k;
182 zaxpy_(&i__2, &ct, &b[k + (k + 1) * b_dim1], ldb, &a[k + (
183 k + 1) * a_dim1], lda);
184 i__2 = *n - k;
185 z__1.r = -1., z__1.i = -0.;
186 zher2_(uplo, &i__2, &z__1, &a[k + (k + 1) * a_dim1], lda,
187 &b[k + (k + 1) * b_dim1], ldb, &a[k + 1 + (k + 1)
188 * a_dim1], lda, (ftnlen)1);
189 i__2 = *n - k;
190 zaxpy_(&i__2, &ct, &b[k + (k + 1) * b_dim1], ldb, &a[k + (
191 k + 1) * a_dim1], lda);
192 i__2 = *n - k;
193 zlacgv_(&i__2, &b[k + (k + 1) * b_dim1], ldb);
194 i__2 = *n - k;
195 ztrsv_(uplo, "Conjugate transpose", "Non-unit", &i__2, &b[
196 k + 1 + (k + 1) * b_dim1], ldb, &a[k + (k + 1) *
197 a_dim1], lda, (ftnlen)1, (ftnlen)19, (ftnlen)8);
198 i__2 = *n - k;
199 zlacgv_(&i__2, &a[k + (k + 1) * a_dim1], lda);
200 }
201 /* L10: */
202 }
203 } else {
204
205 /* Compute inv(L)*A*inv(L') */
206
207 i__1 = *n;
208 for (k = 1; k <= i__1; ++k) {
209
210 /* Update the lower triangle of A(k:n,k:n) */
211
212 i__2 = k + k * a_dim1;
213 akk = a[i__2].r;
214 i__2 = k + k * b_dim1;
215 bkk = b[i__2].r;
216 /* Computing 2nd power */
217 d__1 = bkk;
218 akk /= d__1 * d__1;
219 i__2 = k + k * a_dim1;
220 a[i__2].r = akk, a[i__2].i = 0.;
221 if (k < *n) {
222 i__2 = *n - k;
223 d__1 = 1. / bkk;
224 zdscal_(&i__2, &d__1, &a[k + 1 + k * a_dim1], &c__1);
225 d__1 = akk * -.5;
226 ct.r = d__1, ct.i = 0.;
227 i__2 = *n - k;
228 zaxpy_(&i__2, &ct, &b[k + 1 + k * b_dim1], &c__1, &a[k +
229 1 + k * a_dim1], &c__1);
230 i__2 = *n - k;
231 z__1.r = -1., z__1.i = -0.;
232 zher2_(uplo, &i__2, &z__1, &a[k + 1 + k * a_dim1], &c__1,
233 &b[k + 1 + k * b_dim1], &c__1, &a[k + 1 + (k + 1)
234 * a_dim1], lda, (ftnlen)1);
235 i__2 = *n - k;
236 zaxpy_(&i__2, &ct, &b[k + 1 + k * b_dim1], &c__1, &a[k +
237 1 + k * a_dim1], &c__1);
238 i__2 = *n - k;
239 ztrsv_(uplo, "No transpose", "Non-unit", &i__2, &b[k + 1
240 + (k + 1) * b_dim1], ldb, &a[k + 1 + k * a_dim1],
241 &c__1, (ftnlen)1, (ftnlen)12, (ftnlen)8);
242 }
243 /* L20: */
244 }
245 }
246 } else {
247 if (upper) {
248
249 /* Compute U*A*U' */
250
251 i__1 = *n;
252 for (k = 1; k <= i__1; ++k) {
253
254 /* Update the upper triangle of A(1:k,1:k) */
255
256 i__2 = k + k * a_dim1;
257 akk = a[i__2].r;
258 i__2 = k + k * b_dim1;
259 bkk = b[i__2].r;
260 i__2 = k - 1;
261 ztrmv_(uplo, "No transpose", "Non-unit", &i__2, &b[b_offset],
262 ldb, &a[k * a_dim1 + 1], &c__1, (ftnlen)1, (ftnlen)12,
263 (ftnlen)8);
264 d__1 = akk * .5;
265 ct.r = d__1, ct.i = 0.;
266 i__2 = k - 1;
267 zaxpy_(&i__2, &ct, &b[k * b_dim1 + 1], &c__1, &a[k * a_dim1 +
268 1], &c__1);
269 i__2 = k - 1;
270 zher2_(uplo, &i__2, &c_b1, &a[k * a_dim1 + 1], &c__1, &b[k *
271 b_dim1 + 1], &c__1, &a[a_offset], lda, (ftnlen)1);
272 i__2 = k - 1;
273 zaxpy_(&i__2, &ct, &b[k * b_dim1 + 1], &c__1, &a[k * a_dim1 +
274 1], &c__1);
275 i__2 = k - 1;
276 zdscal_(&i__2, &bkk, &a[k * a_dim1 + 1], &c__1);
277 i__2 = k + k * a_dim1;
278 /* Computing 2nd power */
279 d__2 = bkk;
280 d__1 = akk * (d__2 * d__2);
281 a[i__2].r = d__1, a[i__2].i = 0.;
282 /* L30: */
283 }
284 } else {
285
286 /* Compute L'*A*L */
287
288 i__1 = *n;
289 for (k = 1; k <= i__1; ++k) {
290
291 /* Update the lower triangle of A(1:k,1:k) */
292
293 i__2 = k + k * a_dim1;
294 akk = a[i__2].r;
295 i__2 = k + k * b_dim1;
296 bkk = b[i__2].r;
297 i__2 = k - 1;
298 zlacgv_(&i__2, &a[k + a_dim1], lda);
299 i__2 = k - 1;
300 ztrmv_(uplo, "Conjugate transpose", "Non-unit", &i__2, &b[
301 b_offset], ldb, &a[k + a_dim1], lda, (ftnlen)1, (
302 ftnlen)19, (ftnlen)8);
303 d__1 = akk * .5;
304 ct.r = d__1, ct.i = 0.;
305 i__2 = k - 1;
306 zlacgv_(&i__2, &b[k + b_dim1], ldb);
307 i__2 = k - 1;
308 zaxpy_(&i__2, &ct, &b[k + b_dim1], ldb, &a[k + a_dim1], lda);
309 i__2 = k - 1;
310 zher2_(uplo, &i__2, &c_b1, &a[k + a_dim1], lda, &b[k + b_dim1]
311 , ldb, &a[a_offset], lda, (ftnlen)1);
312 i__2 = k - 1;
313 zaxpy_(&i__2, &ct, &b[k + b_dim1], ldb, &a[k + a_dim1], lda);
314 i__2 = k - 1;
315 zlacgv_(&i__2, &b[k + b_dim1], ldb);
316 i__2 = k - 1;
317 zdscal_(&i__2, &bkk, &a[k + a_dim1], lda);
318 i__2 = k - 1;
319 zlacgv_(&i__2, &a[k + a_dim1], lda);
320 i__2 = k + k * a_dim1;
321 /* Computing 2nd power */
322 d__2 = bkk;
323 d__1 = akk * (d__2 * d__2);
324 a[i__2].r = d__1, a[i__2].i = 0.;
325 /* L40: */
326 }
327 }
328 }
329 return 0;
330
331 /* End of ZHEGS2 */
332
333 } /* zhegs2_ */
334
335