Determinant, Inertia, and Inverse of a UDU-Factored Symmetric Matrix in Packed Storage

Calling Sequence

CALL DSPDI	(DA, N, IPIVOT, DDET, INERT, DWORK, JOB)
CALL SSPDI	(SA, N, IPIVOT, SDET, INERT, SWORK, JOB)
CALL ZSPDI	(ZA, N, IPIVOT, ZDET, ZWORK, JOB)
CALL CSPDI	(CA, N, IPIVOT, CDET, CWORK, JOB)
void dspdi	(double *da, int n, int *ipivot, double *ddet, int *inert, int job)
void sspdi	(float *sa, int n, int *ipivot, float *sdet, int *inert, int job)
void zspdi	(doublecomplex *za, int n, int *ipivot, doublecomplex *zdet, int job)
void cspdi	(complex *ca, int n, int *ipivot, complex *cdet, int job)

Arguments

xA	On entry, the UDU factorization of the matrix A, as computed by xSPCO or SPFA. On exit, if the c digit of JOB 0, then A contains the upper triangle of the inverse of the original matrix A; otherwise unchanged.
N	Order of the original matrix A. N 0.
IPIVOT	Pivot vector as computed by xSPCO or xSPFA.
xDET	On exit, if the b digit of JOB 0, then DET contains the determinant of the matrix A. The determinant is stored as b × 10expon where b is stored in DET(1) and expon is stored in DET(2).
	1.0 |DET(1)| < 10.0 or DET(1) = 0.0. If the b digit of JOB = 0, DET is not referenced.
INERT	On exit, if the a digit of JOB 0 then INERT contains an integer triplet where:
	INERT(1) = number of positive eigenvalues
	INERT(2) = number of negative eigenvalues
	INERT(3) = number of zero eigenvalues
	If the a digit of JOB = 0 then INERT is not referenced.
xWORK	Scratch array with a dimension of N.
JOB	Integer in the form abc for real subroutines or bc for complex subroutines; determines operation the subroutine will perform:
	a 0	Compute the inertia.
	b 0	Compute the determinant.
	c 0	Compute the inverse.
	Note that the inverse should not be computed if xSPCO has set RCOND = 0 or if xSPFA has set INFO 0.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER IDODET, IDOINR, IDOINV, LENGTA, N

PARAMETER (IDODET = 10)

PARAMETER (IDOINR = 100)

PARAMETER (IDOINV = 1)

PARAMETER (N = 3)

PARAMETER (LENGTA = (N * N + N) / 2)

C

DOUBLE PRECISION A(LENGTA), DET(2), RCOND, WORK(N)

INTEGER INERT(3), IPIVOT(N), JOB

C

EXTERNAL DSPCO, DSPDI

C

C Initialize the array A to store in packed symmetric format

C the matrix A shown below.

C

C 1 0 4

C A = 0 2 0

C 4 0 3

C

DATA A / 3.0D0, 2.0D0, 2.0D0, 1.0D0, 1.0D0, 1.0D0 /

C

PRINT 1000

PRINT 1010, A(1), A(2), A(4)

PRINT 1010, A(2), A(3), A(5)

PRINT 1010, A(4), A(5), A(6)

CALL DSPCO (A, N, IPIVOT, RCOND, WORK)

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

PRINT 1020

END IF

JOB = IDODET + IDOINV + IDOINR

CALL DSPDI (A, N, IPIVOT, DET, INERT, WORK, JOB)

PRINT 1030, DET(1) * (10.0D0 ** DET(2))

PRINT 1040, INERT

PRINT 1050

PRINT 1010, A(1), A(2), A(4)

PRINT 1010, A(2), A(3), A(5)

PRINT 1010, A(4), A(5), A(6)

C

1000 FORMAT (1X, 'A:')

1010 FORMAT (3(3X, F5.2))

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

1030 FORMAT (/1X, 'Determinant of A: ', F7.3)

1040 FORMAT (1X, 'Inertia of A: <', I1, ',', I1, ',', I1, '>')

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

C

END

Sample Output

A:
3.00 2.00 1.00
2.00 2.00 1.00
1.00 1.00 1.00
Determinant of A: 1.000
Inertia of A: <3,0,0>
A**(-1):
1.00 -1.00 0.00
-1.00 2.00 -1.00
0.00 -1.00 2.00