UDU Factorization and Condition Number of a Hermitian Matrix

Calling Sequence

CALL ZHICO	(ZA, LDA, N, IPIVOT, DRCOND, ZWORK)
CALL CHICO	(CA, LDA, N, IPIVOT, SRCOND, CWORK)
void zhico	(doublecomplex *za, int lda, int n, int *ipivot, double *drcond)
void chico	(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 = 3)

PARAMETER (LDA = 3)

C

REAL RCOND

COMPLEX A(LDA,N), B(N), WORK(N)

INTEGER ICOL, IPIVOT(N), IROW

C

EXTERNAL CHICO, CHISL

INTRINSIC CONJG

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 1+2i 1+2i 95-180i

C A = 1+2i 6 -2+6i b = 545-118i

C 1+2i -2+6i 11 865+ 62i

C

DATA A / (1.0,0.0), (8E8,8E8), (8E8,8E8),

$ (1.0,-2.0), (6.0,0.0), (8E8,8E8),

$ (1.0,-2.0), (6.0,-2.0), (11.0,0.0) /

DATA B / (95.0,-180.0), (545.0,-118.0), (865.0,62.0) /

C

PRINT 1000

DO 100, IROW = 1, N

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

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

100 CONTINUE

PRINT 1020

DO 110, IROW = 1, N

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

110 CONTINUE

PRINT 1030

PRINT 1040, B

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

PRINT 1050, RCOND

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

PRINT 1060

END IF

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

PRINT 1070

PRINT 1040, B

C

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

1010 FORMAT (4(: 3X, '(', F4.1, ',', F4.1, ')'))

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

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

1040 FORMAT (3X, '(', F6.1, ',', F6.1, ')')

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

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

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

C

END

Sample Output

A in full form:
( 1.0, 0.0) ( 1.0,-2.0) ( 1.0,-2.0)
( 1.0, 2.0) ( 6.0, 0.0) ( 6.0,-2.0)
( 1.0, 2.0) ( 6.0, 2.0) (11.0, 0.0)
A in Hermitian form: (* in unused elements)
( 1.0, 0.0) ( 1.0,-2.0) ( 1.0,-2.0)
(****,****) ( 6.0, 0.0) ( 6.0,-2.0)
(****,****) (****,****) (11.0, 0.0)
b:
( 95.0,-180.0)
( 545.0,-118.0)
( 865.0, 62.0)
Reciprocal condition number of A: 0.001
A**(-1) * b:
( 5.0, 0.0)
( 26.0, 0.0)
( 64.0, 0.0)