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