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