1 /* ./src_f77/csytri.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 complex c_b1 = {1.f,0.f};
11 static complex c_b2 = {0.f,0.f};
12 static integer c__1 = 1;
13
csytri_(char * uplo,integer * n,complex * a,integer * lda,integer * ipiv,complex * work,integer * info,ftnlen uplo_len)14 /* Subroutine */ int csytri_(char *uplo, integer *n, complex *a, integer *lda,
15 integer *ipiv, complex *work, integer *info, ftnlen uplo_len)
16 {
17 /* System generated locals */
18 integer a_dim1, a_offset, i__1, i__2, i__3;
19 complex q__1, q__2, q__3;
20
21 /* Builtin functions */
22 void c_div(complex *, complex *, complex *);
23
24 /* Local variables */
25 static complex d__;
26 static integer k;
27 static complex t, ak;
28 static integer kp;
29 static complex akp1, temp, akkp1;
30 extern logical lsame_(char *, char *, ftnlen, ftnlen);
31 extern /* Subroutine */ int ccopy_(integer *, complex *, integer *,
32 complex *, integer *);
33 extern /* Complex */ VOID cdotu_(complex *, integer *, complex *, integer
34 *, complex *, integer *);
35 extern /* Subroutine */ int cswap_(integer *, complex *, integer *,
36 complex *, integer *);
37 static integer kstep;
38 static logical upper;
39 extern /* Subroutine */ int csymv_(char *, integer *, complex *, complex *
40 , integer *, complex *, integer *, complex *, complex *, integer *
41 , ftnlen), xerbla_(char *, integer *, ftnlen);
42
43
44 /* -- LAPACK routine (version 3.0) -- */
45 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
46 /* Courant Institute, Argonne National Lab, and Rice University */
47 /* September 30, 1994 */
48
49 /* .. Scalar Arguments .. */
50 /* .. */
51 /* .. Array Arguments .. */
52 /* .. */
53
54 /* Purpose */
55 /* ======= */
56
57 /* CSYTRI computes the inverse of a complex symmetric indefinite matrix */
58 /* A using the factorization A = U*D*U**T or A = L*D*L**T computed by */
59 /* CSYTRF. */
60
61 /* Arguments */
62 /* ========= */
63
64 /* UPLO (input) CHARACTER*1 */
65 /* Specifies whether the details of the factorization are stored */
66 /* as an upper or lower triangular matrix. */
67 /* = 'U': Upper triangular, form is A = U*D*U**T; */
68 /* = 'L': Lower triangular, form is A = L*D*L**T. */
69
70 /* N (input) INTEGER */
71 /* The order of the matrix A. N >= 0. */
72
73 /* A (input/output) COMPLEX array, dimension (LDA,N) */
74 /* On entry, the block diagonal matrix D and the multipliers */
75 /* used to obtain the factor U or L as computed by CSYTRF. */
76
77 /* On exit, if INFO = 0, the (symmetric) inverse of the original */
78 /* matrix. If UPLO = 'U', the upper triangular part of the */
79 /* inverse is formed and the part of A below the diagonal is not */
80 /* referenced; if UPLO = 'L' the lower triangular part of the */
81 /* inverse is formed and the part of A above the diagonal is */
82 /* not referenced. */
83
84 /* LDA (input) INTEGER */
85 /* The leading dimension of the array A. LDA >= max(1,N). */
86
87 /* IPIV (input) INTEGER array, dimension (N) */
88 /* Details of the interchanges and the block structure of D */
89 /* as determined by CSYTRF. */
90
91 /* WORK (workspace) COMPLEX array, dimension (2*N) */
92
93 /* INFO (output) INTEGER */
94 /* = 0: successful exit */
95 /* < 0: if INFO = -i, the i-th argument had an illegal value */
96 /* > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its */
97 /* inverse could not be computed. */
98
99 /* ===================================================================== */
100
101 /* .. Parameters .. */
102 /* .. */
103 /* .. Local Scalars .. */
104 /* .. */
105 /* .. External Functions .. */
106 /* .. */
107 /* .. External Subroutines .. */
108 /* .. */
109 /* .. Intrinsic Functions .. */
110 /* .. */
111 /* .. Executable Statements .. */
112
113 /* Test the input parameters. */
114
115 /* Parameter adjustments */
116 a_dim1 = *lda;
117 a_offset = 1 + a_dim1;
118 a -= a_offset;
119 --ipiv;
120 --work;
121
122 /* Function Body */
123 *info = 0;
124 upper = lsame_(uplo, "U", (ftnlen)1, (ftnlen)1);
125 if (! upper && ! lsame_(uplo, "L", (ftnlen)1, (ftnlen)1)) {
126 *info = -1;
127 } else if (*n < 0) {
128 *info = -2;
129 } else if (*lda < max(1,*n)) {
130 *info = -4;
131 }
132 if (*info != 0) {
133 i__1 = -(*info);
134 xerbla_("CSYTRI", &i__1, (ftnlen)6);
135 return 0;
136 }
137
138 /* Quick return if possible */
139
140 if (*n == 0) {
141 return 0;
142 }
143
144 /* Check that the diagonal matrix D is nonsingular. */
145
146 if (upper) {
147
148 /* Upper triangular storage: examine D from bottom to top */
149
150 for (*info = *n; *info >= 1; --(*info)) {
151 i__1 = *info + *info * a_dim1;
152 if (ipiv[*info] > 0 && (a[i__1].r == 0.f && a[i__1].i == 0.f)) {
153 return 0;
154 }
155 /* L10: */
156 }
157 } else {
158
159 /* Lower triangular storage: examine D from top to bottom. */
160
161 i__1 = *n;
162 for (*info = 1; *info <= i__1; ++(*info)) {
163 i__2 = *info + *info * a_dim1;
164 if (ipiv[*info] > 0 && (a[i__2].r == 0.f && a[i__2].i == 0.f)) {
165 return 0;
166 }
167 /* L20: */
168 }
169 }
170 *info = 0;
171
172 if (upper) {
173
174 /* Compute inv(A) from the factorization A = U*D*U'. */
175
176 /* K is the main loop index, increasing from 1 to N in steps of */
177 /* 1 or 2, depending on the size of the diagonal blocks. */
178
179 k = 1;
180 L30:
181
182 /* If K > N, exit from loop. */
183
184 if (k > *n) {
185 goto L40;
186 }
187
188 if (ipiv[k] > 0) {
189
190 /* 1 x 1 diagonal block */
191
192 /* Invert the diagonal block. */
193
194 i__1 = k + k * a_dim1;
195 c_div(&q__1, &c_b1, &a[k + k * a_dim1]);
196 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
197
198 /* Compute column K of the inverse. */
199
200 if (k > 1) {
201 i__1 = k - 1;
202 ccopy_(&i__1, &a[k * a_dim1 + 1], &c__1, &work[1], &c__1);
203 i__1 = k - 1;
204 q__1.r = -1.f, q__1.i = -0.f;
205 csymv_(uplo, &i__1, &q__1, &a[a_offset], lda, &work[1], &c__1,
206 &c_b2, &a[k * a_dim1 + 1], &c__1, (ftnlen)1);
207 i__1 = k + k * a_dim1;
208 i__2 = k + k * a_dim1;
209 i__3 = k - 1;
210 cdotu_(&q__2, &i__3, &work[1], &c__1, &a[k * a_dim1 + 1], &
211 c__1);
212 q__1.r = a[i__2].r - q__2.r, q__1.i = a[i__2].i - q__2.i;
213 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
214 }
215 kstep = 1;
216 } else {
217
218 /* 2 x 2 diagonal block */
219
220 /* Invert the diagonal block. */
221
222 i__1 = k + (k + 1) * a_dim1;
223 t.r = a[i__1].r, t.i = a[i__1].i;
224 c_div(&q__1, &a[k + k * a_dim1], &t);
225 ak.r = q__1.r, ak.i = q__1.i;
226 c_div(&q__1, &a[k + 1 + (k + 1) * a_dim1], &t);
227 akp1.r = q__1.r, akp1.i = q__1.i;
228 c_div(&q__1, &a[k + (k + 1) * a_dim1], &t);
229 akkp1.r = q__1.r, akkp1.i = q__1.i;
230 q__3.r = ak.r * akp1.r - ak.i * akp1.i, q__3.i = ak.r * akp1.i +
231 ak.i * akp1.r;
232 q__2.r = q__3.r - 1.f, q__2.i = q__3.i - 0.f;
233 q__1.r = t.r * q__2.r - t.i * q__2.i, q__1.i = t.r * q__2.i + t.i
234 * q__2.r;
235 d__.r = q__1.r, d__.i = q__1.i;
236 i__1 = k + k * a_dim1;
237 c_div(&q__1, &akp1, &d__);
238 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
239 i__1 = k + 1 + (k + 1) * a_dim1;
240 c_div(&q__1, &ak, &d__);
241 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
242 i__1 = k + (k + 1) * a_dim1;
243 q__2.r = -akkp1.r, q__2.i = -akkp1.i;
244 c_div(&q__1, &q__2, &d__);
245 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
246
247 /* Compute columns K and K+1 of the inverse. */
248
249 if (k > 1) {
250 i__1 = k - 1;
251 ccopy_(&i__1, &a[k * a_dim1 + 1], &c__1, &work[1], &c__1);
252 i__1 = k - 1;
253 q__1.r = -1.f, q__1.i = -0.f;
254 csymv_(uplo, &i__1, &q__1, &a[a_offset], lda, &work[1], &c__1,
255 &c_b2, &a[k * a_dim1 + 1], &c__1, (ftnlen)1);
256 i__1 = k + k * a_dim1;
257 i__2 = k + k * a_dim1;
258 i__3 = k - 1;
259 cdotu_(&q__2, &i__3, &work[1], &c__1, &a[k * a_dim1 + 1], &
260 c__1);
261 q__1.r = a[i__2].r - q__2.r, q__1.i = a[i__2].i - q__2.i;
262 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
263 i__1 = k + (k + 1) * a_dim1;
264 i__2 = k + (k + 1) * a_dim1;
265 i__3 = k - 1;
266 cdotu_(&q__2, &i__3, &a[k * a_dim1 + 1], &c__1, &a[(k + 1) *
267 a_dim1 + 1], &c__1);
268 q__1.r = a[i__2].r - q__2.r, q__1.i = a[i__2].i - q__2.i;
269 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
270 i__1 = k - 1;
271 ccopy_(&i__1, &a[(k + 1) * a_dim1 + 1], &c__1, &work[1], &
272 c__1);
273 i__1 = k - 1;
274 q__1.r = -1.f, q__1.i = -0.f;
275 csymv_(uplo, &i__1, &q__1, &a[a_offset], lda, &work[1], &c__1,
276 &c_b2, &a[(k + 1) * a_dim1 + 1], &c__1, (ftnlen)1);
277 i__1 = k + 1 + (k + 1) * a_dim1;
278 i__2 = k + 1 + (k + 1) * a_dim1;
279 i__3 = k - 1;
280 cdotu_(&q__2, &i__3, &work[1], &c__1, &a[(k + 1) * a_dim1 + 1]
281 , &c__1);
282 q__1.r = a[i__2].r - q__2.r, q__1.i = a[i__2].i - q__2.i;
283 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
284 }
285 kstep = 2;
286 }
287
288 kp = (i__1 = ipiv[k], abs(i__1));
289 if (kp != k) {
290
291 /* Interchange rows and columns K and KP in the leading */
292 /* submatrix A(1:k+1,1:k+1) */
293
294 i__1 = kp - 1;
295 cswap_(&i__1, &a[k * a_dim1 + 1], &c__1, &a[kp * a_dim1 + 1], &
296 c__1);
297 i__1 = k - kp - 1;
298 cswap_(&i__1, &a[kp + 1 + k * a_dim1], &c__1, &a[kp + (kp + 1) *
299 a_dim1], lda);
300 i__1 = k + k * a_dim1;
301 temp.r = a[i__1].r, temp.i = a[i__1].i;
302 i__1 = k + k * a_dim1;
303 i__2 = kp + kp * a_dim1;
304 a[i__1].r = a[i__2].r, a[i__1].i = a[i__2].i;
305 i__1 = kp + kp * a_dim1;
306 a[i__1].r = temp.r, a[i__1].i = temp.i;
307 if (kstep == 2) {
308 i__1 = k + (k + 1) * a_dim1;
309 temp.r = a[i__1].r, temp.i = a[i__1].i;
310 i__1 = k + (k + 1) * a_dim1;
311 i__2 = kp + (k + 1) * a_dim1;
312 a[i__1].r = a[i__2].r, a[i__1].i = a[i__2].i;
313 i__1 = kp + (k + 1) * a_dim1;
314 a[i__1].r = temp.r, a[i__1].i = temp.i;
315 }
316 }
317
318 k += kstep;
319 goto L30;
320 L40:
321
322 ;
323 } else {
324
325 /* Compute inv(A) from the factorization A = L*D*L'. */
326
327 /* K is the main loop index, increasing from 1 to N in steps of */
328 /* 1 or 2, depending on the size of the diagonal blocks. */
329
330 k = *n;
331 L50:
332
333 /* If K < 1, exit from loop. */
334
335 if (k < 1) {
336 goto L60;
337 }
338
339 if (ipiv[k] > 0) {
340
341 /* 1 x 1 diagonal block */
342
343 /* Invert the diagonal block. */
344
345 i__1 = k + k * a_dim1;
346 c_div(&q__1, &c_b1, &a[k + k * a_dim1]);
347 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
348
349 /* Compute column K of the inverse. */
350
351 if (k < *n) {
352 i__1 = *n - k;
353 ccopy_(&i__1, &a[k + 1 + k * a_dim1], &c__1, &work[1], &c__1);
354 i__1 = *n - k;
355 q__1.r = -1.f, q__1.i = -0.f;
356 csymv_(uplo, &i__1, &q__1, &a[k + 1 + (k + 1) * a_dim1], lda,
357 &work[1], &c__1, &c_b2, &a[k + 1 + k * a_dim1], &c__1,
358 (ftnlen)1);
359 i__1 = k + k * a_dim1;
360 i__2 = k + k * a_dim1;
361 i__3 = *n - k;
362 cdotu_(&q__2, &i__3, &work[1], &c__1, &a[k + 1 + k * a_dim1],
363 &c__1);
364 q__1.r = a[i__2].r - q__2.r, q__1.i = a[i__2].i - q__2.i;
365 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
366 }
367 kstep = 1;
368 } else {
369
370 /* 2 x 2 diagonal block */
371
372 /* Invert the diagonal block. */
373
374 i__1 = k + (k - 1) * a_dim1;
375 t.r = a[i__1].r, t.i = a[i__1].i;
376 c_div(&q__1, &a[k - 1 + (k - 1) * a_dim1], &t);
377 ak.r = q__1.r, ak.i = q__1.i;
378 c_div(&q__1, &a[k + k * a_dim1], &t);
379 akp1.r = q__1.r, akp1.i = q__1.i;
380 c_div(&q__1, &a[k + (k - 1) * a_dim1], &t);
381 akkp1.r = q__1.r, akkp1.i = q__1.i;
382 q__3.r = ak.r * akp1.r - ak.i * akp1.i, q__3.i = ak.r * akp1.i +
383 ak.i * akp1.r;
384 q__2.r = q__3.r - 1.f, q__2.i = q__3.i - 0.f;
385 q__1.r = t.r * q__2.r - t.i * q__2.i, q__1.i = t.r * q__2.i + t.i
386 * q__2.r;
387 d__.r = q__1.r, d__.i = q__1.i;
388 i__1 = k - 1 + (k - 1) * a_dim1;
389 c_div(&q__1, &akp1, &d__);
390 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
391 i__1 = k + k * a_dim1;
392 c_div(&q__1, &ak, &d__);
393 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
394 i__1 = k + (k - 1) * a_dim1;
395 q__2.r = -akkp1.r, q__2.i = -akkp1.i;
396 c_div(&q__1, &q__2, &d__);
397 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
398
399 /* Compute columns K-1 and K of the inverse. */
400
401 if (k < *n) {
402 i__1 = *n - k;
403 ccopy_(&i__1, &a[k + 1 + k * a_dim1], &c__1, &work[1], &c__1);
404 i__1 = *n - k;
405 q__1.r = -1.f, q__1.i = -0.f;
406 csymv_(uplo, &i__1, &q__1, &a[k + 1 + (k + 1) * a_dim1], lda,
407 &work[1], &c__1, &c_b2, &a[k + 1 + k * a_dim1], &c__1,
408 (ftnlen)1);
409 i__1 = k + k * a_dim1;
410 i__2 = k + k * a_dim1;
411 i__3 = *n - k;
412 cdotu_(&q__2, &i__3, &work[1], &c__1, &a[k + 1 + k * a_dim1],
413 &c__1);
414 q__1.r = a[i__2].r - q__2.r, q__1.i = a[i__2].i - q__2.i;
415 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
416 i__1 = k + (k - 1) * a_dim1;
417 i__2 = k + (k - 1) * a_dim1;
418 i__3 = *n - k;
419 cdotu_(&q__2, &i__3, &a[k + 1 + k * a_dim1], &c__1, &a[k + 1
420 + (k - 1) * a_dim1], &c__1);
421 q__1.r = a[i__2].r - q__2.r, q__1.i = a[i__2].i - q__2.i;
422 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
423 i__1 = *n - k;
424 ccopy_(&i__1, &a[k + 1 + (k - 1) * a_dim1], &c__1, &work[1], &
425 c__1);
426 i__1 = *n - k;
427 q__1.r = -1.f, q__1.i = -0.f;
428 csymv_(uplo, &i__1, &q__1, &a[k + 1 + (k + 1) * a_dim1], lda,
429 &work[1], &c__1, &c_b2, &a[k + 1 + (k - 1) * a_dim1],
430 &c__1, (ftnlen)1);
431 i__1 = k - 1 + (k - 1) * a_dim1;
432 i__2 = k - 1 + (k - 1) * a_dim1;
433 i__3 = *n - k;
434 cdotu_(&q__2, &i__3, &work[1], &c__1, &a[k + 1 + (k - 1) *
435 a_dim1], &c__1);
436 q__1.r = a[i__2].r - q__2.r, q__1.i = a[i__2].i - q__2.i;
437 a[i__1].r = q__1.r, a[i__1].i = q__1.i;
438 }
439 kstep = 2;
440 }
441
442 kp = (i__1 = ipiv[k], abs(i__1));
443 if (kp != k) {
444
445 /* Interchange rows and columns K and KP in the trailing */
446 /* submatrix A(k-1:n,k-1:n) */
447
448 if (kp < *n) {
449 i__1 = *n - kp;
450 cswap_(&i__1, &a[kp + 1 + k * a_dim1], &c__1, &a[kp + 1 + kp *
451 a_dim1], &c__1);
452 }
453 i__1 = kp - k - 1;
454 cswap_(&i__1, &a[k + 1 + k * a_dim1], &c__1, &a[kp + (k + 1) *
455 a_dim1], lda);
456 i__1 = k + k * a_dim1;
457 temp.r = a[i__1].r, temp.i = a[i__1].i;
458 i__1 = k + k * a_dim1;
459 i__2 = kp + kp * a_dim1;
460 a[i__1].r = a[i__2].r, a[i__1].i = a[i__2].i;
461 i__1 = kp + kp * a_dim1;
462 a[i__1].r = temp.r, a[i__1].i = temp.i;
463 if (kstep == 2) {
464 i__1 = k + (k - 1) * a_dim1;
465 temp.r = a[i__1].r, temp.i = a[i__1].i;
466 i__1 = k + (k - 1) * a_dim1;
467 i__2 = kp + (k - 1) * a_dim1;
468 a[i__1].r = a[i__2].r, a[i__1].i = a[i__2].i;
469 i__1 = kp + (k - 1) * a_dim1;
470 a[i__1].r = temp.r, a[i__1].i = temp.i;
471 }
472 }
473
474 k -= kstep;
475 goto L50;
476 L60:
477 ;
478 }
479
480 return 0;
481
482 /* End of CSYTRI */
483
484 } /* csytri_ */
485
486