Solution to a Linear System in a Symmetric Positive Definite Tridiagonal Matrix (Expert Driver)

Calling Sequence

CALL DPTSVX	(FACT, N, NRHS, DDIAG, DSUB, DDIAGF, DSUBF, DB, LDB, DX, LDX, DRCOND, DFERR, DBERR, DWORK, INFO)
CALL SPTSVX	(FACT, N, NRHS, SDIAG, SSUB, SDIAGF, SSUBF, SB, LDB, SX, LDX, SRCOND, SFERR, SBERR, SWORK, INFO)
CALL ZPTSVX	(FACT, N, NRHS, DDIAG, ZSUB, DDIAGF, ZSUBF, ZB, LDB, ZX, LDX, DRCOND, DFERR, DBERR, ZWORK, DWORK2, INFO)
CALL CPTSVX	(FACT, N, NRHS, SDIAG, CSUB, SDIAGF, CSUBF, CB, LDB, CX, LDX, SRCOND, SFERR, SBERR, CWORK, SWORK2, INFO)
void dptsvx	(char fact, int n, int nrhs, double *ddiag, double *dsub, double *ddiagf, double *dsubf, double *db, int ldb, double *dx, int ldx, double *drcond, double *dferr, double *dberr, int *info)
void sptsvx	(char fact, int n, int nrhs, float *sdiag, float *ssub, float *sdiagf, float *ssubf, float *sb, int ldb, float *sx, int ldx, float *srcond, float *sferr, float *sberr, int *info)
void zptsvx	(char fact, int n, int nrhs, double *ddiag, doublecomplex *zsub, double *ddiagf, doublecomplex *zsubf, doublecomplex *zb, int ldb, doublecomplex *zx, int ldx, double *drcond, double *dferr, double *dberr, int *info)
void cptsvx	(char fact, int n, int nrhs, float *sdiag, complex *csub, float *sdiagf, complex *csubf, complex *cb, int ldb, complex *cx, int ldx, float *srcond, float *sferr, float *sberr, int *info)

Arguments

FACT	Indicates whether or not the factored form of A has been supplied on entry.
	'F' or 'f'	On entry, DIAGF and SUBF contain the factored form of A. DIAG, SUB, DIAGF, and SUBF will not be modified.
	'N' or 'n'	The matrix A will be copied to DIAGF and SUBF and factored.
N	Order of the matrix A. N 0.
NRHS	Number of right-hand sides, equal to the number of columns of the matrices B and X. NRHS 0.
xDIAG	The N diagonal elements of the matrix A.
xSUB	The N-1 subdiagonal elements of the matrix A.
xDIAGF	On entry, if FACT = 'F' or 'f', then DIAGF contains the N diagonal elements of the diagonal matrix D from the LDL factorization of A as computed by xPTTRF. On exit, if FACT = 'N' or 'n', then DIAGF contains the N diagonal elements of the diagonal matrix D from the LDL factorization of A as computed by xPTTRF.
xSUBF	On entry, if FACT = 'F' or 'f', then SUBF contains the N-1 subdiagonal elements of the unit bidiagonal factor L from the LDL factorization of A, as computed by xPTTRF. On exit, if FACT = 'N' or 'n', then SUBF contains the N-1 subdiagonal elements of the unit bidiagonal factor L from the LDL factorization of A, as computed by xPTTRF.
xB	The N×NRHS right-hand side matrix B.
LDB	Leading dimension of the array B as specified in a dimension or type statement. LDB max(1, N).
xX	On exit, the N×NRHS solution matrix X.
LDX	Leading dimension of the array X as specified in a dimension or type statement. LDX max(1, N).
xRCOND	On exit, an estimate of the reciprocal condition number of the matrix A, where the reciprocal condition number is defined to be 1 / (||A|| × ||A-1||). The reciprocal of the condition number is estimated instead of the condition number itself to avoid overflow or division by zero. If RCOND is less than machine precision (in particular, if RCOND = 0) then A is singular to working precision. In this case, INFO > 0 is returned and the solution and error bounds are not computed.
xFERR	On exit, the estimated forward error bound for each solution vector X(*, j) for 1 j NRHS. If X' is the true solution corresponding to X(*, j) then FERR(j) is an upper bound on the magnitude of the largest element in X(*, j) - X' divided by the magnitude of the largest element in X(*, j).
xBERR	On exit, BERR(j) is the smallest relative change in any element of A or B(*, j) that makes X(*, j) an exact solution to AX(*, j) = B(*, j) for 1 j NRHS.
xWORK	Scratch array with a dimension of 3 × N for real subroutines or 2 × N for complex subroutines.
xWORK2	Scratch array with a dimension of N for complex subroutines.
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO < 0	The ith argument, where i = |INFO|, had an illegal value.
	1 INFO N	The leading minor of order i of A, where i = INFO, is not positive definite. The factorization has not been completed, unless i = N, and the solution and error bounds could not be computed.
	INFO = N + 1	RCOND is less than machine precision. The factorization has been completed, but the matrix is singular to working precision, and the solution and error bounds could not be computed.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER LDB, LDWORK, LDX, N, NRHS

