LU Factorization and Condition Number of a General Matrix

Calling Sequence

CALL DGECO	(DA, LDA, N, IPIVOT, DRCOND, DWORK)
CALL SGECO	(SA, LDA, N, IPIVOT, SRCOND, SWORK)
CALL ZGECO	(ZA, LDA, N, IPIVOT, DRCOND, ZWORK)
CALL CGECO	(CA, LDA, N, IPIVOT, SRCOND, CWORK)
void dgeco	(double *da, int lda, int n, int *ipivot, double *drcond)
void sgeco	(float *sa, int lda, int n, int *ipivot, float *srcond)
void zgeco	(doublecomplex *za, int lda, int n, int *ipivot, double *drcond)
void cgeco	(complex *ca, int lda, int n, int *ipivot, float *srcond)

Arguments

xA	On entry, the matrix A. On exit, an LU factorization of A.
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	On exit, a vector of pivot indices.
xRCOND	On exit, an estimate of the reciprocal condition number of A. 0.0 RCOND 1.0. As the value of RCOND gets smaller, operations with A such as solving Ax = b may become less stable. If RCOND satisfies RCOND + 1.0 = 1.0 then A may be singular to working precision.
xWORK	Scratch array with a dimension of N.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER IAXEQB, LDA, LDB, N

PARAMETER (IAXEQB = 0)

PARAMETER (N = 3)

PARAMETER (LDA = N)

PARAMETER (LDB = LDA)

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

INTEGER ICOL, IPIVOT(N), IROW, JOB

C

EXTERNAL DGECO, 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 /

PRINT 1000

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

PRINT 1020

PRINT 1030, B

CALL DGECO (A, LDA, N, IPIVOT, RCOND, WORK)

PRINT 1040, RCOND

IF ((RCOND + 1.0D0) .EQ. 1.0D0) THEN

PRINT 1060

END IF

JOB = IAXEQB

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

PRINT 1050

PRINT 1030, B

C

1000 FORMAT (1X, 'A:')

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

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

1030 FORMAT (3X, F4.1)

1040 FORMAT (/1X, 'Estimated reciprocal condition number:', F7.4)

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

1060 FORMAT (1X, 'A may be singular to working precision.')

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
Estimated reciprocal condition number: 0.1429
A**(-1) * b:
3.0
3.0
3.0