001:       SUBROUTINE CPBTF2( 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            AB( LDAB, * )
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  CPBTF2 computes the Cholesky factorization of a complex Hermitian
019: *  positive definite band matrix A.
020: *
021: *  The factorization has the form
022: *     A = U' * U ,  if UPLO = 'U', or
023: *     A = L  * L',  if UPLO = 'L',
024: *  where U is an upper triangular matrix, U' is the conjugate transpose
025: *  of U, and L is lower triangular.
026: *
027: *  This is the unblocked version of the algorithm, calling Level 2 BLAS.
028: *
029: *  Arguments
030: *  =========
031: *
032: *  UPLO    (input) CHARACTER*1
033: *          Specifies whether the upper or lower triangular part of the
034: *          Hermitian matrix A is stored:
035: *          = 'U':  Upper triangular
036: *          = 'L':  Lower triangular
037: *
038: *  N       (input) INTEGER
039: *          The order of the matrix A.  N >= 0.
040: *
041: *  KD      (input) INTEGER
042: *          The number of super-diagonals of the matrix A if UPLO = 'U',
043: *          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0.
044: *
045: *  AB      (input/output) COMPLEX array, dimension (LDAB,N)
046: *          On entry, the upper or lower triangle of the Hermitian band
047: *          matrix A, stored in the first KD+1 rows of the array.  The
048: *          j-th column of A is stored in the j-th column of the array AB
049: *          as follows:
050: *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
051: *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
052: *
053: *          On exit, if INFO = 0, the triangular factor U or L from the
054: *          Cholesky factorization A = U'*U or A = L*L' of the band
055: *          matrix A, in the same storage format as A.
056: *
057: *  LDAB    (input) INTEGER
058: *          The leading dimension of the array AB.  LDAB >= KD+1.
059: *
060: *  INFO    (output) INTEGER
061: *          = 0: successful exit
062: *          < 0: if INFO = -k, the k-th argument had an illegal value
063: *          > 0: if INFO = k, the leading minor of order k is not
064: *               positive definite, and the factorization could not be
065: *               completed.
066: *
067: *  Further Details
068: *  ===============
069: *
070: *  The band storage scheme is illustrated by the following example, when
071: *  N = 6, KD = 2, and UPLO = 'U':
072: *
073: *  On entry:                       On exit:
074: *
075: *      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
076: *      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
077: *     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
078: *
079: *  Similarly, if UPLO = 'L' the format of A is as follows:
080: *
081: *  On entry:                       On exit:
082: *
083: *     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
084: *     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
085: *     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *
086: *
087: *  Array elements marked * are not used by the routine.
088: *
089: *  =====================================================================
090: *
091: *     .. Parameters ..
092:       REAL               ONE, ZERO
093:       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
094: *     ..
095: *     .. Local Scalars ..
096:       LOGICAL            UPPER
097:       INTEGER            J, KLD, KN
098:       REAL               AJJ
099: *     ..
100: *     .. External Functions ..
101:       LOGICAL            LSAME
102:       EXTERNAL           LSAME
103: *     ..
104: *     .. External Subroutines ..
105:       EXTERNAL           CHER, CLACGV, CSSCAL, XERBLA
106: *     ..
107: *     .. Intrinsic Functions ..
108:       INTRINSIC          MAX, MIN, REAL, SQRT
109: *     ..
110: *     .. Executable Statements ..
111: *
112: *     Test the input parameters.
113: *
114:       INFO = 0
115:       UPPER = LSAME( UPLO, 'U' )
116:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
117:          INFO = -1
118:       ELSE IF( N.LT.0 ) THEN
119:          INFO = -2
120:       ELSE IF( KD.LT.0 ) THEN
121:          INFO = -3
122:       ELSE IF( LDAB.LT.KD+1 ) THEN
123:          INFO = -5
124:       END IF
125:       IF( INFO.NE.0 ) THEN
126:          CALL XERBLA( 'CPBTF2', -INFO )
127:          RETURN
128:       END IF
129: *
130: *     Quick return if possible
131: *
132:       IF( N.EQ.0 )
133:      $   RETURN
134: *
135:       KLD = MAX( 1, LDAB-1 )
136: *
137:       IF( UPPER ) THEN
138: *
139: *        Compute the Cholesky factorization A = U'*U.
140: *
141:          DO 10 J = 1, N
142: *
143: *           Compute U(J,J) and test for non-positive-definiteness.
144: *
145:             AJJ = REAL( AB( KD+1, J ) )
146:             IF( AJJ.LE.ZERO ) THEN
147:                AB( KD+1, J ) = AJJ
148:                GO TO 30
149:             END IF
150:             AJJ = SQRT( AJJ )
151:             AB( KD+1, J ) = AJJ
152: *
153: *           Compute elements J+1:J+KN of row J and update the
154: *           trailing submatrix within the band.
155: *
156:             KN = MIN( KD, N-J )
157:             IF( KN.GT.0 ) THEN
158:                CALL CSSCAL( KN, ONE / AJJ, AB( KD, J+1 ), KLD )
159:                CALL CLACGV( KN, AB( KD, J+1 ), KLD )
160:                CALL CHER( 'Upper', KN, -ONE, AB( KD, J+1 ), KLD,
161:      $                    AB( KD+1, J+1 ), KLD )
162:                CALL CLACGV( KN, AB( KD, J+1 ), KLD )
163:             END IF
164:    10    CONTINUE
165:       ELSE
166: *
167: *        Compute the Cholesky factorization A = L*L'.
168: *
169:          DO 20 J = 1, N
170: *
171: *           Compute L(J,J) and test for non-positive-definiteness.
172: *
173:             AJJ = REAL( AB( 1, J ) )
174:             IF( AJJ.LE.ZERO ) THEN
175:                AB( 1, J ) = AJJ
176:                GO TO 30
177:             END IF
178:             AJJ = SQRT( AJJ )
179:             AB( 1, J ) = AJJ
180: *
181: *           Compute elements J+1:J+KN of column J and update the
182: *           trailing submatrix within the band.
183: *
184:             KN = MIN( KD, N-J )
185:             IF( KN.GT.0 ) THEN
186:                CALL CSSCAL( KN, ONE / AJJ, AB( 2, J ), 1 )
187:                CALL CHER( 'Lower', KN, -ONE, AB( 2, J ), 1,
188:      $                    AB( 1, J+1 ), KLD )
189:             END IF
190:    20    CONTINUE
191:       END IF
192:       RETURN
193: *
194:    30 CONTINUE
195:       INFO = J
196:       RETURN
197: *
198: *     End of CPBTF2
199: *
200:       END
201: