UDU Factorization of a Symmetric Matrix

Calling Sequence

CALL DSIFA	(DA, LDA, N, IPIVOT, INFO)
CALL SSIFA	(SA, LDA, N, IPIVOT, INFO)
CALL ZSIFA	(ZA, LDA, N, IPIVOT, INFO)
CALL CSIFA	(CA, LDA, N, IPIVOT, INFO)
void dsifa	(double *da, int lda, int n, int *ipivot, int *info)
void ssifa	(float *sa, int lda, int n, int *ipivot, int *info)
void zsifa	(doublecomplex *za, int lda, int n, int *ipivot, int *info)
void csifa	(complex *ca, int lda, 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. 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.
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO > 0	Returns a value k if the kth pivot block is singular to indicate that xSISL or xSIDI will divide by zero if called.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER LDA, N

PARAMETER (N = 4)

PARAMETER (LDA = N)

C

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

INTEGER ICOL, INFO, IPIVOT(N), IROW

C

EXTERNAL DSIFA, 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 6

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

C -.5 -1.5 -2.5 -2.5 6

C -.5 -1.5 -2.5 -3.5 12

C

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

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

$ -3.5D0 /

DATA B / 6.0D0, 1.2D1, 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 DSIFA (A, LDA, N, IPIVOT, INFO)

IF (INFO .EQ. 0) THEN

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

PRINT 1050

PRINT 1040, B

ELSE

PRINT 1060

END IF

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, 'A**(-1) * b:')

1060 FORMAT (/1X, 'A is too poorly conditioned.')

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:
6.0
12.0
6.0
12.0
A**(-1) * b:
-6.0
-12.0
12.0
-6.0