Previous Next Contents Generated Index Doc Set Home



Solution to a Linear System in an LU-Factored General Tridiagonal Matrix

The subroutines described in this section solve a linear system
AX = B, ATX = B, or AHX = B for a tridiagonal matrix A, which has been LU-factored by xGTTRF, and general matrices B and X.

Calling Sequence

CALL DGTTRS 
(TRANSA, N, NRHS, DLOW, DDIAG, DUP1, DUP2, IPIVOT, DB, 
LDB, INFO)
CALL SGTTRS 
(TRANSA, N, NRHS, SLOW, SDIAG, SUP1, SUP2, IPIVOT, SB, 
LDB, INFO)
CALL ZGTTRS 
(TRANSA, N, NRHS, ZLOW, ZDIAG, ZUP1, ZUP2, IPIVOT, ZB, 
LDB, INFO)
CALL CGTTRS 
(TRANSA, N, NRHS, CLOW, CDIAG, CUP1, CUP2, IPIVOT, CB, 
LDB, INFO)






void dgttrs 
(char trans, int n, int nrhs, double *dlow, double 
*ddiag, double *dup1, double *dup2, int *ipivot, double 
*db, int ldb, int *info)
void sgttrs 
(char trans, int n, int nrhs, float *dlow, float 
*ddiag, float *dup1, float *dup2, int *ipivot, float 
*sb, int ldb, int *info)
void zgttrs 
(char trans, int n, int nrhs, doublecomplex *zlow, 
doublecomplex *zdiag, doublecomplex *zup1, 
doublecomplex *zup2, int *ipivot, doublecomplex *zb, 
int ldb, int *info)
void cgttrs 
(char trans, int n, int nrhs, complex *clow, complex 
*cdiag, complex *cup1, complex *cup2, int *ipivot, 
complex *cb, int ldb, int *info)

Arguments

TRANSA

Indicates the form of the linear system to solve. The legal values for TRANSA are listed below. Any values not listed below are illegal.

'N' or 'n'

No transpose, solve AX = B.

'T' or 't'

Transpose, solve ATX = B.

'C' or 'c'

Conjugate transpose, solve AHX = B.

Note that AT = AH for real matrices.

N

Order of the matrix A. N 0.

NRHS

Number of right-hand sides, equal to the number of columns of the matrix B. NRHS 0.

xLOW

The N-1 multipliers that define the matrix L from the LU factorization of A, as computed by xGTTRF.

xDIAG

The N diagonal elements of the upper triangular matrix U from the LU factorization of A, as computed by xGTTRF.

xUP1

The N-1 elements of the first superdiagonal of U from the LU factorization of A, as computed by xGTTRF.

xUP2

The N-2 elements of the second superdiagonal of U from the LU factorization of A, as computed by xGTTRF.

IPIVOT

Pivot indices as computed by xGTTRF.

xB

On entry, the N×NRHS right-hand side matrix B.
On exit, the N×NRHS solution matrix X.

LDB

Leading dimension of the array B as specified in a dimension or type statement. LDB max(1, N).

INFO

On exit:

INFO = 0

Subroutine completed normally.

INFO < 0

The ith argument, where i = |INFO|, had an illegal value.

Sample Program




      PROGRAM TEST
      IMPLICIT NONE
C
      INTEGER           LDB, LDIWRK, LDWORK, N, NRHS
      PARAMETER        (N = 5)
      PARAMETER        (LDB = N)
      PARAMETER        (LDIWRK = N)
      PARAMETER        (LDWORK = 2 * N)
      PARAMETER        (NRHS = 1)
C
      DOUBLE PRECISION  ANORM, B(LDB,NRHS), DIAG(N), DLOWER(N-1)
      DOUBLE PRECISION  DUP1(N-1), DUP2(N-2), RCOND, WORK(LDWORK)
      INTEGER           ICOL, INFO, IPIVOT(N), IROW, IWORK(LDIWRK)
C
      EXTERNAL          DGTCON, DGTTRF, DGTTRS
      INTRINSIC         ABS, MAX
C
C     Initialize the arrays DLOWER, DIAG, and DUP1 to store 
C     the first subdiagonal, the diagonal, and the first
C     superdiagonal of the coefficient matrix A shown below. 
C     Initialize the array B to store the right hand side 
C     matrix b shown below.
C
C         -3   1                      6
C          1  -2   1                 12
C     A =      1  -2   1         b = 18
C                  1  -2   1         12
C                      1  -1          6
C
      DATA DLOWER / 1.0D0, 1.0D0, 1.0D0, 1.0D0 /
      DATA DIAG / -3.0D0, -2.0D0, -2.0D0, -2.0D0, -1.0D0 /
      DATA DUP1 / 1.0D0, 1.0D0, 1.0D0, 1.0D0 /
      DATA B / 6.0D0, 1.2D1, 1.8D1, 1.2D1, 6.0D0 /
C
C     Print the initial value of A.
C
      PRINT 1000
      DO 100, IROW = 1, N
        PRINT 1010, (0.0D0, ICOL = 1, IROW - 2),
     $              (DLOWER(ICOL + 1), ICOL = ABS(IROW - 2),
     $              IROW - 2), DIAG(IROW),
     $              (DUP1(IROW), ICOL = 1, MIN(1, N - IROW)),
     $              (0.0D0, ICOL = IROW + 2, N)
  100 CONTINUE
      PRINT 1020
      PRINT 1030, B
C
C     LU factor A.
C
      CALL DGTTRF (N, DLOWER, DIAG, DUP1, DUP2, IPIVOT, INFO)
      IF (INFO .NE. 0) THEN
        PRINT 1040, INFO
        STOP 1
      END IF
      ANORM = 4.0D0
C
C     Estimate and print the condition number of A.
C
      CALL DGTCON ('ONE-NORM', N, DLOWER, DIAG, DUP1, DUP2,
     $             IPIVOT, ANORM, RCOND, WORK, IWORK, INFO)
      IF (INFO .NE. 0) THEN
        PRINT 1050, ABS (INFO)
        STOP 2
      END IF
      PRINT 1060, 1.0D0 / RCOND
C
C     Solve the system and print the solution.
C
      CALL DGTTRS ('NO TRANSPOSE A', N, NRHS, DLOWER, DIAG, DUP1,
     $             DUP2, IPIVOT, B, LDB, INFO)
      IF (INFO .NE. 0) THEN
        PRINT 1070, INFO
        STOP 3
      END IF
      PRINT 1080
      PRINT 1030, B
C
 1000 FORMAT (1X, 'A:')
 1010 FORMAT (5(3X, F5.2))
 1020 FORMAT (/1X, 'b:')
 1030 FORMAT (1X, F8.2)
 1040 FORMAT (1X, 'Error factoring A, INFO = ', I5)
 1050 FORMAT (1X, 'Illegal argument to DGTCON, argument #', I2)
 1060 FORMAT (/1X, 'Estimated condition number of A: ', F6.2)
 1070 FORMAT (1X, 'Error solving Ax=b, INFO = ', I5)
 1080 FORMAT (/1X, 'x:')
C
      END
 

Sample Output

 
 A:
   -3.00    1.00    0.00    0.00    0.00
    1.00   -2.00    1.00    0.00    0.00
    0.00    1.00   -2.00    1.00    0.00
    0.00    0.00    1.00   -2.00    1.00
    0.00    0.00    0.00    1.00   -1.00



 b:
     6.00
    12.00
    18.00
    12.00
     6.00



 Estimated condition number of A:  50.00



 x:
   -27.00
   -75.00
  -111.00
  -129.00
  -135.00






Previous Next Contents Generated Index Doc Set Home