001:       SUBROUTINE ZGBTF2( M, N, KL, KU, AB, LDAB, IPIV, 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:       INTEGER            INFO, KL, KU, LDAB, M, N
009: *     ..
010: *     .. Array Arguments ..
011:       INTEGER            IPIV( * )
012:       COMPLEX*16         AB( LDAB, * )
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  ZGBTF2 computes an LU factorization of a complex m-by-n band matrix
019: *  A using partial pivoting with row interchanges.
020: *
021: *  This is the unblocked version of the algorithm, calling Level 2 BLAS.
022: *
023: *  Arguments
024: *  =========
025: *
026: *  M       (input) INTEGER
027: *          The number of rows of the matrix A.  M >= 0.
028: *
029: *  N       (input) INTEGER
030: *          The number of columns of the matrix A.  N >= 0.
031: *
032: *  KL      (input) INTEGER
033: *          The number of subdiagonals within the band of A.  KL >= 0.
034: *
035: *  KU      (input) INTEGER
036: *          The number of superdiagonals within the band of A.  KU >= 0.
037: *
038: *  AB      (input/output) COMPLEX*16 array, dimension (LDAB,N)
039: *          On entry, the matrix A in band storage, in rows KL+1 to
040: *          2*KL+KU+1; rows 1 to KL of the array need not be set.
041: *          The j-th column of A is stored in the j-th column of the
042: *          array AB as follows:
043: *          AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
044: *
045: *          On exit, details of the factorization: U is stored as an
046: *          upper triangular band matrix with KL+KU superdiagonals in
047: *          rows 1 to KL+KU+1, and the multipliers used during the
048: *          factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
049: *          See below for further details.
050: *
051: *  LDAB    (input) INTEGER
052: *          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
053: *
054: *  IPIV    (output) INTEGER array, dimension (min(M,N))
055: *          The pivot indices; for 1 <= i <= min(M,N), row i of the
056: *          matrix was interchanged with row IPIV(i).
057: *
058: *  INFO    (output) INTEGER
059: *          = 0: successful exit
060: *          < 0: if INFO = -i, the i-th argument had an illegal value
061: *          > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
062: *               has been completed, but the factor U is exactly
063: *               singular, and division by zero will occur if it is used
064: *               to solve a system of equations.
065: *
066: *  Further Details
067: *  ===============
068: *
069: *  The band storage scheme is illustrated by the following example, when
070: *  M = N = 6, KL = 2, KU = 1:
071: *
072: *  On entry:                       On exit:
073: *
074: *      *    *    *    +    +    +       *    *    *   u14  u25  u36
075: *      *    *    +    +    +    +       *    *   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: *     a21  a32  a43  a54  a65   *      m21  m32  m43  m54  m65   *
079: *     a31  a42  a53  a64   *    *      m31  m42  m53  m64   *    *
080: *
081: *  Array elements marked * are not used by the routine; elements marked
082: *  + need not be set on entry, but are required by the routine to store
083: *  elements of U, because of fill-in resulting from the row
084: *  interchanges.
085: *
086: *  =====================================================================
087: *
088: *     .. Parameters ..
089:       COMPLEX*16         ONE, ZERO
090:       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
091:      $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
092: *     ..
093: *     .. Local Scalars ..
094:       INTEGER            I, J, JP, JU, KM, KV
095: *     ..
096: *     .. External Functions ..
097:       INTEGER            IZAMAX
098:       EXTERNAL           IZAMAX
099: *     ..
100: *     .. External Subroutines ..
101:       EXTERNAL           XERBLA, ZGERU, ZSCAL, ZSWAP
102: *     ..
103: *     .. Intrinsic Functions ..
104:       INTRINSIC          MAX, MIN
105: *     ..
106: *     .. Executable Statements ..
107: *
108: *     KV is the number of superdiagonals in the factor U, allowing for
109: *     fill-in.
110: *
111:       KV = KU + KL
112: *
113: *     Test the input parameters.
114: *
115:       INFO = 0
116:       IF( M.LT.0 ) THEN
117:          INFO = -1
118:       ELSE IF( N.LT.0 ) THEN
119:          INFO = -2
120:       ELSE IF( KL.LT.0 ) THEN
121:          INFO = -3
122:       ELSE IF( KU.LT.0 ) THEN
123:          INFO = -4
124:       ELSE IF( LDAB.LT.KL+KV+1 ) THEN
125:          INFO = -6
126:       END IF
127:       IF( INFO.NE.0 ) THEN
128:          CALL XERBLA( 'ZGBTF2', -INFO )
129:          RETURN
130:       END IF
131: *
132: *     Quick return if possible
133: *
134:       IF( M.EQ.0 .OR. N.EQ.0 )
135:      $   RETURN
136: *
137: *     Gaussian elimination with partial pivoting
138: *
139: *     Set fill-in elements in columns KU+2 to KV to zero.
140: *
141:       DO 20 J = KU + 2, MIN( KV, N )
142:          DO 10 I = KV - J + 2, KL
143:             AB( I, J ) = ZERO
144:    10    CONTINUE
145:    20 CONTINUE
146: *
147: *     JU is the index of the last column affected by the current stage
148: *     of the factorization.
149: *
150:       JU = 1
151: *
152:       DO 40 J = 1, MIN( M, N )
153: *
154: *        Set fill-in elements in column J+KV to zero.
155: *
156:          IF( J+KV.LE.N ) THEN
157:             DO 30 I = 1, KL
158:                AB( I, J+KV ) = ZERO
159:    30       CONTINUE
160:          END IF
161: *
162: *        Find pivot and test for singularity. KM is the number of
163: *        subdiagonal elements in the current column.
164: *
165:          KM = MIN( KL, M-J )
166:          JP = IZAMAX( KM+1, AB( KV+1, J ), 1 )
167:          IPIV( J ) = JP + J - 1
168:          IF( AB( KV+JP, J ).NE.ZERO ) THEN
169:             JU = MAX( JU, MIN( J+KU+JP-1, N ) )
170: *
171: *           Apply interchange to columns J to JU.
172: *
173:             IF( JP.NE.1 )
174:      $         CALL ZSWAP( JU-J+1, AB( KV+JP, J ), LDAB-1,
175:      $                     AB( KV+1, J ), LDAB-1 )
176:             IF( KM.GT.0 ) THEN
177: *
178: *              Compute multipliers.
179: *
180:                CALL ZSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 )
181: *
182: *              Update trailing submatrix within the band.
183: *
184:                IF( JU.GT.J )
185:      $            CALL ZGERU( KM, JU-J, -ONE, AB( KV+2, J ), 1,
186:      $                        AB( KV, J+1 ), LDAB-1, AB( KV+1, J+1 ),
187:      $                        LDAB-1 )
188:             END IF
189:          ELSE
190: *
191: *           If pivot is zero, set INFO to the index of the pivot
192: *           unless a zero pivot has already been found.
193: *
194:             IF( INFO.EQ.0 )
195:      $         INFO = J
196:          END IF
197:    40 CONTINUE
198:       RETURN
199: *
200: *     End of ZGBTF2
201: *
202:       END
203: