Return to zgttrf.f CVS log

Up to [local] / rpl / lapack / lapack

Initial revision

    1:       SUBROUTINE ZGTTRF( 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:       COMPLEX*16         D( * ), DL( * ), DU( * ), DU2( * )
   14: *     ..
   15: *
   16: *  Purpose
   17: *  =======
   18: *
   19: *  ZGTTRF computes an LU factorization of a complex 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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:       COMPLEX*16         FACT, TEMP, ZDUM
   81: *     ..
   82: *     .. External Subroutines ..
   83:       EXTERNAL           XERBLA
   84: *     ..
   85: *     .. Intrinsic Functions ..
   86:       INTRINSIC          ABS, DBLE, DIMAG
   87: *     ..
   88: *     .. Statement Functions ..
   89:       DOUBLE PRECISION   CABS1
   90: *     ..
   91: *     .. Statement Function definitions ..
   92:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
   93: *     ..
   94: *     .. Executable Statements ..
   95: *
   96:       INFO = 0
   97:       IF( N.LT.0 ) THEN
   98:          INFO = -1
   99:          CALL XERBLA( 'ZGTTRF', -INFO )
  100:          RETURN
  101:       END IF
  102: *
  103: *     Quick return if possible
  104: *
  105:       IF( N.EQ.0 )
  106:      $   RETURN
  107: *
  108: *     Initialize IPIV(i) = i and DU2(i) = 0
  109: *
  110:       DO 10 I = 1, N
  111:          IPIV( I ) = I
  112:    10 CONTINUE
  113:       DO 20 I = 1, N - 2
  114:          DU2( I ) = ZERO
  115:    20 CONTINUE
  116: *
  117:       DO 30 I = 1, N - 2
  118:          IF( CABS1( D( I ) ).GE.CABS1( DL( I ) ) ) THEN
  119: *
  120: *           No row interchange required, eliminate DL(I)
  121: *
  122:             IF( CABS1( D( I ) ).NE.ZERO ) THEN
  123:                FACT = DL( I ) / D( I )
  124:                DL( I ) = FACT
  125:                D( I+1 ) = D( I+1 ) - FACT*DU( I )
  126:             END IF
  127:          ELSE
  128: *
  129: *           Interchange rows I and I+1, eliminate DL(I)
  130: *
  131:             FACT = D( I ) / DL( I )
  132:             D( I ) = DL( I )
  133:             DL( I ) = FACT
  134:             TEMP = DU( I )
  135:             DU( I ) = D( I+1 )
  136:             D( I+1 ) = TEMP - FACT*D( I+1 )
  137:             DU2( I ) = DU( I+1 )
  138:             DU( I+1 ) = -FACT*DU( I+1 )
  139:             IPIV( I ) = I + 1
  140:          END IF
  141:    30 CONTINUE
  142:       IF( N.GT.1 ) THEN
  143:          I = N - 1
  144:          IF( CABS1( D( I ) ).GE.CABS1( DL( I ) ) ) THEN
  145:             IF( CABS1( D( I ) ).NE.ZERO ) THEN
  146:                FACT = DL( I ) / D( I )
  147:                DL( I ) = FACT
  148:                D( I+1 ) = D( I+1 ) - FACT*DU( I )
  149:             END IF
  150:          ELSE
  151:             FACT = D( I ) / DL( I )
  152:             D( I ) = DL( I )
  153:             DL( I ) = FACT
  154:             TEMP = DU( I )
  155:             DU( I ) = D( I+1 )
  156:             D( I+1 ) = TEMP - FACT*D( I+1 )
  157:             IPIV( I ) = I + 1
  158:          END IF
  159:       END IF
  160: *
  161: *     Check for a zero on the diagonal of U.
  162: *
  163:       DO 40 I = 1, N
  164:          IF( CABS1( D( I ) ).EQ.ZERO ) THEN
  165:             INFO = I
  166:             GO TO 50
  167:          END IF
  168:    40 CONTINUE
  169:    50 CONTINUE
  170: *
  171:       RETURN
  172: *
  173: *     End of ZGTTRF
  174: *
  175:       END