001:       SUBROUTINE DPBTRF( UPLO, N, KD, AB, LDAB, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          UPLO
010:       INTEGER            INFO, KD, LDAB, N
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   AB( LDAB, * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  DPBTRF computes the Cholesky factorization of a real symmetric
020: *  positive definite band matrix A.
021: *
022: *  The factorization has the form
023: *     A = U**T * U,  if UPLO = 'U', or
024: *     A = L  * L**T,  if UPLO = 'L',
025: *  where U is an upper triangular matrix and L is lower triangular.
026: *
027: *  Arguments
028: *  =========
029: *
030: *  UPLO    (input) CHARACTER*1
031: *          = 'U':  Upper triangle of A is stored;
032: *          = 'L':  Lower triangle of A is stored.
033: *
034: *  N       (input) INTEGER
035: *          The order of the matrix A.  N >= 0.
036: *
037: *  KD      (input) INTEGER
038: *          The number of superdiagonals of the matrix A if UPLO = 'U',
039: *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
040: *
041: *  AB      (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
042: *          On entry, the upper or lower triangle of the symmetric band
043: *          matrix A, stored in the first KD+1 rows of the array.  The
044: *          j-th column of A is stored in the j-th column of the array AB
045: *          as follows:
046: *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
047: *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
048: *
049: *          On exit, if INFO = 0, the triangular factor U or L from the
050: *          Cholesky factorization A = U**T*U or A = L*L**T of the band
051: *          matrix A, in the same storage format as A.
052: *
053: *  LDAB    (input) INTEGER
054: *          The leading dimension of the array AB.  LDAB >= KD+1.
055: *
056: *  INFO    (output) INTEGER
057: *          = 0:  successful exit
058: *          < 0:  if INFO = -i, the i-th argument had an illegal value
059: *          > 0:  if INFO = i, the leading minor of order i is not
060: *                positive definite, and the factorization could not be
061: *                completed.
062: *
063: *  Further Details
064: *  ===============
065: *
066: *  The band storage scheme is illustrated by the following example, when
067: *  N = 6, KD = 2, and UPLO = 'U':
068: *
069: *  On entry:                       On exit:
070: *
071: *      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
072: *      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
073: *     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
074: *
075: *  Similarly, if UPLO = 'L' the format of A is as follows:
076: *
077: *  On entry:                       On exit:
078: *
079: *     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
080: *     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
081: *     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *
082: *
083: *  Array elements marked * are not used by the routine.
084: *
085: *  Contributed by
086: *  Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989
087: *
088: *  =====================================================================
089: *
090: *     .. Parameters ..
091:       DOUBLE PRECISION   ONE, ZERO
092:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
093:       INTEGER            NBMAX, LDWORK
094:       PARAMETER          ( NBMAX = 32, LDWORK = NBMAX+1 )
095: *     ..
096: *     .. Local Scalars ..
097:       INTEGER            I, I2, I3, IB, II, J, JJ, NB
098: *     ..
099: *     .. Local Arrays ..
100:       DOUBLE PRECISION   WORK( LDWORK, NBMAX )
101: *     ..
102: *     .. External Functions ..
103:       LOGICAL            LSAME
104:       INTEGER            ILAENV
105:       EXTERNAL           LSAME, ILAENV
106: *     ..
107: *     .. External Subroutines ..
108:       EXTERNAL           DGEMM, DPBTF2, DPOTF2, DSYRK, DTRSM, XERBLA
109: *     ..
110: *     .. Intrinsic Functions ..
111:       INTRINSIC          MIN
112: *     ..
113: *     .. Executable Statements ..
114: *
115: *     Test the input parameters.
116: *
117:       INFO = 0
118:       IF( ( .NOT.LSAME( UPLO, 'U' ) ) .AND.
119:      $    ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN
120:          INFO = -1
121:       ELSE IF( N.LT.0 ) THEN
122:          INFO = -2
123:       ELSE IF( KD.LT.0 ) THEN
124:          INFO = -3
125:       ELSE IF( LDAB.LT.KD+1 ) THEN
126:          INFO = -5
127:       END IF
128:       IF( INFO.NE.0 ) THEN
129:          CALL XERBLA( 'DPBTRF', -INFO )
130:          RETURN
131:       END IF
132: *
133: *     Quick return if possible
134: *
135:       IF( N.EQ.0 )
136:      $   RETURN
137: *
138: *     Determine the block size for this environment
139: *
140:       NB = ILAENV( 1, 'DPBTRF', UPLO, N, KD, -1, -1 )
141: *
142: *     The block size must not exceed the semi-bandwidth KD, and must not
143: *     exceed the limit set by the size of the local array WORK.
144: *
145:       NB = MIN( NB, NBMAX )
146: *
147:       IF( NB.LE.1 .OR. NB.GT.KD ) THEN
148: *
149: *        Use unblocked code
150: *
151:          CALL DPBTF2( UPLO, N, KD, AB, LDAB, INFO )
152:       ELSE
153: *
154: *        Use blocked code
155: *
156:          IF( LSAME( UPLO, 'U' ) ) THEN
157: *
158: *           Compute the Cholesky factorization of a symmetric band
159: *           matrix, given the upper triangle of the matrix in band
160: *           storage.
161: *
162: *           Zero the upper triangle of the work array.
163: *
164:             DO 20 J = 1, NB
165:                DO 10 I = 1, J - 1
166:                   WORK( I, J ) = ZERO
167:    10          CONTINUE
168:    20       CONTINUE
169: *
170: *           Process the band matrix one diagonal block at a time.
171: *
172:             DO 70 I = 1, N, NB
173:                IB = MIN( NB, N-I+1 )
174: *
175: *              Factorize the diagonal block
176: *
177:                CALL DPOTF2( UPLO, IB, AB( KD+1, I ), LDAB-1, II )
178:                IF( II.NE.0 ) THEN
179:                   INFO = I + II - 1
180:                   GO TO 150
181:                END IF
182:                IF( I+IB.LE.N ) THEN
183: *
184: *                 Update the relevant part of the trailing submatrix.
185: *                 If A11 denotes the diagonal block which has just been
186: *                 factorized, then we need to update the remaining
187: *                 blocks in the diagram:
188: *
189: *                    A11   A12   A13
190: *                          A22   A23
191: *                                A33
192: *
193: *                 The numbers of rows and columns in the partitioning
194: *                 are IB, I2, I3 respectively. The blocks A12, A22 and
195: *                 A23 are empty if IB = KD. The upper triangle of A13
196: *                 lies outside the band.
197: *
198:                   I2 = MIN( KD-IB, N-I-IB+1 )
199:                   I3 = MIN( IB, N-I-KD+1 )
200: *
201:                   IF( I2.GT.0 ) THEN
202: *
203: *                    Update A12
204: *
205:                      CALL DTRSM( 'Left', 'Upper', 'Transpose',
206:      $                           'Non-unit', IB, I2, ONE, AB( KD+1, I ),
207:      $                           LDAB-1, AB( KD+1-IB, I+IB ), LDAB-1 )
208: *
209: *                    Update A22
210: *
211:                      CALL DSYRK( 'Upper', 'Transpose', I2, IB, -ONE,
212:      $                           AB( KD+1-IB, I+IB ), LDAB-1, ONE,
213:      $                           AB( KD+1, I+IB ), LDAB-1 )
214:                   END IF
215: *
216:                   IF( I3.GT.0 ) THEN
217: *
218: *                    Copy the lower triangle of A13 into the work array.
219: *
220:                      DO 40 JJ = 1, I3
221:                         DO 30 II = JJ, IB
222:                            WORK( II, JJ ) = AB( II-JJ+1, JJ+I+KD-1 )
223:    30                   CONTINUE
224:    40                CONTINUE
225: *
226: *                    Update A13 (in the work array).
227: *
228:                      CALL DTRSM( 'Left', 'Upper', 'Transpose',
229:      $                           'Non-unit', IB, I3, ONE, AB( KD+1, I ),
230:      $                           LDAB-1, WORK, LDWORK )
231: *
232: *                    Update A23
233: *
234:                      IF( I2.GT.0 )
235:      $                  CALL DGEMM( 'Transpose', 'No Transpose', I2, I3,
236:      $                              IB, -ONE, AB( KD+1-IB, I+IB ),
237:      $                              LDAB-1, WORK, LDWORK, ONE,
238:      $                              AB( 1+IB, I+KD ), LDAB-1 )
239: *
240: *                    Update A33
241: *
242:                      CALL DSYRK( 'Upper', 'Transpose', I3, IB, -ONE,
243:      $                           WORK, LDWORK, ONE, AB( KD+1, I+KD ),
244:      $                           LDAB-1 )
245: *
246: *                    Copy the lower triangle of A13 back into place.
247: *
248:                      DO 60 JJ = 1, I3
249:                         DO 50 II = JJ, IB
250:                            AB( II-JJ+1, JJ+I+KD-1 ) = WORK( II, JJ )
251:    50                   CONTINUE
252:    60                CONTINUE
253:                   END IF
254:                END IF
255:    70       CONTINUE
256:          ELSE
257: *
258: *           Compute the Cholesky factorization of a symmetric band
259: *           matrix, given the lower triangle of the matrix in band
260: *           storage.
261: *
262: *           Zero the lower triangle of the work array.
263: *
264:             DO 90 J = 1, NB
265:                DO 80 I = J + 1, NB
266:                   WORK( I, J ) = ZERO
267:    80          CONTINUE
268:    90       CONTINUE
269: *
270: *           Process the band matrix one diagonal block at a time.
271: *
272:             DO 140 I = 1, N, NB
273:                IB = MIN( NB, N-I+1 )
274: *
275: *              Factorize the diagonal block
276: *
277:                CALL DPOTF2( UPLO, IB, AB( 1, I ), LDAB-1, II )
278:                IF( II.NE.0 ) THEN
279:                   INFO = I + II - 1
280:                   GO TO 150
281:                END IF
282:                IF( I+IB.LE.N ) THEN
283: *
284: *                 Update the relevant part of the trailing submatrix.
285: *                 If A11 denotes the diagonal block which has just been
286: *                 factorized, then we need to update the remaining
287: *                 blocks in the diagram:
288: *
289: *                    A11
290: *                    A21   A22
291: *                    A31   A32   A33
292: *
293: *                 The numbers of rows and columns in the partitioning
294: *                 are IB, I2, I3 respectively. The blocks A21, A22 and
295: *                 A32 are empty if IB = KD. The lower triangle of A31
296: *                 lies outside the band.
297: *
298:                   I2 = MIN( KD-IB, N-I-IB+1 )
299:                   I3 = MIN( IB, N-I-KD+1 )
300: *
301:                   IF( I2.GT.0 ) THEN
302: *
303: *                    Update A21
304: *
305:                      CALL DTRSM( 'Right', 'Lower', 'Transpose',
306:      $                           'Non-unit', I2, IB, ONE, AB( 1, I ),
307:      $                           LDAB-1, AB( 1+IB, I ), LDAB-1 )
308: *
309: *                    Update A22
310: *
311:                      CALL DSYRK( 'Lower', 'No Transpose', I2, IB, -ONE,
312:      $                           AB( 1+IB, I ), LDAB-1, ONE,
313:      $                           AB( 1, I+IB ), LDAB-1 )
314:                   END IF
315: *
316:                   IF( I3.GT.0 ) THEN
317: *
318: *                    Copy the upper triangle of A31 into the work array.
319: *
320:                      DO 110 JJ = 1, IB
321:                         DO 100 II = 1, MIN( JJ, I3 )
322:                            WORK( II, JJ ) = AB( KD+1-JJ+II, JJ+I-1 )
323:   100                   CONTINUE
324:   110                CONTINUE
325: *
326: *                    Update A31 (in the work array).
327: *
328:                      CALL DTRSM( 'Right', 'Lower', 'Transpose',
329:      $                           'Non-unit', I3, IB, ONE, AB( 1, I ),
330:      $                           LDAB-1, WORK, LDWORK )
331: *
332: *                    Update A32
333: *
334:                      IF( I2.GT.0 )
335:      $                  CALL DGEMM( 'No transpose', 'Transpose', I3, I2,
336:      $                              IB, -ONE, WORK, LDWORK,
337:      $                              AB( 1+IB, I ), LDAB-1, ONE,
338:      $                              AB( 1+KD-IB, I+IB ), LDAB-1 )
339: *
340: *                    Update A33
341: *
342:                      CALL DSYRK( 'Lower', 'No Transpose', I3, IB, -ONE,
343:      $                           WORK, LDWORK, ONE, AB( 1, I+KD ),
344:      $                           LDAB-1 )
345: *
346: *                    Copy the upper triangle of A31 back into place.
347: *
348:                      DO 130 JJ = 1, IB
349:                         DO 120 II = 1, MIN( JJ, I3 )
350:                            AB( KD+1-JJ+II, JJ+I-1 ) = WORK( II, JJ )
351:   120                   CONTINUE
352:   130                CONTINUE
353:                   END IF
354:                END IF
355:   140       CONTINUE
356:          END IF
357:       END IF
358:       RETURN
359: *
360:   150 CONTINUE
361:       RETURN
362: *
363: *     End of DPBTRF
364: *
365:       END
366: