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