Solution to a Linear System in a Tridiagonal Matrix

Calling Sequence

CALL DGTSL	(N, DSUB, DDIAG, DSUPER, DB, INFO)
CALL SGTSL	(N, SSUB, SDIAG, SSUPER, SB, INFO)
CALL ZGTSL	(N, ZSUB, ZDIAG, ZSUPER, ZB, INFO)
CALL CGTSL	(N, CSUB, CDIAG, CSUPER, CB, INFO)
void dgtsl	(int n, double *dsub, double *ddiag, double *dsuper, double *db, int *info)
void sgtsl	(int n, float *ssub, float *sdiag, float *ssuper, float *sb, int *info)
void zgtsl	(int n, doublecomplex *zsub, doublecomplex *zdiag, doublecomplex *zsuper, doublecomplex *zb, int *info)
void cgtsl	(int n, complex *csub, complex *cdiag, complex *csuper, complex *cb, int *info)

Arguments

N	Order of the matrix A. N 0.
xSUB	Subdiagonal of A. SUB(2) through SUB(N) contain the subdiagonal and SUB(1) is not referenced.
xDIAG	Diagonal of A.
xSUPER	Superdiagonal of A. SUPER(1) through SUPER(N-1) contain the superdiagonal and SUPER(N) is not referenced.
xB	On entry, the right-hand side vector b. On exit, the solution vector x.
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO > 0	Returns a value k when the kth element of the diagonal becomes exactly zero.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER N

PARAMETER (N = 5)

C

DOUBLE PRECISION B(N), SUBD(N), DIAG(N), SUPERD(N)

INTEGER INFO

C

EXTERNAL DGTSL

C

C Initialize the arrays SUBD, DIAG, and SUPERD to store

C the subdiagonal, diagonal, and superdiagonal of the

C tridiagonal matrix A shown below. Initialize the array

C B to store the right hand side of Ax=b.

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 SUBD / 8D8, 4*1.0D0 /

DATA DIAG / -3.0D0, -2.0D0, -2.0D0, -2.0D0, -1.0D0 /

DATA SUPERD / 4*1.0D0, 8D8 /

DATA B / 6.0D0, 1.2D1, 1.8D1, 1.2D1, 6.0D0 /

C

PRINT 1000

PRINT 1010, DIAG(1), SUPERD(1)

PRINT 1020, SUBD(2), DIAG(2), SUPERD(2)

PRINT 1030, SUBD(3), DIAG(3), SUPERD(3)

PRINT 1040, SUBD(4), DIAG(4), SUPERD(4)

PRINT 1050, SUBD(5), DIAG(5)

PRINT 1060

PRINT 1070, B

CALL DGTSL (N, SUBD, DIAG, SUPERD, B, INFO)

IF (INFO .EQ. 0) THEN

PRINT 1080

PRINT 1070, B

ELSE

PRINT 1090, INFO

END IF

C

1000 FORMAT (1X, 'A:')

1010 FORMAT (1X, 2(2X, F4.1))

1020 FORMAT (1X, 3(2X, F4.1))

1030 FORMAT (1X, 6X, 3(2X, F4.1))

1040 FORMAT (1X, 12X, 3(2X, F4.1))

1050 FORMAT (1X, 18X, 2(2X, F4.1))

1060 FORMAT (/1X, 'b:')

1070 FORMAT (3X, F6.1)

1080 FORMAT (/1X, 'A**(-1) * b:')

1090 FORMAT (/1X, 'A appears singular at ', I1)

C

END

Sample Output

A:
-3.0 1.0
1.0 -2.0 1.0
1.0 -2.0 1.0
1.0 -2.0 1.0
1.0 -1.0
b:
6.0
12.0
18.0
12.0
6.0
A**(-1) * b:
-27.0
-75.0
-111.0
-129.0
-135.0