Solution to a Linear System in a UDU-Factored Symmetric Matrix

Calling Sequence

CALL DSISL	(DA, LDA, N, IPIVOT, DB)
CALL SSISL	(SA, LDA, N, IPIVOT, SB)
CALL ZSISL	(ZA, LDA, N, IPIVOT, ZB)
CALL CSISL	(CA, LDA, N, IPIVOT, CB)
void dsisl	(double *da, int lda, int n, int *ipivot, double *db)
void ssisl	(float *sa, int lda, int n, int *ipivot, float *sb)
void zsisl	(doublecomplex *za, int lda, int n, int *ipivot, doublecomplex *zb)
void csisl	(complex *ca, int lda, int n, int *ipivot, complex *cb)

Arguments

xA	UDU factorization of the matrix A, as computed by xSICO or xSIFA.
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	Pivot vector as computed by xSICO or xSIFA.
xB	On entry, the right-hand side vector b. On exit, the solution vector x.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER LDA, N

PARAMETER (N = 4)

PARAMETER (LDA = N)

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 / 4*-5.0D-1, -5.0D-1, 3*-1.5D0, -5.0D-1,

$ -1.5D0, 2*-2.5D0, -5.0D-1, -1.5D0, -2.5D0, -3.5D0 /

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)

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

PRINT 1050

END IF

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

PRINT 1060

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, 'A may be singular to working precision.')

1060 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
-0.5 -1.5 -1.5 -1.5
-0.5 -1.5 -2.5 -2.5
-0.5 -1.5 -2.5 -3.5
b:
12.0
6.0
6.0
12.0
A**(-1) * b:
-30.0
6.0
6.0
-6.0