001:       SUBROUTINE ZPBTRF( 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:       COMPLEX*16         AB( LDAB, * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  ZPBTRF computes the Cholesky factorization of a complex Hermitian
020: *  positive definite band matrix A.
021: *
022: *  The factorization has the form
023: *     A = U**H * U,  if UPLO = 'U', or
024: *     A = L  * L**H,  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) COMPLEX*16 array, dimension (LDAB,N)
042: *          On entry, the upper or lower triangle of the Hermitian 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**H*U or A = L*L**H 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:       COMPLEX*16         CONE
094:       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
095:       INTEGER            NBMAX, LDWORK
096:       PARAMETER          ( NBMAX = 32, LDWORK = NBMAX+1 )
097: *     ..
098: *     .. Local Scalars ..
099:       INTEGER            I, I2, I3, IB, II, J, JJ, NB
100: *     ..
101: *     .. Local Arrays ..
102:       COMPLEX*16         WORK( LDWORK, NBMAX )
103: *     ..
104: *     .. External Functions ..
105:       LOGICAL            LSAME
106:       INTEGER            ILAENV
107:       EXTERNAL           LSAME, ILAENV
108: *     ..
109: *     .. External Subroutines ..
110:       EXTERNAL           XERBLA, ZGEMM, ZHERK, ZPBTF2, ZPOTF2, ZTRSM
111: *     ..
112: *     .. Intrinsic Functions ..
113:       INTRINSIC          MIN
114: *     ..
115: *     .. Executable Statements ..
116: *
117: *     Test the input parameters.
118: *
119:       INFO = 0
120:       IF( ( .NOT.LSAME( UPLO, 'U' ) ) .AND.
121:      $    ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN
122:          INFO = -1
123:       ELSE IF( N.LT.0 ) THEN
124:          INFO = -2
125:       ELSE IF( KD.LT.0 ) THEN
126:          INFO = -3
127:       ELSE IF( LDAB.LT.KD+1 ) THEN
128:          INFO = -5
129:       END IF
130:       IF( INFO.NE.0 ) THEN
131:          CALL XERBLA( 'ZPBTRF', -INFO )
132:          RETURN
133:       END IF
134: *
135: *     Quick return if possible
136: *
137:       IF( N.EQ.0 )
138:      $   RETURN
139: *
140: *     Determine the block size for this environment
141: *
142:       NB = ILAENV( 1, 'ZPBTRF', UPLO, N, KD, -1, -1 )
143: *
144: *     The block size must not exceed the semi-bandwidth KD, and must not
145: *     exceed the limit set by the size of the local array WORK.
146: *
147:       NB = MIN( NB, NBMAX )
148: *
149:       IF( NB.LE.1 .OR. NB.GT.KD ) THEN
150: *
151: *        Use unblocked code
152: *
153:          CALL ZPBTF2( UPLO, N, KD, AB, LDAB, INFO )
154:       ELSE
155: *
156: *        Use blocked code
157: *
158:          IF( LSAME( UPLO, 'U' ) ) THEN
159: *
160: *           Compute the Cholesky factorization of a Hermitian band
161: *           matrix, given the upper triangle of the matrix in band
162: *           storage.
163: *
164: *           Zero the upper triangle of the work array.
165: *
166:             DO 20 J = 1, NB
167:                DO 10 I = 1, J - 1
168:                   WORK( I, J ) = ZERO
169:    10          CONTINUE
170:    20       CONTINUE
171: *
172: *           Process the band matrix one diagonal block at a time.
173: *
174:             DO 70 I = 1, N, NB
175:                IB = MIN( NB, N-I+1 )
176: *
177: *              Factorize the diagonal block
178: *
179:                CALL ZPOTF2( UPLO, IB, AB( KD+1, I ), LDAB-1, II )
180:                IF( II.NE.0 ) THEN
181:                   INFO = I + II - 1
182:                   GO TO 150
183:                END IF
184:                IF( I+IB.LE.N ) THEN
185: *
186: *                 Update the relevant part of the trailing submatrix.
187: *                 If A11 denotes the diagonal block which has just been
188: *                 factorized, then we need to update the remaining
189: *                 blocks in the diagram:
190: *
191: *                    A11   A12   A13
192: *                          A22   A23
193: *                                A33
194: *
195: *                 The numbers of rows and columns in the partitioning
196: *                 are IB, I2, I3 respectively. The blocks A12, A22 and
197: *                 A23 are empty if IB = KD. The upper triangle of A13
198: *                 lies outside the band.
199: *
200:                   I2 = MIN( KD-IB, N-I-IB+1 )
201:                   I3 = MIN( IB, N-I-KD+1 )
202: *
203:                   IF( I2.GT.0 ) THEN
204: *
205: *                    Update A12
206: *
207:                      CALL ZTRSM( 'Left', 'Upper', 'Conjugate transpose',
208:      $                           'Non-unit', IB, I2, CONE,
209:      $                           AB( KD+1, I ), LDAB-1,
210:      $                           AB( KD+1-IB, I+IB ), LDAB-1 )
211: *
212: *                    Update A22
213: *
214:                      CALL ZHERK( 'Upper', 'Conjugate transpose', I2, IB,
215:      $                           -ONE, AB( KD+1-IB, I+IB ), LDAB-1, ONE,
216:      $                           AB( KD+1, I+IB ), LDAB-1 )
217:                   END IF
218: *
219:                   IF( I3.GT.0 ) THEN
220: *
221: *                    Copy the lower triangle of A13 into the work array.
222: *
223:                      DO 40 JJ = 1, I3
224:                         DO 30 II = JJ, IB
225:                            WORK( II, JJ ) = AB( II-JJ+1, JJ+I+KD-1 )
226:    30                   CONTINUE
227:    40                CONTINUE
228: *
229: *                    Update A13 (in the work array).
230: *
231:                      CALL ZTRSM( 'Left', 'Upper', 'Conjugate transpose',
232:      $                           'Non-unit', IB, I3, CONE,
233:      $                           AB( KD+1, I ), LDAB-1, WORK, LDWORK )
234: *
235: *                    Update A23
236: *
237:                      IF( I2.GT.0 )
238:      $                  CALL ZGEMM( 'Conjugate transpose',
239:      $                              'No transpose', I2, I3, IB, -CONE,
240:      $                              AB( KD+1-IB, I+IB ), LDAB-1, WORK,
241:      $                              LDWORK, CONE, AB( 1+IB, I+KD ),
242:      $                              LDAB-1 )
243: *
244: *                    Update A33
245: *
246:                      CALL ZHERK( 'Upper', 'Conjugate transpose', I3, IB,
247:      $                           -ONE, WORK, LDWORK, ONE,
248:      $                           AB( KD+1, I+KD ), LDAB-1 )
249: *
250: *                    Copy the lower triangle of A13 back into place.
251: *
252:                      DO 60 JJ = 1, I3
253:                         DO 50 II = JJ, IB
254:                            AB( II-JJ+1, JJ+I+KD-1 ) = WORK( II, JJ )
255:    50                   CONTINUE
256:    60                CONTINUE
257:                   END IF
258:                END IF
259:    70       CONTINUE
260:          ELSE
261: *
262: *           Compute the Cholesky factorization of a Hermitian band
263: *           matrix, given the lower triangle of the matrix in band
264: *           storage.
265: *
266: *           Zero the lower triangle of the work array.
267: *
268:             DO 90 J = 1, NB
269:                DO 80 I = J + 1, NB
270:                   WORK( I, J ) = ZERO
271:    80          CONTINUE
272:    90       CONTINUE
273: *
274: *           Process the band matrix one diagonal block at a time.
275: *
276:             DO 140 I = 1, N, NB
277:                IB = MIN( NB, N-I+1 )
278: *
279: *              Factorize the diagonal block
280: *
281:                CALL ZPOTF2( UPLO, IB, AB( 1, I ), LDAB-1, II )
282:                IF( II.NE.0 ) THEN
283:                   INFO = I + II - 1
284:                   GO TO 150
285:                END IF
286:                IF( I+IB.LE.N ) THEN
287: *
288: *                 Update the relevant part of the trailing submatrix.
289: *                 If A11 denotes the diagonal block which has just been
290: *                 factorized, then we need to update the remaining
291: *                 blocks in the diagram:
292: *
293: *                    A11
294: *                    A21   A22
295: *                    A31   A32   A33
296: *
297: *                 The numbers of rows and columns in the partitioning
298: *                 are IB, I2, I3 respectively. The blocks A21, A22 and
299: *                 A32 are empty if IB = KD. The lower triangle of A31
300: *                 lies outside the band.
301: *
302:                   I2 = MIN( KD-IB, N-I-IB+1 )
303:                   I3 = MIN( IB, N-I-KD+1 )
304: *
305:                   IF( I2.GT.0 ) THEN
306: *
307: *                    Update A21
308: *
309:                      CALL ZTRSM( 'Right', 'Lower',
310:      $                           'Conjugate transpose', 'Non-unit', I2,
311:      $                           IB, CONE, AB( 1, I ), LDAB-1,
312:      $                           AB( 1+IB, I ), LDAB-1 )
313: *
314: *                    Update A22
315: *
316:                      CALL ZHERK( 'Lower', 'No transpose', I2, IB, -ONE,
317:      $                           AB( 1+IB, I ), LDAB-1, ONE,
318:      $                           AB( 1, I+IB ), LDAB-1 )
319:                   END IF
320: *
321:                   IF( I3.GT.0 ) THEN
322: *
323: *                    Copy the upper triangle of A31 into the work array.
324: *
325:                      DO 110 JJ = 1, IB
326:                         DO 100 II = 1, MIN( JJ, I3 )
327:                            WORK( II, JJ ) = AB( KD+1-JJ+II, JJ+I-1 )
328:   100                   CONTINUE
329:   110                CONTINUE
330: *
331: *                    Update A31 (in the work array).
332: *
333:                      CALL ZTRSM( 'Right', 'Lower',
334:      $                           'Conjugate transpose', 'Non-unit', I3,
335:      $                           IB, CONE, AB( 1, I ), LDAB-1, WORK,
336:      $                           LDWORK )
337: *
338: *                    Update A32
339: *
340:                      IF( I2.GT.0 )
341:      $                  CALL ZGEMM( 'No transpose',
342:      $                              'Conjugate transpose', I3, I2, IB,
343:      $                              -CONE, WORK, LDWORK, AB( 1+IB, I ),
344:      $                              LDAB-1, CONE, AB( 1+KD-IB, I+IB ),
345:      $                              LDAB-1 )
346: *
347: *                    Update A33
348: *
349:                      CALL ZHERK( 'Lower', 'No transpose', I3, IB, -ONE,
350:      $                           WORK, LDWORK, ONE, AB( 1, I+KD ),
351:      $                           LDAB-1 )
352: *
353: *                    Copy the upper triangle of A31 back into place.
354: *
355:                      DO 130 JJ = 1, IB
356:                         DO 120 II = 1, MIN( JJ, I3 )
357:                            AB( KD+1-JJ+II, JJ+I-1 ) = WORK( II, JJ )
358:   120                   CONTINUE
359:   130                CONTINUE
360:                   END IF
361:                END IF
362:   140       CONTINUE
363:          END IF
364:       END IF
365:       RETURN
366: *
367:   150 CONTINUE
368:       RETURN
369: *
370: *     End of ZPBTRF
371: *
372:       END
373: