1 /* ./src_f77/cpbtrf.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 integer c__1 = 1;
12 static integer c_n1 = -1;
13 static real c_b21 = -1.f;
14 static real c_b22 = 1.f;
15 static integer c__33 = 33;
16
cpbtrf_(char * uplo,integer * n,integer * kd,complex * ab,integer * ldab,integer * info,ftnlen uplo_len)17 /* Subroutine */ int cpbtrf_(char *uplo, integer *n, integer *kd, complex *ab,
18 integer *ldab, integer *info, ftnlen uplo_len)
19 {
20 /* System generated locals */
21 integer ab_dim1, ab_offset, i__1, i__2, i__3, i__4, i__5, i__6;
22 complex q__1;
23
24 /* Local variables */
25 static integer i__, j, i2, i3, ib, nb, ii, jj;
26 static complex work[1056] /* was [33][32] */;
27 extern /* Subroutine */ int cgemm_(char *, char *, integer *, integer *,
28 integer *, complex *, complex *, integer *, complex *, integer *,
29 complex *, complex *, integer *, ftnlen, ftnlen), cherk_(char *,
30 char *, integer *, integer *, real *, complex *, integer *, real *
31 , complex *, integer *, ftnlen, ftnlen);
32 extern logical lsame_(char *, char *, ftnlen, ftnlen);
33 extern /* Subroutine */ int ctrsm_(char *, char *, char *, char *,
34 integer *, integer *, complex *, complex *, integer *, complex *,
35 integer *, ftnlen, ftnlen, ftnlen, ftnlen), cpbtf2_(char *,
36 integer *, integer *, complex *, integer *, integer *, ftnlen),
37 cpotf2_(char *, integer *, complex *, integer *, integer *,
38 ftnlen), xerbla_(char *, integer *, ftnlen);
39 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
40 integer *, integer *, ftnlen, ftnlen);
41
42
43 /* -- LAPACK routine (version 3.0) -- */
44 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
45 /* Courant Institute, Argonne National Lab, and Rice University */
46 /* September 30, 1994 */
47
48 /* .. Scalar Arguments .. */
49 /* .. */
50 /* .. Array Arguments .. */
51 /* .. */
52
53 /* Purpose */
54 /* ======= */
55
56 /* CPBTRF computes the Cholesky factorization of a complex Hermitian */
57 /* positive definite band matrix A. */
58
59 /* The factorization has the form */
60 /* A = U**H * U, if UPLO = 'U', or */
61 /* A = L * L**H, if UPLO = 'L', */
62 /* where U is an upper triangular matrix and L is lower triangular. */
63
64 /* Arguments */
65 /* ========= */
66
67 /* UPLO (input) CHARACTER*1 */
68 /* = 'U': Upper triangle of A is stored; */
69 /* = 'L': Lower triangle of A is stored. */
70
71 /* N (input) INTEGER */
72 /* The order of the matrix A. N >= 0. */
73
74 /* KD (input) INTEGER */
75 /* The number of superdiagonals of the matrix A if UPLO = 'U', */
76 /* or the number of subdiagonals if UPLO = 'L'. KD >= 0. */
77
78 /* AB (input/output) COMPLEX array, dimension (LDAB,N) */
79 /* On entry, the upper or lower triangle of the Hermitian band */
80 /* matrix A, stored in the first KD+1 rows of the array. The */
81 /* j-th column of A is stored in the j-th column of the array AB */
82 /* as follows: */
83 /* if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; */
84 /* if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). */
85
86 /* On exit, if INFO = 0, the triangular factor U or L from the */
87 /* Cholesky factorization A = U**H*U or A = L*L**H of the band */
88 /* matrix A, in the same storage format as A. */
89
90 /* LDAB (input) INTEGER */
91 /* The leading dimension of the array AB. LDAB >= KD+1. */
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, the leading minor of order i is not */
97 /* positive definite, and the factorization could not be */
98 /* completed. */
99
100 /* Further Details */
101 /* =============== */
102
103 /* The band storage scheme is illustrated by the following example, when */
104 /* N = 6, KD = 2, and UPLO = 'U': */
105
106 /* On entry: On exit: */
107
108 /* * * a13 a24 a35 a46 * * u13 u24 u35 u46 */
109 /* * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 */
110 /* a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 */
111
112 /* Similarly, if UPLO = 'L' the format of A is as follows: */
113
114 /* On entry: On exit: */
115
116 /* a11 a22 a33 a44 a55 a66 l11 l22 l33 l44 l55 l66 */
117 /* a21 a32 a43 a54 a65 * l21 l32 l43 l54 l65 * */
118 /* a31 a42 a53 a64 * * l31 l42 l53 l64 * * */
119
120 /* Array elements marked * are not used by the routine. */
121
122 /* Contributed by */
123 /* Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989 */
124
125 /* ===================================================================== */
126
127 /* .. Parameters .. */
128 /* .. */
129 /* .. Local Scalars .. */
130 /* .. */
131 /* .. Local Arrays .. */
132 /* .. */
133 /* .. External Functions .. */
134 /* .. */
135 /* .. External Subroutines .. */
136 /* .. */
137 /* .. Intrinsic Functions .. */
138 /* .. */
139 /* .. Executable Statements .. */
140
141 /* Test the input parameters. */
142
143 /* Parameter adjustments */
144 ab_dim1 = *ldab;
145 ab_offset = 1 + ab_dim1;
146 ab -= ab_offset;
147
148 /* Function Body */
149 *info = 0;
150 if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
151 ftnlen)1, (ftnlen)1)) {
152 *info = -1;
153 } else if (*n < 0) {
154 *info = -2;
155 } else if (*kd < 0) {
156 *info = -3;
157 } else if (*ldab < *kd + 1) {
158 *info = -5;
159 }
160 if (*info != 0) {
161 i__1 = -(*info);
162 xerbla_("CPBTRF", &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 /* Determine the block size for this environment */
173
174 nb = ilaenv_(&c__1, "CPBTRF", uplo, n, kd, &c_n1, &c_n1, (ftnlen)6, (
175 ftnlen)1);
176
177 /* The block size must not exceed the semi-bandwidth KD, and must not */
178 /* exceed the limit set by the size of the local array WORK. */
179
180 nb = min(nb,32);
181
182 if (nb <= 1 || nb > *kd) {
183
184 /* Use unblocked code */
185
186 cpbtf2_(uplo, n, kd, &ab[ab_offset], ldab, info, (ftnlen)1);
187 } else {
188
189 /* Use blocked code */
190
191 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
192
193 /* Compute the Cholesky factorization of a Hermitian band */
194 /* matrix, given the upper triangle of the matrix in band */
195 /* storage. */
196
197 /* Zero the upper triangle of the work array. */
198
199 i__1 = nb;
200 for (j = 1; j <= i__1; ++j) {
201 i__2 = j - 1;
202 for (i__ = 1; i__ <= i__2; ++i__) {
203 i__3 = i__ + j * 33 - 34;
204 work[i__3].r = 0.f, work[i__3].i = 0.f;
205 /* L10: */
206 }
207 /* L20: */
208 }
209
210 /* Process the band matrix one diagonal block at a time. */
211
212 i__1 = *n;
213 i__2 = nb;
214 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
215 /* Computing MIN */
216 i__3 = nb, i__4 = *n - i__ + 1;
217 ib = min(i__3,i__4);
218
219 /* Factorize the diagonal block */
220
221 i__3 = *ldab - 1;
222 cpotf2_(uplo, &ib, &ab[*kd + 1 + i__ * ab_dim1], &i__3, &ii, (
223 ftnlen)1);
224 if (ii != 0) {
225 *info = i__ + ii - 1;
226 goto L150;
227 }
228 if (i__ + ib <= *n) {
229
230 /* Update the relevant part of the trailing submatrix. */
231 /* If A11 denotes the diagonal block which has just been */
232 /* factorized, then we need to update the remaining */
233 /* blocks in the diagram: */
234
235 /* A11 A12 A13 */
236 /* A22 A23 */
237 /* A33 */
238
239 /* The numbers of rows and columns in the partitioning */
240 /* are IB, I2, I3 respectively. The blocks A12, A22 and */
241 /* A23 are empty if IB = KD. The upper triangle of A13 */
242 /* lies outside the band. */
243
244 /* Computing MIN */
245 i__3 = *kd - ib, i__4 = *n - i__ - ib + 1;
246 i2 = min(i__3,i__4);
247 /* Computing MIN */
248 i__3 = ib, i__4 = *n - i__ - *kd + 1;
249 i3 = min(i__3,i__4);
250
251 if (i2 > 0) {
252
253 /* Update A12 */
254
255 i__3 = *ldab - 1;
256 i__4 = *ldab - 1;
257 ctrsm_("Left", "Upper", "Conjugate transpose", "Non-"
258 "unit", &ib, &i2, &c_b1, &ab[*kd + 1 + i__ *
259 ab_dim1], &i__3, &ab[*kd + 1 - ib + (i__ + ib)
260 * ab_dim1], &i__4, (ftnlen)4, (ftnlen)5, (
261 ftnlen)19, (ftnlen)8);
262
263 /* Update A22 */
264
265 i__3 = *ldab - 1;
266 i__4 = *ldab - 1;
267 cherk_("Upper", "Conjugate transpose", &i2, &ib, &
268 c_b21, &ab[*kd + 1 - ib + (i__ + ib) *
269 ab_dim1], &i__3, &c_b22, &ab[*kd + 1 + (i__ +
270 ib) * ab_dim1], &i__4, (ftnlen)5, (ftnlen)19);
271 }
272
273 if (i3 > 0) {
274
275 /* Copy the lower triangle of A13 into the work array. */
276
277 i__3 = i3;
278 for (jj = 1; jj <= i__3; ++jj) {
279 i__4 = ib;
280 for (ii = jj; ii <= i__4; ++ii) {
281 i__5 = ii + jj * 33 - 34;
282 i__6 = ii - jj + 1 + (jj + i__ + *kd - 1) *
283 ab_dim1;
284 work[i__5].r = ab[i__6].r, work[i__5].i = ab[
285 i__6].i;
286 /* L30: */
287 }
288 /* L40: */
289 }
290
291 /* Update A13 (in the work array). */
292
293 i__3 = *ldab - 1;
294 ctrsm_("Left", "Upper", "Conjugate transpose", "Non-"
295 "unit", &ib, &i3, &c_b1, &ab[*kd + 1 + i__ *
296 ab_dim1], &i__3, work, &c__33, (ftnlen)4, (
297 ftnlen)5, (ftnlen)19, (ftnlen)8);
298
299 /* Update A23 */
300
301 if (i2 > 0) {
302 q__1.r = -1.f, q__1.i = -0.f;
303 i__3 = *ldab - 1;
304 i__4 = *ldab - 1;
305 cgemm_("Conjugate transpose", "No transpose", &i2,
306 &i3, &ib, &q__1, &ab[*kd + 1 - ib + (i__
307 + ib) * ab_dim1], &i__3, work, &c__33, &
308 c_b1, &ab[ib + 1 + (i__ + *kd) * ab_dim1],
309 &i__4, (ftnlen)19, (ftnlen)12);
310 }
311
312 /* Update A33 */
313
314 i__3 = *ldab - 1;
315 cherk_("Upper", "Conjugate transpose", &i3, &ib, &
316 c_b21, work, &c__33, &c_b22, &ab[*kd + 1 + (
317 i__ + *kd) * ab_dim1], &i__3, (ftnlen)5, (
318 ftnlen)19);
319
320 /* Copy the lower triangle of A13 back into place. */
321
322 i__3 = i3;
323 for (jj = 1; jj <= i__3; ++jj) {
324 i__4 = ib;
325 for (ii = jj; ii <= i__4; ++ii) {
326 i__5 = ii - jj + 1 + (jj + i__ + *kd - 1) *
327 ab_dim1;
328 i__6 = ii + jj * 33 - 34;
329 ab[i__5].r = work[i__6].r, ab[i__5].i = work[
330 i__6].i;
331 /* L50: */
332 }
333 /* L60: */
334 }
335 }
336 }
337 /* L70: */
338 }
339 } else {
340
341 /* Compute the Cholesky factorization of a Hermitian band */
342 /* matrix, given the lower triangle of the matrix in band */
343 /* storage. */
344
345 /* Zero the lower triangle of the work array. */
346
347 i__2 = nb;
348 for (j = 1; j <= i__2; ++j) {
349 i__1 = nb;
350 for (i__ = j + 1; i__ <= i__1; ++i__) {
351 i__3 = i__ + j * 33 - 34;
352 work[i__3].r = 0.f, work[i__3].i = 0.f;
353 /* L80: */
354 }
355 /* L90: */
356 }
357
358 /* Process the band matrix one diagonal block at a time. */
359
360 i__2 = *n;
361 i__1 = nb;
362 for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) {
363 /* Computing MIN */
364 i__3 = nb, i__4 = *n - i__ + 1;
365 ib = min(i__3,i__4);
366
367 /* Factorize the diagonal block */
368
369 i__3 = *ldab - 1;
370 cpotf2_(uplo, &ib, &ab[i__ * ab_dim1 + 1], &i__3, &ii, (
371 ftnlen)1);
372 if (ii != 0) {
373 *info = i__ + ii - 1;
374 goto L150;
375 }
376 if (i__ + ib <= *n) {
377
378 /* Update the relevant part of the trailing submatrix. */
379 /* If A11 denotes the diagonal block which has just been */
380 /* factorized, then we need to update the remaining */
381 /* blocks in the diagram: */
382
383 /* A11 */
384 /* A21 A22 */
385 /* A31 A32 A33 */
386
387 /* The numbers of rows and columns in the partitioning */
388 /* are IB, I2, I3 respectively. The blocks A21, A22 and */
389 /* A32 are empty if IB = KD. The lower triangle of A31 */
390 /* lies outside the band. */
391
392 /* Computing MIN */
393 i__3 = *kd - ib, i__4 = *n - i__ - ib + 1;
394 i2 = min(i__3,i__4);
395 /* Computing MIN */
396 i__3 = ib, i__4 = *n - i__ - *kd + 1;
397 i3 = min(i__3,i__4);
398
399 if (i2 > 0) {
400
401 /* Update A21 */
402
403 i__3 = *ldab - 1;
404 i__4 = *ldab - 1;
405 ctrsm_("Right", "Lower", "Conjugate transpose", "Non"
406 "-unit", &i2, &ib, &c_b1, &ab[i__ * ab_dim1 +
407 1], &i__3, &ab[ib + 1 + i__ * ab_dim1], &i__4,
408 (ftnlen)5, (ftnlen)5, (ftnlen)19, (ftnlen)8);
409
410 /* Update A22 */
411
412 i__3 = *ldab - 1;
413 i__4 = *ldab - 1;
414 cherk_("Lower", "No transpose", &i2, &ib, &c_b21, &ab[
415 ib + 1 + i__ * ab_dim1], &i__3, &c_b22, &ab[(
416 i__ + ib) * ab_dim1 + 1], &i__4, (ftnlen)5, (
417 ftnlen)12);
418 }
419
420 if (i3 > 0) {
421
422 /* Copy the upper triangle of A31 into the work array. */
423
424 i__3 = ib;
425 for (jj = 1; jj <= i__3; ++jj) {
426 i__4 = min(jj,i3);
427 for (ii = 1; ii <= i__4; ++ii) {
428 i__5 = ii + jj * 33 - 34;
429 i__6 = *kd + 1 - jj + ii + (jj + i__ - 1) *
430 ab_dim1;
431 work[i__5].r = ab[i__6].r, work[i__5].i = ab[
432 i__6].i;
433 /* L100: */
434 }
435 /* L110: */
436 }
437
438 /* Update A31 (in the work array). */
439
440 i__3 = *ldab - 1;
441 ctrsm_("Right", "Lower", "Conjugate transpose", "Non"
442 "-unit", &i3, &ib, &c_b1, &ab[i__ * ab_dim1 +
443 1], &i__3, work, &c__33, (ftnlen)5, (ftnlen)5,
444 (ftnlen)19, (ftnlen)8);
445
446 /* Update A32 */
447
448 if (i2 > 0) {
449 q__1.r = -1.f, q__1.i = -0.f;
450 i__3 = *ldab - 1;
451 i__4 = *ldab - 1;
452 cgemm_("No transpose", "Conjugate transpose", &i3,
453 &i2, &ib, &q__1, work, &c__33, &ab[ib +
454 1 + i__ * ab_dim1], &i__3, &c_b1, &ab[*kd
455 + 1 - ib + (i__ + ib) * ab_dim1], &i__4, (
456 ftnlen)12, (ftnlen)19);
457 }
458
459 /* Update A33 */
460
461 i__3 = *ldab - 1;
462 cherk_("Lower", "No transpose", &i3, &ib, &c_b21,
463 work, &c__33, &c_b22, &ab[(i__ + *kd) *
464 ab_dim1 + 1], &i__3, (ftnlen)5, (ftnlen)12);
465
466 /* Copy the upper triangle of A31 back into place. */
467
468 i__3 = ib;
469 for (jj = 1; jj <= i__3; ++jj) {
470 i__4 = min(jj,i3);
471 for (ii = 1; ii <= i__4; ++ii) {
472 i__5 = *kd + 1 - jj + ii + (jj + i__ - 1) *
473 ab_dim1;
474 i__6 = ii + jj * 33 - 34;
475 ab[i__5].r = work[i__6].r, ab[i__5].i = work[
476 i__6].i;
477 /* L120: */
478 }
479 /* L130: */
480 }
481 }
482 }
483 /* L140: */
484 }
485 }
486 }
487 return 0;
488
489 L150:
490 return 0;
491
492 /* End of CPBTRF */
493
494 } /* cpbtrf_ */
495
496