File:  [local] / rpl / lapack / lapack / dgbtf2.f
Revision 1.18: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:38:47 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    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: *> \ingroup doubleGBcomputational
  119: *
  120: *> \par Further Details:
  121: *  =====================
  122: *>
  123: *> \verbatim
  124: *>
  125: *>  The band storage scheme is illustrated by the following example, when
  126: *>  M = N = 6, KL = 2, KU = 1:
  127: *>
  128: *>  On entry:                       On exit:
  129: *>
  130: *>      *    *    *    +    +    +       *    *    *   u14  u25  u36
  131: *>      *    *    +    +    +    +       *    *   u13  u24  u35  u46
  132: *>      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
  133: *>     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
  134: *>     a21  a32  a43  a54  a65   *      m21  m32  m43  m54  m65   *
  135: *>     a31  a42  a53  a64   *    *      m31  m42  m53  m64   *    *
  136: *>
  137: *>  Array elements marked * are not used by the routine; elements marked
  138: *>  + need not be set on entry, but are required by the routine to store
  139: *>  elements of U, because of fill-in resulting from the row
  140: *>  interchanges.
  141: *> \endverbatim
  142: *>
  143: *  =====================================================================
  144:       SUBROUTINE DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO )
  145: *
  146: *  -- LAPACK computational routine --
  147: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  148: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  149: *
  150: *     .. Scalar Arguments ..
  151:       INTEGER            INFO, KL, KU, LDAB, M, N
  152: *     ..
  153: *     .. Array Arguments ..
  154:       INTEGER            IPIV( * )
  155:       DOUBLE PRECISION   AB( LDAB, * )
  156: *     ..
  157: *
  158: *  =====================================================================
  159: *
  160: *     .. Parameters ..
  161:       DOUBLE PRECISION   ONE, ZERO
  162:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  163: *     ..
  164: *     .. Local Scalars ..
  165:       INTEGER            I, J, JP, JU, KM, KV
  166: *     ..
  167: *     .. External Functions ..
  168:       INTEGER            IDAMAX
  169:       EXTERNAL           IDAMAX
  170: *     ..
  171: *     .. External Subroutines ..
  172:       EXTERNAL           DGER, DSCAL, DSWAP, XERBLA
  173: *     ..
  174: *     .. Intrinsic Functions ..
  175:       INTRINSIC          MAX, MIN
  176: *     ..
  177: *     .. Executable Statements ..
  178: *
  179: *     KV is the number of superdiagonals in the factor U, allowing for
  180: *     fill-in.
  181: *
  182:       KV = KU + KL
  183: *
  184: *     Test the input parameters.
  185: *
  186:       INFO = 0
  187:       IF( M.LT.0 ) THEN
  188:          INFO = -1
  189:       ELSE IF( N.LT.0 ) THEN
  190:          INFO = -2
  191:       ELSE IF( KL.LT.0 ) THEN
  192:          INFO = -3
  193:       ELSE IF( KU.LT.0 ) THEN
  194:          INFO = -4
  195:       ELSE IF( LDAB.LT.KL+KV+1 ) THEN
  196:          INFO = -6
  197:       END IF
  198:       IF( INFO.NE.0 ) THEN
  199:          CALL XERBLA( 'DGBTF2', -INFO )
  200:          RETURN
  201:       END IF
  202: *
  203: *     Quick return if possible
  204: *
  205:       IF( M.EQ.0 .OR. N.EQ.0 )
  206:      $   RETURN
  207: *
  208: *     Gaussian elimination with partial pivoting
  209: *
  210: *     Set fill-in elements in columns KU+2 to KV to zero.
  211: *
  212:       DO 20 J = KU + 2, MIN( KV, N )
  213:          DO 10 I = KV - J + 2, KL
  214:             AB( I, J ) = ZERO
  215:    10    CONTINUE
  216:    20 CONTINUE
  217: *
  218: *     JU is the index of the last column affected by the current stage
  219: *     of the factorization.
  220: *
  221:       JU = 1
  222: *
  223:       DO 40 J = 1, MIN( M, N )
  224: *
  225: *        Set fill-in elements in column J+KV to zero.
  226: *
  227:          IF( J+KV.LE.N ) THEN
  228:             DO 30 I = 1, KL
  229:                AB( I, J+KV ) = ZERO
  230:    30       CONTINUE
  231:          END IF
  232: *
  233: *        Find pivot and test for singularity. KM is the number of
  234: *        subdiagonal elements in the current column.
  235: *
  236:          KM = MIN( KL, M-J )
  237:          JP = IDAMAX( KM+1, AB( KV+1, J ), 1 )
  238:          IPIV( J ) = JP + J - 1
  239:          IF( AB( KV+JP, J ).NE.ZERO ) THEN
  240:             JU = MAX( JU, MIN( J+KU+JP-1, N ) )
  241: *
  242: *           Apply interchange to columns J to JU.
  243: *
  244:             IF( JP.NE.1 )
  245:      $         CALL DSWAP( JU-J+1, AB( KV+JP, J ), LDAB-1,
  246:      $                     AB( KV+1, J ), LDAB-1 )
  247: *
  248:             IF( KM.GT.0 ) THEN
  249: *
  250: *              Compute multipliers.
  251: *
  252:                CALL DSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 )
  253: *
  254: *              Update trailing submatrix within the band.
  255: *
  256:                IF( JU.GT.J )
  257:      $            CALL DGER( KM, JU-J, -ONE, AB( KV+2, J ), 1,
  258:      $                       AB( KV, J+1 ), LDAB-1, AB( KV+1, J+1 ),
  259:      $                       LDAB-1 )
  260:             END IF
  261:          ELSE
  262: *
  263: *           If pivot is zero, set INFO to the index of the pivot
  264: *           unless a zero pivot has already been found.
  265: *
  266:             IF( INFO.EQ.0 )
  267:      $         INFO = J
  268:          END IF
  269:    40 CONTINUE
  270:       RETURN
  271: *
  272: *     End of DGBTF2
  273: *
  274:       END

CVSweb interface <joel.bertrand@systella.fr>