UDU Factorization and Condition Number of a Symmetric Matrix

Calling Sequence

CALL DSICO	(DA, LDA, N, IPIVOT, DRCOND, DWORK)
CALL SSICO	(SA, LDA, N, IPIVOT, SRCOND, SWORK)
CALL ZSICO	(ZA, LDA, N, IPIVOT, DRCOND, ZWORK)
CALL CSICO	(CA, LDA, N, IPIVOT, SRCOND, CWORK)
void dsico	(double *da, int lda, int n, int *ipivot, double *drcond)
void ssico	(float *sa, int lda, int n, int *ipivot, float *srcond)
void zsico	(doublecomplex *za, int lda, int n, int *ipivot, double *drcond)
void csico	(complex *ca, int lda, int n, int *ipivot, float *srcond)

Arguments

xA	On entry, the upper triangle of the matrix A. On exit, a UDU factorization of the matrix A. The strict lower triangle of A is not referenced.
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 LDA, N

PARAMETER (N = 4)

PARAMETER (LDA = 5)

C

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

INTEGER ICOL, IPIVOT(N), IROW

C

EXTERNAL DSICO, DSISL

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 -.5 -.5 -.5 -.5 12

C A = -.5 -1.5 -1.5 -1.5 b = 6

C -.5 -1.5 -2.5 -2.5 6

C -.5 -1.5 -2.5 -3.5 12

C

DATA A / -5.0D-1, 4*8D8, -5.0D-1, -1.5D0, 3*8D8, -5.0D-1,

$ -1.5D0, -2.5D0, 2*8D8, -5.0D-1, -1.5D0, -2.5D0,

$ -3.5D0, 8D8 /

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

C

PRINT 1000

DO 100, IROW = 1, N

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

$ (A(IROW,ICOL), ICOL = IROW + 1, N)

100 CONTINUE

PRINT 1020

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

PRINT 1030

PRINT 1040, B

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

PRINT 1050, RCOND

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

PRINT 1060

END IF

PRINT 1070, 1.0D0 / RCOND

CALL DSISL (A, LDA, N, IPIVOT, B)

PRINT 1080

PRINT 1040, B

C

1000 FORMAT (1X, 'A in full form:')

1010 FORMAT (4(3X, F5.1))

1020 FORMAT (/1X, 'A in symmetric form: (* in unused elements)')

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

1040 FORMAT (3X, F5.1)

1050 FORMAT (/1X, 'Reciprocal condition number of A:', F6.3)

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

1070 FORMAT (1X, 'Condition number of A: ', F6.3)

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

C

END

Sample Output

A in full form:
-0.5 -0.5 -0.5 -0.5
-0.5 -1.5 -1.5 -1.5
-0.5 -1.5 -2.5 -2.5
-0.5 -1.5 -2.5 -3.5
A in symmetric form: (* in unused elements)
-0.5 -0.5 -0.5 -0.5
***** -1.5 -1.5 -1.5
***** ***** -2.5 -2.5
***** ***** ***** -3.5
b:
12.0
6.0
6.0
12.0
Reciprocal condition number of A: 0.031
Condition number of A: 32.000
A**(-1) * b:
-30.0
6.0
6.0
-6.0