File:  [local] / rpl / lapack / lapack / dgttrf.f
Revision 1.18: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:38:52 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 DGTTRF
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download DGTTRF + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgttrf.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgttrf.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgttrf.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DGTTRF( N, DL, D, DU, DU2, IPIV, INFO )
   22: *
   23: *       .. Scalar Arguments ..
   24: *       INTEGER            INFO, N
   25: *       ..
   26: *       .. Array Arguments ..
   27: *       INTEGER            IPIV( * )
   28: *       DOUBLE PRECISION   D( * ), DL( * ), DU( * ), DU2( * )
   29: *       ..
   30: *
   31: *
   32: *> \par Purpose:
   33: *  =============
   34: *>
   35: *> \verbatim
   36: *>
   37: *> DGTTRF computes an LU factorization of a real tridiagonal matrix A
   38: *> using elimination with partial pivoting and row interchanges.
   39: *>
   40: *> The factorization has the form
   41: *>    A = L * U
   42: *> where L is a product of permutation and unit lower bidiagonal
   43: *> matrices and U is upper triangular with nonzeros in only the main
   44: *> diagonal and first two superdiagonals.
   45: *> \endverbatim
   46: *
   47: *  Arguments:
   48: *  ==========
   49: *
   50: *> \param[in] N
   51: *> \verbatim
   52: *>          N is INTEGER
   53: *>          The order of the matrix A.
   54: *> \endverbatim
   55: *>
   56: *> \param[in,out] DL
   57: *> \verbatim
   58: *>          DL is DOUBLE PRECISION array, dimension (N-1)
   59: *>          On entry, DL must contain the (n-1) sub-diagonal elements of
   60: *>          A.
   61: *>
   62: *>          On exit, DL is overwritten by the (n-1) multipliers that
   63: *>          define the matrix L from the LU factorization of A.
   64: *> \endverbatim
   65: *>
   66: *> \param[in,out] D
   67: *> \verbatim
   68: *>          D is DOUBLE PRECISION array, dimension (N)
   69: *>          On entry, D must contain the diagonal elements of A.
   70: *>
   71: *>          On exit, D is overwritten by the n diagonal elements of the
   72: *>          upper triangular matrix U from the LU factorization of A.
   73: *> \endverbatim
   74: *>
   75: *> \param[in,out] DU
   76: *> \verbatim
   77: *>          DU is DOUBLE PRECISION array, dimension (N-1)
   78: *>          On entry, DU must contain the (n-1) super-diagonal elements
   79: *>          of A.
   80: *>
   81: *>          On exit, DU is overwritten by the (n-1) elements of the first
   82: *>          super-diagonal of U.
   83: *> \endverbatim
   84: *>
   85: *> \param[out] DU2
   86: *> \verbatim
   87: *>          DU2 is DOUBLE PRECISION array, dimension (N-2)
   88: *>          On exit, DU2 is overwritten by the (n-2) elements of the
   89: *>          second super-diagonal of U.
   90: *> \endverbatim
   91: *>
   92: *> \param[out] IPIV
   93: *> \verbatim
   94: *>          IPIV is INTEGER array, dimension (N)
   95: *>          The pivot indices; for 1 <= i <= n, row i of the matrix was
   96: *>          interchanged with row IPIV(i).  IPIV(i) will always be either
   97: *>          i or i+1; IPIV(i) = i indicates a row interchange was not
   98: *>          required.
   99: *> \endverbatim
  100: *>
  101: *> \param[out] INFO
  102: *> \verbatim
  103: *>          INFO is INTEGER
  104: *>          = 0:  successful exit
  105: *>          < 0:  if INFO = -k, the k-th argument had an illegal value
  106: *>          > 0:  if INFO = k, U(k,k) is exactly zero. The factorization
  107: *>                has been completed, but the factor U is exactly
  108: *>                singular, and division by zero will occur if it is used
  109: *>                to solve a system of equations.
  110: *> \endverbatim
  111: *
  112: *  Authors:
  113: *  ========
  114: *
  115: *> \author Univ. of Tennessee
  116: *> \author Univ. of California Berkeley
  117: *> \author Univ. of Colorado Denver
  118: *> \author NAG Ltd.
  119: *
  120: *> \ingroup doubleGTcomputational
  121: *
  122: *  =====================================================================
  123:       SUBROUTINE DGTTRF( N, DL, D, DU, DU2, IPIV, INFO )
  124: *
  125: *  -- LAPACK computational routine --
  126: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  127: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  128: *
  129: *     .. Scalar Arguments ..
  130:       INTEGER            INFO, N
  131: *     ..
  132: *     .. Array Arguments ..
  133:       INTEGER            IPIV( * )
  134:       DOUBLE PRECISION   D( * ), DL( * ), DU( * ), DU2( * )
  135: *     ..
  136: *
  137: *  =====================================================================
  138: *
  139: *     .. Parameters ..
  140:       DOUBLE PRECISION   ZERO
  141:       PARAMETER          ( ZERO = 0.0D+0 )
  142: *     ..
  143: *     .. Local Scalars ..
  144:       INTEGER            I
  145:       DOUBLE PRECISION   FACT, TEMP
  146: *     ..
  147: *     .. Intrinsic Functions ..
  148:       INTRINSIC          ABS
  149: *     ..
  150: *     .. External Subroutines ..
  151:       EXTERNAL           XERBLA
  152: *     ..
  153: *     .. Executable Statements ..
  154: *
  155:       INFO = 0
  156:       IF( N.LT.0 ) THEN
  157:          INFO = -1
  158:          CALL XERBLA( 'DGTTRF', -INFO )
  159:          RETURN
  160:       END IF
  161: *
  162: *     Quick return if possible
  163: *
  164:       IF( N.EQ.0 )
  165:      $   RETURN
  166: *
  167: *     Initialize IPIV(i) = i and DU2(I) = 0
  168: *
  169:       DO 10 I = 1, N
  170:          IPIV( I ) = I
  171:    10 CONTINUE
  172:       DO 20 I = 1, N - 2
  173:          DU2( I ) = ZERO
  174:    20 CONTINUE
  175: *
  176:       DO 30 I = 1, N - 2
  177:          IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
  178: *
  179: *           No row interchange required, eliminate DL(I)
  180: *
  181:             IF( D( I ).NE.ZERO ) THEN
  182:                FACT = DL( I ) / D( I )
  183:                DL( I ) = FACT
  184:                D( I+1 ) = D( I+1 ) - FACT*DU( I )
  185:             END IF
  186:          ELSE
  187: *
  188: *           Interchange rows I and I+1, eliminate DL(I)
  189: *
  190:             FACT = D( I ) / DL( I )
  191:             D( I ) = DL( I )
  192:             DL( I ) = FACT
  193:             TEMP = DU( I )
  194:             DU( I ) = D( I+1 )
  195:             D( I+1 ) = TEMP - FACT*D( I+1 )
  196:             DU2( I ) = DU( I+1 )
  197:             DU( I+1 ) = -FACT*DU( I+1 )
  198:             IPIV( I ) = I + 1
  199:          END IF
  200:    30 CONTINUE
  201:       IF( N.GT.1 ) THEN
  202:          I = N - 1
  203:          IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
  204:             IF( D( I ).NE.ZERO ) THEN
  205:                FACT = DL( I ) / D( I )
  206:                DL( I ) = FACT
  207:                D( I+1 ) = D( I+1 ) - FACT*DU( I )
  208:             END IF
  209:          ELSE
  210:             FACT = D( I ) / DL( I )
  211:             D( I ) = DL( I )
  212:             DL( I ) = FACT
  213:             TEMP = DU( I )
  214:             DU( I ) = D( I+1 )
  215:             D( I+1 ) = TEMP - FACT*D( I+1 )
  216:             IPIV( I ) = I + 1
  217:          END IF
  218:       END IF
  219: *
  220: *     Check for a zero on the diagonal of U.
  221: *
  222:       DO 40 I = 1, N
  223:          IF( D( I ).EQ.ZERO ) THEN
  224:             INFO = I
  225:             GO TO 50
  226:          END IF
  227:    40 CONTINUE
  228:    50 CONTINUE
  229: *
  230:       RETURN
  231: *
  232: *     End of DGTTRF
  233: *
  234:       END

CVSweb interface <joel.bertrand@systella.fr>