Return to dgttrf.f CVS log

Up to [local] / rpl / lapack / lapack

Initial revision

    1:       SUBROUTINE DGTTRF( N, DL, D, DU, DU2, 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, N
   10: *     ..
   11: *     .. Array Arguments ..
   12:       INTEGER            IPIV( * )
   13:       DOUBLE PRECISION   D( * ), DL( * ), DU( * ), DU2( * )
   14: *     ..
   15: *
   16: *  Purpose
   17: *  =======
   18: *
   19: *  DGTTRF computes an LU factorization of a real tridiagonal matrix A
   20: *  using elimination with partial pivoting and row interchanges.
   21: *
   22: *  The factorization has the form
   23: *     A = L * U
   24: *  where L is a product of permutation and unit lower bidiagonal
   25: *  matrices and U is upper triangular with nonzeros in only the main
   26: *  diagonal and first two superdiagonals.
   27: *
   28: *  Arguments
   29: *  =========
   30: *
   31: *  N       (input) INTEGER
   32: *          The order of the matrix A.
   33: *
   34: *  DL      (input/output) DOUBLE PRECISION array, dimension (N-1)
   35: *          On entry, DL must contain the (n-1) sub-diagonal elements of
   36: *          A.
   37: *
   38: *          On exit, DL is overwritten by the (n-1) multipliers that
   39: *          define the matrix L from the LU factorization of A.
   40: *
   41: *  D       (input/output) DOUBLE PRECISION array, dimension (N)
   42: *          On entry, D must contain the diagonal elements of A.
   43: *
   44: *          On exit, D is overwritten by the n diagonal elements of the
   45: *          upper triangular matrix U from the LU factorization of A.
   46: *
   47: *  DU      (input/output) DOUBLE PRECISION array, dimension (N-1)
   48: *          On entry, DU must contain the (n-1) super-diagonal elements
   49: *          of A.
   50: *
   51: *          On exit, DU is overwritten by the (n-1) elements of the first
   52: *          super-diagonal of U.
   53: *
   54: *  DU2     (output) DOUBLE PRECISION array, dimension (N-2)
   55: *          On exit, DU2 is overwritten by the (n-2) elements of the
   56: *          second super-diagonal of U.
   57: *
   58: *  IPIV    (output) INTEGER array, dimension (N)
   59: *          The pivot indices; for 1 <= i <= n, row i of the matrix was
   60: *          interchanged with row IPIV(i).  IPIV(i) will always be either
   61: *          i or i+1; IPIV(i) = i indicates a row interchange was not
   62: *          required.
   63: *
   64: *  INFO    (output) INTEGER
   65: *          = 0:  successful exit
   66: *          < 0:  if INFO = -k, the k-th argument had an illegal value
   67: *          > 0:  if INFO = k, U(k,k) is exactly zero. The factorization
   68: *                has been completed, but the factor U is exactly
   69: *                singular, and division by zero will occur if it is used
   70: *                to solve a system of equations.
   71: *
   72: *  =====================================================================
   73: *
   74: *     .. Parameters ..
   75:       DOUBLE PRECISION   ZERO
   76:       PARAMETER          ( ZERO = 0.0D+0 )
   77: *     ..
   78: *     .. Local Scalars ..
   79:       INTEGER            I
   80:       DOUBLE PRECISION   FACT, TEMP
   81: *     ..
   82: *     .. Intrinsic Functions ..
   83:       INTRINSIC          ABS
   84: *     ..
   85: *     .. External Subroutines ..
   86:       EXTERNAL           XERBLA
   87: *     ..
   88: *     .. Executable Statements ..
   89: *
   90:       INFO = 0
   91:       IF( N.LT.0 ) THEN
   92:          INFO = -1
   93:          CALL XERBLA( 'DGTTRF', -INFO )
   94:          RETURN
   95:       END IF
   96: *
   97: *     Quick return if possible
   98: *
   99:       IF( N.EQ.0 )
  100:      $   RETURN
  101: *
  102: *     Initialize IPIV(i) = i and DU2(I) = 0
  103: *
  104:       DO 10 I = 1, N
  105:          IPIV( I ) = I
  106:    10 CONTINUE
  107:       DO 20 I = 1, N - 2
  108:          DU2( I ) = ZERO
  109:    20 CONTINUE
  110: *
  111:       DO 30 I = 1, N - 2
  112:          IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
  113: *
  114: *           No row interchange required, eliminate DL(I)
  115: *
  116:             IF( D( I ).NE.ZERO ) THEN
  117:                FACT = DL( I ) / D( I )
  118:                DL( I ) = FACT
  119:                D( I+1 ) = D( I+1 ) - FACT*DU( I )
  120:             END IF
  121:          ELSE
  122: *
  123: *           Interchange rows I and I+1, eliminate DL(I)
  124: *
  125:             FACT = D( I ) / DL( I )
  126:             D( I ) = DL( I )
  127:             DL( I ) = FACT
  128:             TEMP = DU( I )
  129:             DU( I ) = D( I+1 )
  130:             D( I+1 ) = TEMP - FACT*D( I+1 )
  131:             DU2( I ) = DU( I+1 )
  132:             DU( I+1 ) = -FACT*DU( I+1 )
  133:             IPIV( I ) = I + 1
  134:          END IF
  135:    30 CONTINUE
  136:       IF( N.GT.1 ) THEN
  137:          I = N - 1
  138:          IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
  139:             IF( D( I ).NE.ZERO ) THEN
  140:                FACT = DL( I ) / D( I )
  141:                DL( I ) = FACT
  142:                D( I+1 ) = D( I+1 ) - FACT*DU( I )
  143:             END IF
  144:          ELSE
  145:             FACT = D( I ) / DL( I )
  146:             D( I ) = DL( I )
  147:             DL( I ) = FACT
  148:             TEMP = DU( I )
  149:             DU( I ) = D( I+1 )
  150:             D( I+1 ) = TEMP - FACT*D( I+1 )
  151:             IPIV( I ) = I + 1
  152:          END IF
  153:       END IF
  154: *
  155: *     Check for a zero on the diagonal of U.
  156: *
  157:       DO 40 I = 1, N
  158:          IF( D( I ).EQ.ZERO ) THEN
  159:             INFO = I
  160:             GO TO 50
  161:          END IF
  162:    40 CONTINUE
  163:    50 CONTINUE
  164: *
  165:       RETURN
  166: *
  167: *     End of DGTTRF
  168: *
  169:       END