PARAMETER (N = 4)

PARAMETER (NRHS = 1)

PARAMETER (LDB = N)

PARAMETER (LDX = N)

PARAMETER (LDWORK = 2 * N)

C

DOUBLE PRECISION B(LDB,NRHS), BERR(NRHS), DIAG(N), DIAGF(N)

DOUBLE PRECISION DLOW(N-1), DLOWF(N-1), FERR(NRHS), RCOND

DOUBLE PRECISION WORK(LDWORK), X(LDX)

INTEGER ICOL, INFO, IROW

C

EXTERNAL DPTSVX

C

C Initialize the arrays DIAG and DLOWER to store in symmetric

C tridiagonal form the 4x4 symmetric positive definite

C coefficient matrix A shown below. Initialize the array B

C to store the right hand side vector b shown below.

C

C 2 -1 0 0 6

C A = -1 2 -1 0 b = 12

C 0 -1 2 -1 12

C 0 0 -1 2 6

C

DATA DIAG / 2.0D0, 2.0D0, 2.0D0, 2.0D0 /

DATA DLOW / -1.0D0, -1.0D0, -1.0D0 /

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

C

C Print the initial values of the arrays.

C

PRINT 1000

DO 100, IROW = 1, N

PRINT 1010, (0.0D0, ICOL = 1, IROW - 2),

$ (DLOW(ICOL + 1), ICOL = ABS(IROW - 2), IROW - 2),

$ DIAG(IROW),

$ (DLOW(IROW), ICOL = 1, MIN(1, N - IROW)),

$ (0.0D0, ICOL = IROW + 2, N)

100 CONTINUE

PRINT 1020

PRINT 1030, B

C

C Solve the system and print the results.

C

CALL DPTSVX ('NOT FACTORED A', N, NRHS, DIAG, DLOW, DIAGF,

$ DLOWF, B, LDB, X, LDX, RCOND, FERR, BERR,

$ WORK, INFO)

IF (INFO .NE. 0) THEN

PRINT 1040, INFO

IF (INFO .EQ. (N + 1)) THEN

PRINT 1050

STOP 1

END IF

STOP 2

END IF

PRINT 1060

PRINT 1030, X

PRINT 1070, 1.0D0 / RCOND

PRINT 1080, (IROW, BERR(IROW), IROW = 1, NRHS)

PRINT 1090, (IROW, FERR(IROW), IROW = 1, NRHS)

C

1000 FORMAT (1X, 'A:')

1010 FORMAT (4(3X, F6.3))

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

1030 FORMAT (1X, F6.2)

1040 FORMAT (1X, 'Error solving Ax=b, INFO = ', I5)

1050 FORMAT (1X, 'Matrix is singular to working precision.')

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

1070 FORMAT (/1X, 'Estimated condition number of A: ', F7.2)

1080 FORMAT (/1X, 'Backward error in system #', I1, ': ', E12.6)

1090 FORMAT (1X, 'Forward error in system #', I1, ': ', E12.6)

C

END

Sample Output

A:
2.000 -1.000 0.000 0.000
-1.000 2.000 -1.000 0.000
0.000 -1.000 2.000 -1.000
0.000 0.000 -1.000 2.000
b:
6.00
12.00
12.00
6.00
x:
18.00
30.00
30.00
18.00
Estimated condition number of A: 12.00
Backward error in system #1: 0.296059E-16
Forward error in system #1: 0.568434E-14