SRC/zgbtf2.f

   1 *> \brief \b ZGBTF2 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 ZGBTF2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgbtf2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgbtf2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgbtf2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZGBTF2( 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 *       COMPLEX*16         AB( LDAB, * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> ZGBTF2 computes an LU factorization of a complex m-by-n band matrix
  38 *> A 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 COMPLEX*16 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 complex16GBcomputational
 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 ZGBTF2( 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       COMPLEX*16         AB( LDAB, * )
 159 *     ..
 160 *
 161 *  =====================================================================
 162 *
 163 *     .. Parameters ..
 164       COMPLEX*16         ONE, ZERO
 165       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
 166      $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
 167 *     ..
 168 *     .. Local Scalars ..
 169       INTEGER            I, J, JP, JU, KM, KV
 170 *     ..
 171 *     .. External Functions ..
 172       INTEGER            IZAMAX
 173       EXTERNAL           IZAMAX
 174 *     ..
 175 *     .. External Subroutines ..
 176       EXTERNAL           XERBLA, ZGERU, ZSCAL, ZSWAP
 177 *     ..
 178 *     .. Intrinsic Functions ..
 179       INTRINSIC          MAX, MIN
 180 *     ..
 181 *     .. Executable Statements ..
 182 *
 183 *     KV is the number of superdiagonals in the factor U, allowing for
 184 *     fill-in.
 185 *
 186       KV = KU + KL
 187 *
 188 *     Test the input parameters.
 189 *
 190       INFO = 0
 191       IF( M.LT.0 ) THEN
 192          INFO = -1
 193       ELSE IF( N.LT.0 ) THEN
 194          INFO = -2
 195       ELSE IF( KL.LT.0 ) THEN
 196          INFO = -3
 197       ELSE IF( KU.LT.0 ) THEN
 198          INFO = -4
 199       ELSE IF( LDAB.LT.KL+KV+1 ) THEN
 200          INFO = -6
 201       END IF
 202       IF( INFO.NE.0 ) THEN
 203          CALL XERBLA( 'ZGBTF2', -INFO )
 204          RETURN
 205       END IF
 206 *
 207 *     Quick return if possible
 208 *
 209       IF( M.EQ.0 .OR. N.EQ.0 )
 210      $   RETURN
 211 *
 212 *     Gaussian elimination with partial pivoting
 213 *
 214 *     Set fill-in elements in columns KU+2 to KV to zero.
 215 *
 216       DO 20 J = KU + 2, MIN( KV, N )
 217          DO 10 I = KV - J + 2, KL
 218             AB( I, J ) = ZERO
 219    10    CONTINUE
 220    20 CONTINUE
 221 *
 222 *     JU is the index of the last column affected by the current stage
 223 *     of the factorization.
 224 *
 225       JU = 1
 226 *
 227       DO 40 J = 1, MIN( M, N )
 228 *
 229 *        Set fill-in elements in column J+KV to zero.
 230 *
 231          IF( J+KV.LE.N ) THEN
 232             DO 30 I = 1, KL
 233                AB( I, J+KV ) = ZERO
 234    30       CONTINUE
 235          END IF
 236 *
 237 *        Find pivot and test for singularity. KM is the number of
 238 *        subdiagonal elements in the current column.
 239 *
 240          KM = MIN( KL, M-J )
 241          JP = IZAMAX( KM+1, AB( KV+1, J ), 1 )
 242          IPIV( J ) = JP + J - 1
 243          IF( AB( KV+JP, J ).NE.ZERO ) THEN
 244             JU = MAX( JU, MIN( J+KU+JP-1, N ) )
 245 *
 246 *           Apply interchange to columns J to JU.
 247 *
 248             IF( JP.NE.1 )
 249      $         CALL ZSWAP( JU-J+1, AB( KV+JP, J ), LDAB-1,
 250      $                     AB( KV+1, J ), LDAB-1 )
 251             IF( KM.GT.0 ) THEN
 252 *
 253 *              Compute multipliers.
 254 *
 255                CALL ZSCAL( 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 ZGERU( 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 ZGBTF2
 276 *
 277       END