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