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