1*> \brief \b DGBTF2 computes the LU factorization of a general band matrix using the unblocked version of the algorithm.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGBTF2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgbtf2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgbtf2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgbtf2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO )
22*
23*       .. Scalar Arguments ..
24*       INTEGER            INFO, KL, KU, LDAB, M, N
25*       ..
26*       .. Array Arguments ..
27*       INTEGER            IPIV( * )
28*       DOUBLE PRECISION   AB( LDAB, * )
29*       ..
30*
31*
32*> \par Purpose:
33*  =============
34*>
35*> \verbatim
36*>
37*> DGBTF2 computes an LU factorization of a real m-by-n band matrix A
38*> using partial pivoting with row interchanges.
39*>
40*> This is the unblocked version of the algorithm, calling Level 2 BLAS.
41*> \endverbatim
42*
43*  Arguments:
44*  ==========
45*
46*> \param[in] M
47*> \verbatim
48*>          M is INTEGER
49*>          The number of rows of the matrix A.  M >= 0.
50*> \endverbatim
51*>
52*> \param[in] N
53*> \verbatim
54*>          N is INTEGER
55*>          The number of columns of the matrix A.  N >= 0.
56*> \endverbatim
57*>
58*> \param[in] KL
59*> \verbatim
60*>          KL is INTEGER
61*>          The number of subdiagonals within the band of A.  KL >= 0.
62*> \endverbatim
63*>
64*> \param[in] KU
65*> \verbatim
66*>          KU is INTEGER
67*>          The number of superdiagonals within the band of A.  KU >= 0.
68*> \endverbatim
69*>
70*> \param[in,out] AB
71*> \verbatim
72*>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
73*>          On entry, the matrix A in band storage, in rows KL+1 to
74*>          2*KL+KU+1; rows 1 to KL of the array need not be set.
75*>          The j-th column of A is stored in the j-th column of the
76*>          array AB as follows:
77*>          AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
78*>
79*>          On exit, details of the factorization: U is stored as an
80*>          upper triangular band matrix with KL+KU superdiagonals in
81*>          rows 1 to KL+KU+1, and the multipliers used during the
82*>          factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
83*>          See below for further details.
84*> \endverbatim
85*>
86*> \param[in] LDAB
87*> \verbatim
88*>          LDAB is INTEGER
89*>          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
90*> \endverbatim
91*>
92*> \param[out] IPIV
93*> \verbatim
94*>          IPIV is INTEGER array, dimension (min(M,N))
95*>          The pivot indices; for 1 <= i <= min(M,N), row i of the
96*>          matrix was interchanged with row IPIV(i).
97*> \endverbatim
98*>
99*> \param[out] INFO
100*> \verbatim
101*>          INFO is INTEGER
102*>          = 0: successful exit
103*>          < 0: if INFO = -i, the i-th argument had an illegal value
104*>          > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
105*>               has been completed, but the factor U is exactly
106*>               singular, and division by zero will occur if it is used
107*>               to solve a system of equations.
108*> \endverbatim
109*
110*  Authors:
111*  ========
112*
113*> \author Univ. of Tennessee
114*> \author Univ. of California Berkeley
115*> \author Univ. of Colorado Denver
116*> \author NAG Ltd.
117*
118*> \date September 2012
119*
120*> \ingroup doubleGBcomputational
121*
122*> \par Further Details:
123*  =====================
124*>
125*> \verbatim
126*>
127*>  The band storage scheme is illustrated by the following example, when
128*>  M = N = 6, KL = 2, KU = 1:
129*>
130*>  On entry:                       On exit:
131*>
132*>      *    *    *    +    +    +       *    *    *   u14  u25  u36
133*>      *    *    +    +    +    +       *    *   u13  u24  u35  u46
134*>      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
135*>     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
136*>     a21  a32  a43  a54  a65   *      m21  m32  m43  m54  m65   *
137*>     a31  a42  a53  a64   *    *      m31  m42  m53  m64   *    *
138*>
139*>  Array elements marked * are not used by the routine; elements marked
140*>  + need not be set on entry, but are required by the routine to store
141*>  elements of U, because of fill-in resulting from the row
142*>  interchanges.
143*> \endverbatim
144*>
145*  =====================================================================
146      SUBROUTINE DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO )
147*
148*  -- LAPACK computational routine (version 3.4.2) --
149*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
150*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
151*     September 2012
152*
153*     .. Scalar Arguments ..
154      INTEGER            INFO, KL, KU, LDAB, M, N
155*     ..
156*     .. Array Arguments ..
157      INTEGER            IPIV( * )
158      DOUBLE PRECISION   AB( LDAB, * )
159*     ..
160*
161*  =====================================================================
162*
163*     .. Parameters ..
164      DOUBLE PRECISION   ONE, ZERO
165      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
166*     ..
167*     .. Local Scalars ..
168      INTEGER            I, J, JP, JU, KM, KV
169*     ..
170*     .. External Functions ..
171      INTEGER            IDAMAX
172      EXTERNAL           IDAMAX
173*     ..
174*     .. External Subroutines ..
175      EXTERNAL           DGER, DSCAL, DSWAP, XERBLA
176*     ..
177*     .. Intrinsic Functions ..
178      INTRINSIC          MAX, MIN
179*     ..
180*     .. Executable Statements ..
181*
182*     KV is the number of superdiagonals in the factor U, allowing for
183*     fill-in.
184*
185      KV = KU + KL
186*
187*     Test the input parameters.
188*
189      INFO = 0
190      IF( M.LT.0 ) THEN
191         INFO = -1
192      ELSE IF( N.LT.0 ) THEN
193         INFO = -2
194      ELSE IF( KL.LT.0 ) THEN
195         INFO = -3
196      ELSE IF( KU.LT.0 ) THEN
197         INFO = -4
198      ELSE IF( LDAB.LT.KL+KV+1 ) THEN
199         INFO = -6
200      END IF
201      IF( INFO.NE.0 ) THEN
202         CALL XERBLA( 'DGBTF2', -INFO )
203         RETURN
204      END IF
205*
206*     Quick return if possible
207*
208      IF( M.EQ.0 .OR. N.EQ.0 )
209     $   RETURN
210*
211*     Gaussian elimination with partial pivoting
212*
213*     Set fill-in elements in columns KU+2 to KV to zero.
214*
215      DO 20 J = KU + 2, MIN( KV, N )
216         DO 10 I = KV - J + 2, KL
217            AB( I, J ) = ZERO
218   10    CONTINUE
219   20 CONTINUE
220*
221*     JU is the index of the last column affected by the current stage
222*     of the factorization.
223*
224      JU = 1
225*
226      DO 40 J = 1, MIN( M, N )
227*
228*        Set fill-in elements in column J+KV to zero.
229*
230         IF( J+KV.LE.N ) THEN
231            DO 30 I = 1, KL
232               AB( I, J+KV ) = ZERO
233   30       CONTINUE
234         END IF
235*
236*        Find pivot and test for singularity. KM is the number of
237*        subdiagonal elements in the current column.
238*
239         KM = MIN( KL, M-J )
240         JP = IDAMAX( KM+1, AB( KV+1, J ), 1 )
241         IPIV( J ) = JP + J - 1
242         IF( AB( KV+JP, J ).NE.ZERO ) THEN
243            JU = MAX( JU, MIN( J+KU+JP-1, N ) )
244*
245*           Apply interchange to columns J to JU.
246*
247            IF( JP.NE.1 )
248     $         CALL DSWAP( JU-J+1, AB( KV+JP, J ), LDAB-1,
249     $                     AB( KV+1, J ), LDAB-1 )
250*
251            IF( KM.GT.0 ) THEN
252*
253*              Compute multipliers.
254*
255               CALL DSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 )
256*
257*              Update trailing submatrix within the band.
258*
259               IF( JU.GT.J )
260     $            CALL DGER( KM, JU-J, -ONE, AB( KV+2, J ), 1,
261     $                       AB( KV, J+1 ), LDAB-1, AB( KV+1, J+1 ),
262     $                       LDAB-1 )
263            END IF
264         ELSE
265*
266*           If pivot is zero, set INFO to the index of the pivot
267*           unless a zero pivot has already been found.
268*
269            IF( INFO.EQ.0 )
270     $         INFO = J
271         END IF
272   40 CONTINUE
273      RETURN
274*
275*     End of DGBTF2
276*
277      END
278