1*> \brief \b DPBTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite band matrix (unblocked algorithm).
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DPBTF2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpbtf2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpbtf2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpbtf2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DPBTF2( UPLO, N, KD, AB, LDAB, INFO )
22*
23*       .. Scalar Arguments ..
24*       CHARACTER          UPLO
25*       INTEGER            INFO, KD, LDAB, N
26*       ..
27*       .. Array Arguments ..
28*       DOUBLE PRECISION   AB( LDAB, * )
29*       ..
30*
31*
32*> \par Purpose:
33*  =============
34*>
35*> \verbatim
36*>
37*> DPBTF2 computes the Cholesky factorization of a real symmetric
38*> positive definite band matrix A.
39*>
40*> The factorization has the form
41*>    A = U**T * U ,  if UPLO = 'U', or
42*>    A = L  * L**T,  if UPLO = 'L',
43*> where U is an upper triangular matrix, U**T is the transpose of U, and
44*> L is lower triangular.
45*>
46*> This is the unblocked version of the algorithm, calling Level 2 BLAS.
47*> \endverbatim
48*
49*  Arguments:
50*  ==========
51*
52*> \param[in] UPLO
53*> \verbatim
54*>          UPLO is CHARACTER*1
55*>          Specifies whether the upper or lower triangular part of the
56*>          symmetric matrix A is stored:
57*>          = 'U':  Upper triangular
58*>          = 'L':  Lower triangular
59*> \endverbatim
60*>
61*> \param[in] N
62*> \verbatim
63*>          N is INTEGER
64*>          The order of the matrix A.  N >= 0.
65*> \endverbatim
66*>
67*> \param[in] KD
68*> \verbatim
69*>          KD is INTEGER
70*>          The number of super-diagonals of the matrix A if UPLO = 'U',
71*>          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0.
72*> \endverbatim
73*>
74*> \param[in,out] AB
75*> \verbatim
76*>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
77*>          On entry, the upper or lower triangle of the symmetric band
78*>          matrix A, stored in the first KD+1 rows of the array.  The
79*>          j-th column of A is stored in the j-th column of the array AB
80*>          as follows:
81*>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
82*>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
83*>
84*>          On exit, if INFO = 0, the triangular factor U or L from the
85*>          Cholesky factorization A = U**T*U or A = L*L**T of the band
86*>          matrix A, in the same storage format as A.
87*> \endverbatim
88*>
89*> \param[in] LDAB
90*> \verbatim
91*>          LDAB is INTEGER
92*>          The leading dimension of the array AB.  LDAB >= KD+1.
93*> \endverbatim
94*>
95*> \param[out] INFO
96*> \verbatim
97*>          INFO is INTEGER
98*>          = 0: successful exit
99*>          < 0: if INFO = -k, the k-th argument had an illegal value
100*>          > 0: if INFO = k, the leading minor of order k is not
101*>               positive definite, and the factorization could not be
102*>               completed.
103*> \endverbatim
104*
105*  Authors:
106*  ========
107*
108*> \author Univ. of Tennessee
109*> \author Univ. of California Berkeley
110*> \author Univ. of Colorado Denver
111*> \author NAG Ltd.
112*
113*> \ingroup doubleOTHERcomputational
114*
115*> \par Further Details:
116*  =====================
117*>
118*> \verbatim
119*>
120*>  The band storage scheme is illustrated by the following example, when
121*>  N = 6, KD = 2, and UPLO = 'U':
122*>
123*>  On entry:                       On exit:
124*>
125*>      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
126*>      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
127*>     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
128*>
129*>  Similarly, if UPLO = 'L' the format of A is as follows:
130*>
131*>  On entry:                       On exit:
132*>
133*>     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
134*>     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
135*>     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *
136*>
137*>  Array elements marked * are not used by the routine.
138*> \endverbatim
139*>
140*  =====================================================================
141      SUBROUTINE DPBTF2( UPLO, N, KD, AB, LDAB, INFO )
142*
143*  -- LAPACK computational routine --
144*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
145*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
146*
147*     .. Scalar Arguments ..
148      CHARACTER          UPLO
149      INTEGER            INFO, KD, LDAB, N
150*     ..
151*     .. Array Arguments ..
152      DOUBLE PRECISION   AB( LDAB, * )
153*     ..
154*
155*  =====================================================================
156*
157*     .. Parameters ..
158      DOUBLE PRECISION   ONE, ZERO
159      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
160*     ..
161*     .. Local Scalars ..
162      LOGICAL            UPPER
163      INTEGER            J, KLD, KN
164      DOUBLE PRECISION   AJJ
165*     ..
166*     .. External Functions ..
167      LOGICAL            LSAME
168      EXTERNAL           LSAME
169*     ..
170*     .. External Subroutines ..
171      EXTERNAL           DSCAL, DSYR, XERBLA
172*     ..
173*     .. Intrinsic Functions ..
174      INTRINSIC          MAX, MIN, SQRT
175*     ..
176*     .. Executable Statements ..
177*
178*     Test the input parameters.
179*
180      INFO = 0
181      UPPER = LSAME( UPLO, 'U' )
182      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
183         INFO = -1
184      ELSE IF( N.LT.0 ) THEN
185         INFO = -2
186      ELSE IF( KD.LT.0 ) THEN
187         INFO = -3
188      ELSE IF( LDAB.LT.KD+1 ) THEN
189         INFO = -5
190      END IF
191      IF( INFO.NE.0 ) THEN
192         CALL XERBLA( 'DPBTF2', -INFO )
193         RETURN
194      END IF
195*
196*     Quick return if possible
197*
198      IF( N.EQ.0 )
199     $   RETURN
200*
201      KLD = MAX( 1, LDAB-1 )
202*
203      IF( UPPER ) THEN
204*
205*        Compute the Cholesky factorization A = U**T*U.
206*
207         DO 10 J = 1, N
208*
209*           Compute U(J,J) and test for non-positive-definiteness.
210*
211            AJJ = AB( KD+1, J )
212            IF( AJJ.LE.ZERO )
213     $         GO TO 30
214            AJJ = SQRT( AJJ )
215            AB( KD+1, J ) = AJJ
216*
217*           Compute elements J+1:J+KN of row J and update the
218*           trailing submatrix within the band.
219*
220            KN = MIN( KD, N-J )
221            IF( KN.GT.0 ) THEN
222               CALL DSCAL( KN, ONE / AJJ, AB( KD, J+1 ), KLD )
223               CALL DSYR( 'Upper', KN, -ONE, AB( KD, J+1 ), KLD,
224     $                    AB( KD+1, J+1 ), KLD )
225            END IF
226   10    CONTINUE
227      ELSE
228*
229*        Compute the Cholesky factorization A = L*L**T.
230*
231         DO 20 J = 1, N
232*
233*           Compute L(J,J) and test for non-positive-definiteness.
234*
235            AJJ = AB( 1, J )
236            IF( AJJ.LE.ZERO )
237     $         GO TO 30
238            AJJ = SQRT( AJJ )
239            AB( 1, J ) = AJJ
240*
241*           Compute elements J+1:J+KN of column J and update the
242*           trailing submatrix within the band.
243*
244            KN = MIN( KD, N-J )
245            IF( KN.GT.0 ) THEN
246               CALL DSCAL( KN, ONE / AJJ, AB( 2, J ), 1 )
247               CALL DSYR( 'Lower', KN, -ONE, AB( 2, J ), 1,
248     $                    AB( 1, J+1 ), KLD )
249            END IF
250   20    CONTINUE
251      END IF
252      RETURN
253*
254   30 CONTINUE
255      INFO = J
256      RETURN
257*
258*     End of DPBTF2
259*
260      END
261