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

Calling Sequence

CALL DGESL	(DA, LDA, N, IPIVOT, DB, JOB)
CALL SGESL	(SA, LDA, N, IPIVOT, SB, JOB)
CALL ZGESL	(ZA, LDA, N, IPIVOT, ZB, JOB)
CALL CGESL	(CA, LDA, N, IPIVOT, CB, JOB)
void dgesl	(double *da, int lda, int n, int *ipivot, double *db, int job)
void sgesl	(float *sa, int lda, int n, int *ipivot, float *sb, int job)
void zgesl	(doublecomplex *za, int lda, int n, int *ipivot, doublecomplex *zb, int job)
void cgesl	(complex *ca, int lda, int n, int *ipivot, complex *cb, int job)

Arguments

xA	LU factorization of the matrix A, as computed by xGECO or xGEFA.
LDA	Leading dimension of the array A as specified in a dimension or type statement. LDA max(1,N).
N	Order of the matrix A. N 0.
IPIVOT	Pivot vector as computed by xGECO or xGEFA.
xB	On entry, the right-hand side vector b. On exit, the solution vector x.
JOB	Determines which operation the subroutine will perform:
	0	solve the system Ax = b.
	not 0	solve the linear system AHx = b
		Note that ATx = AHx for real matrices.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER IAXEQB, LDA, LDB, N

PARAMETER (IAXEQB = 0)

PARAMETER (N = 3)

PARAMETER (LDA = N)

PARAMETER (LDB = LDA)

C

DOUBLE PRECISION A(LDA,N), B(LDB)

INTEGER ICOL, INFO, IPIVOT(N), IROW, JOB

C

EXTERNAL DGEFA, DGESL

C

C Initialize the array A to store the matrix A shown below.

C Initialize the array B to store the vector b shown below.

C

C 1 2 2 15

C A = 2 1 2 b = 15

C 2 2 1 15

C

DATA A / 1.0D0, 3*2.0D0, 1.0D0, 3*2.0D0, 1.0D0 /

DATA B / 3*1.5D1 /

C

PRINT 1000

PRINT 1010, ((A(IROW,ICOL), ICOL = 1, N), IROW = 1, N)

PRINT 1020

PRINT 1030, B

CALL DGEFA (A, LDA, N, IPIVOT, INFO)

IF (INFO .EQ. 0) THEN

JOB = IAXEQB

CALL DGESL (A, LDA, N, IPIVOT, B, JOB)

PRINT 1040

PRINT 1030, B

ELSE

PRINT 1050, INFO

END IF

C

1000 FORMAT (1X, 'A:')

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

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

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

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

1050 FORMAT (1X, 'A appears singular at ', I2)

C

END

Sample Output

A:
1.0 2.0 2.0
2.0 1.0 2.0
2.0 2.0 1.0
b:
15.0
15.0
15.0
A**(-1)*b
3.0
3.0
3.0