UDU Factorization of a Symmetric Matrix in Packed Storage

Calling Sequence

CALL DSPFA	(DA, N, IPIVOT, INFO)
CALL SSPFA	(SA, N, IPIVOT, INFO)
CALL ZSPFA	(ZA, N, IPIVOT, INFO)
CALL CSPFA	(CA, N, IPIVOT, INFO)
void dspfa	(double *da, int n, int *ipivot, int *info)
void sspfa	(float *sa, int n, int *ipivot, int *info)
void zspfa	(doublecomplex *za, int n, int *ipivot, int *info)
void cspfa	(complex *ca, int n, int *ipivot, int *info)

Arguments

xA	On entry, the upper triangle of the matrix A. On exit, a UDU factorization of the matrix A.
N	Order of the matrix A. N 0.
IPIVOT	On exit, a vector of pivot indices.
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO > 0	Returns a value k if the kth pivot block is singular to indicate that xSPSL or xSPDI will divide by zero if called.

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), WORK(N)

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

C

EXTERNAL DSPDI, DSPFA

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 DSPFA (A, N, IPIVOT, INFO)

IF (INFO .EQ. 0) THEN

JOB = 111

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

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

PRINT 1030, INERT

PRINT 1040

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

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

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

ELSE

PRINT 1050

END IF

C

1000 FORMAT (1X, 'A:')

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

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

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

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

1050 FORMAT (1X, 'A is too ill-conditioned.')

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