Return to zgtsv.f CVS log

Up to [local] / rpl / lapack / lapack

Commit initial.

    1:       SUBROUTINE ZGTSV( N, NRHS, DL, D, DU, B, LDB, 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, LDB, N, NRHS
   10: *     ..
   11: *     .. Array Arguments ..
   12:       COMPLEX*16         B( LDB, * ), D( * ), DL( * ), DU( * )
   13: *     ..
   14: *
   15: *  Purpose
   16: *  =======
   17: *
   18: *  ZGTSV  solves the equation
   19: *
   20: *     A*X = B,
   21: *
   22: *  where A is an N-by-N tridiagonal matrix, by Gaussian elimination with
   23: *  partial pivoting.
   24: *
   25: *  Note that the equation  A'*X = B  may be solved by interchanging the
   26: *  order of the arguments DU and DL.
   27: *
   28: *  Arguments
   29: *  =========
   30: *
   31: *  N       (input) INTEGER
   32: *          The order of the matrix A.  N >= 0.
   33: *
   34: *  NRHS    (input) INTEGER
   35: *          The number of right hand sides, i.e., the number of columns
   36: *          of the matrix B.  NRHS >= 0.
   37: *
   38: *  DL      (input/output) COMPLEX*16 array, dimension (N-1)
   39: *          On entry, DL must contain the (n-1) subdiagonal elements of
   40: *          A.
   41: *          On exit, DL is overwritten by the (n-2) elements of the
   42: *          second superdiagonal of the upper triangular matrix U from
   43: *          the LU factorization of A, in DL(1), ..., DL(n-2).
   44: *
   45: *  D       (input/output) COMPLEX*16 array, dimension (N)
   46: *          On entry, D must contain the diagonal elements of A.
   47: *          On exit, D is overwritten by the n diagonal elements of U.
   48: *
   49: *  DU      (input/output) COMPLEX*16 array, dimension (N-1)
   50: *          On entry, DU must contain the (n-1) superdiagonal elements
   51: *          of A.
   52: *          On exit, DU is overwritten by the (n-1) elements of the first
   53: *          superdiagonal of U.
   54: *
   55: *  B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
   56: *          On entry, the N-by-NRHS right hand side matrix B.
   57: *          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
   58: *
   59: *  LDB     (input) INTEGER
   60: *          The leading dimension of the array B.  LDB >= max(1,N).
   61: *
   62: *  INFO    (output) INTEGER
   63: *          = 0:  successful exit
   64: *          < 0:  if INFO = -i, the i-th argument had an illegal value
   65: *          > 0:  if INFO = i, U(i,i) is exactly zero, and the solution
   66: *                has not been computed.  The factorization has not been
   67: *                completed unless i = N.
   68: *
   69: *  =====================================================================
   70: *
   71: *     .. Parameters ..
   72:       COMPLEX*16         ZERO
   73:       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
   74: *     ..
   75: *     .. Local Scalars ..
   76:       INTEGER            J, K
   77:       COMPLEX*16         MULT, TEMP, ZDUM
   78: *     ..
   79: *     .. Intrinsic Functions ..
   80:       INTRINSIC          ABS, DBLE, DIMAG, MAX
   81: *     ..
   82: *     .. External Subroutines ..
   83:       EXTERNAL           XERBLA
   84: *     ..
   85: *     .. Statement Functions ..
   86:       DOUBLE PRECISION   CABS1
   87: *     ..
   88: *     .. Statement Function definitions ..
   89:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
   90: *     ..
   91: *     .. Executable Statements ..
   92: *
   93:       INFO = 0
   94:       IF( N.LT.0 ) THEN
   95:          INFO = -1
   96:       ELSE IF( NRHS.LT.0 ) THEN
   97:          INFO = -2
   98:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
   99:          INFO = -7
  100:       END IF
  101:       IF( INFO.NE.0 ) THEN
  102:          CALL XERBLA( 'ZGTSV ', -INFO )
  103:          RETURN
  104:       END IF
  105: *
  106:       IF( N.EQ.0 )
  107:      $   RETURN
  108: *
  109:       DO 30 K = 1, N - 1
  110:          IF( DL( K ).EQ.ZERO ) THEN
  111: *
  112: *           Subdiagonal is zero, no elimination is required.
  113: *
  114:             IF( D( K ).EQ.ZERO ) THEN
  115: *
  116: *              Diagonal is zero: set INFO = K and return; a unique
  117: *              solution can not be found.
  118: *
  119:                INFO = K
  120:                RETURN
  121:             END IF
  122:          ELSE IF( CABS1( D( K ) ).GE.CABS1( DL( K ) ) ) THEN
  123: *
  124: *           No row interchange required
  125: *
  126:             MULT = DL( K ) / D( K )
  127:             D( K+1 ) = D( K+1 ) - MULT*DU( K )
  128:             DO 10 J = 1, NRHS
  129:                B( K+1, J ) = B( K+1, J ) - MULT*B( K, J )
  130:    10       CONTINUE
  131:             IF( K.LT.( N-1 ) )
  132:      $         DL( K ) = ZERO
  133:          ELSE
  134: *
  135: *           Interchange rows K and K+1
  136: *
  137:             MULT = D( K ) / DL( K )
  138:             D( K ) = DL( K )
  139:             TEMP = D( K+1 )
  140:             D( K+1 ) = DU( K ) - MULT*TEMP
  141:             IF( K.LT.( N-1 ) ) THEN
  142:                DL( K ) = DU( K+1 )
  143:                DU( K+1 ) = -MULT*DL( K )
  144:             END IF
  145:             DU( K ) = TEMP
  146:             DO 20 J = 1, NRHS
  147:                TEMP = B( K, J )
  148:                B( K, J ) = B( K+1, J )
  149:                B( K+1, J ) = TEMP - MULT*B( K+1, J )
  150:    20       CONTINUE
  151:          END IF
  152:    30 CONTINUE
  153:       IF( D( N ).EQ.ZERO ) THEN
  154:          INFO = N
  155:          RETURN
  156:       END IF
  157: *
  158: *     Back solve with the matrix U from the factorization.
  159: *
  160:       DO 50 J = 1, NRHS
  161:          B( N, J ) = B( N, J ) / D( N )
  162:          IF( N.GT.1 )
  163:      $      B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / D( N-1 )
  164:          DO 40 K = N - 2, 1, -1
  165:             B( K, J ) = ( B( K, J )-DU( K )*B( K+1, J )-DL( K )*
  166:      $                  B( K+2, J ) ) / D( K )
  167:    40    CONTINUE
  168:    50 CONTINUE
  169: *
  170:       RETURN
  171: *
  172: *     End of ZGTSV
  173: *
  174:       END