001:       SUBROUTINE DPBTF2( 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:       DOUBLE PRECISION   AB( LDAB, * )
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  DPBTF2 computes the Cholesky factorization of a real symmetric
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 transpose of U, and
025: *  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: *          symmetric 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) DOUBLE PRECISION array, dimension (LDAB,N)
046: *          On entry, the upper or lower triangle of the symmetric 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:       DOUBLE PRECISION   ONE, ZERO
093:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
094: *     ..
095: *     .. Local Scalars ..
096:       LOGICAL            UPPER
097:       INTEGER            J, KLD, KN
098:       DOUBLE PRECISION   AJJ
099: *     ..
100: *     .. External Functions ..
101:       LOGICAL            LSAME
102:       EXTERNAL           LSAME
103: *     ..
104: *     .. External Subroutines ..
105:       EXTERNAL           DSCAL, DSYR, XERBLA
106: *     ..
107: *     .. Intrinsic Functions ..
108:       INTRINSIC          MAX, MIN, 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( 'DPBTF2', -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 = AB( KD+1, J )
146:             IF( AJJ.LE.ZERO )
147:      $         GO TO 30
148:             AJJ = SQRT( AJJ )
149:             AB( KD+1, J ) = AJJ
150: *
151: *           Compute elements J+1:J+KN of row J and update the
152: *           trailing submatrix within the band.
153: *
154:             KN = MIN( KD, N-J )
155:             IF( KN.GT.0 ) THEN
156:                CALL DSCAL( KN, ONE / AJJ, AB( KD, J+1 ), KLD )
157:                CALL DSYR( 'Upper', KN, -ONE, AB( KD, J+1 ), KLD,
158:      $                    AB( KD+1, J+1 ), KLD )
159:             END IF
160:    10    CONTINUE
161:       ELSE
162: *
163: *        Compute the Cholesky factorization A = L*L'.
164: *
165:          DO 20 J = 1, N
166: *
167: *           Compute L(J,J) and test for non-positive-definiteness.
168: *
169:             AJJ = AB( 1, J )
170:             IF( AJJ.LE.ZERO )
171:      $         GO TO 30
172:             AJJ = SQRT( AJJ )
173:             AB( 1, J ) = AJJ
174: *
175: *           Compute elements J+1:J+KN of column J and update the
176: *           trailing submatrix within the band.
177: *
178:             KN = MIN( KD, N-J )
179:             IF( KN.GT.0 ) THEN
180:                CALL DSCAL( KN, ONE / AJJ, AB( 2, J ), 1 )
181:                CALL DSYR( 'Lower', KN, -ONE, AB( 2, J ), 1,
182:      $                    AB( 1, J+1 ), KLD )
183:             END IF
184:    20    CONTINUE
185:       END IF
186:       RETURN
187: *
188:    30 CONTINUE
189:       INFO = J
190:       RETURN
191: *
192: *     End of DPBTF2
193: *
194:       END
195: