Return to dgbtf2.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour globale de Lapack 3.2.2.

    1:       SUBROUTINE DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO )
    2: *
    3: *  -- LAPACK routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     .. Scalar Arguments ..
    9:       INTEGER            INFO, KL, KU, LDAB, M, N
   10: *     ..
   11: *     .. Array Arguments ..
   12:       INTEGER            IPIV( * )
   13:       DOUBLE PRECISION   AB( LDAB, * )
   14: *     ..
   15: *
   16: *  Purpose
   17: *  =======
   18: *
   19: *  DGBTF2 computes an LU factorization of a real m-by-n band matrix A
   20: *  using partial pivoting with row interchanges.
   21: *
   22: *  This is the unblocked version of the algorithm, calling Level 2 BLAS.
   23: *
   24: *  Arguments
   25: *  =========
   26: *
   27: *  M       (input) INTEGER
   28: *          The number of rows of the matrix A.  M >= 0.
   29: *
   30: *  N       (input) INTEGER
   31: *          The number of columns of the matrix A.  N >= 0.
   32: *
   33: *  KL      (input) INTEGER
   34: *          The number of subdiagonals within the band of A.  KL >= 0.
   35: *
   36: *  KU      (input) INTEGER
   37: *          The number of superdiagonals within the band of A.  KU >= 0.
   38: *
   39: *  AB      (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
   40: *          On entry, the matrix A in band storage, in rows KL+1 to
   41: *          2*KL+KU+1; rows 1 to KL of the array need not be set.
   42: *          The j-th column of A is stored in the j-th column of the
   43: *          array AB as follows:
   44: *          AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
   45: *
   46: *          On exit, details of the factorization: U is stored as an
   47: *          upper triangular band matrix with KL+KU superdiagonals in
   48: *          rows 1 to KL+KU+1, and the multipliers used during the
   49: *          factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
   50: *          See below for further details.
   51: *
   52: *  LDAB    (input) INTEGER
   53: *          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
   54: *
   55: *  IPIV    (output) INTEGER array, dimension (min(M,N))
   56: *          The pivot indices; for 1 <= i <= min(M,N), row i of the
   57: *          matrix was interchanged with row IPIV(i).
   58: *
   59: *  INFO    (output) INTEGER
   60: *          = 0: successful exit
   61: *          < 0: if INFO = -i, the i-th argument had an illegal value
   62: *          > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
   63: *               has been completed, but the factor U is exactly
   64: *               singular, and division by zero will occur if it is used
   65: *               to solve a system of equations.
   66: *
   67: *  Further Details
   68: *  ===============
   69: *
   70: *  The band storage scheme is illustrated by the following example, when
   71: *  M = N = 6, KL = 2, KU = 1:
   72: *
   73: *  On entry:                       On exit:
   74: *
   75: *      *    *    *    +    +    +       *    *    *   u14  u25  u36
   76: *      *    *    +    +    +    +       *    *   u13  u24  u35  u46
   77: *      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
   78: *     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
   79: *     a21  a32  a43  a54  a65   *      m21  m32  m43  m54  m65   *
   80: *     a31  a42  a53  a64   *    *      m31  m42  m53  m64   *    *
   81: *
   82: *  Array elements marked * are not used by the routine; elements marked
   83: *  + need not be set on entry, but are required by the routine to store
   84: *  elements of U, because of fill-in resulting from the row
   85: *  interchanges.
   86: *
   87: *  =====================================================================
   88: *
   89: *     .. Parameters ..
   90:       DOUBLE PRECISION   ONE, ZERO
   91:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
   92: *     ..
   93: *     .. Local Scalars ..
   94:       INTEGER            I, J, JP, JU, KM, KV
   95: *     ..
   96: *     .. External Functions ..
   97:       INTEGER            IDAMAX
   98:       EXTERNAL           IDAMAX
   99: *     ..
  100: *     .. External Subroutines ..
  101:       EXTERNAL           DGER, DSCAL, DSWAP, XERBLA
  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( 'DGBTF2', -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 = IDAMAX( 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 DSWAP( JU-J+1, AB( KV+JP, J ), LDAB-1,
  175:      $                     AB( KV+1, J ), LDAB-1 )
  176: *
  177:             IF( KM.GT.0 ) THEN
  178: *
  179: *              Compute multipliers.
  180: *
  181:                CALL DSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 )
  182: *
  183: *              Update trailing submatrix within the band.
  184: *
  185:                IF( JU.GT.J )
  186:      $            CALL DGER( KM, JU-J, -ONE, AB( KV+2, J ), 1,
  187:      $                       AB( KV, J+1 ), LDAB-1, AB( KV+1, J+1 ),
  188:      $                       LDAB-1 )
  189:             END IF
  190:          ELSE
  191: *
  192: *           If pivot is zero, set INFO to the index of the pivot
  193: *           unless a zero pivot has already been found.
  194: *
  195:             IF( INFO.EQ.0 )
  196:      $         INFO = J
  197:          END IF
  198:    40 CONTINUE
  199:       RETURN
  200: *
  201: *     End of DGBTF2
  202: *
  203:       END