LFSXD

Solves a real sparse symmetric positive definite system of linear equations, given the Cholesky factorization of the coefficient matrix.

Required Arguments

N — Number of equations. (Input)

MAXSUB — Number of subscripts contained in array NZSUB as output from subroutine LSCXD/DLSCXD. (Input)

NZSUB — Vector of length MAXSUB containing the row subscripts for the off-diagonal nonzeros in the factor as output from subroutine LSCXD/DLSCXD. (Input)

INZSUB — Vector of length N + 1 containing pointers for NZSUB as output from subroutine LSCXD/DLSCXD. (Input)
The row subscripts of column J are stored from location INZSUB(J) to
INZSUB(J + 1) - 1.

MAXNZ — Total number of off-diagonal nonzeros in the Cholesky factor as output from subroutine LSCXD/DLSCXD. (Input)

RLNZ — Vector of length MAXNZ containing the off-diagonal nonzeros in the factor in column ordered format as output from subroutine LNFXD/DLNFXD. (Input)

ILNZ — Vector of length N + 1 containing pointers to RLNZ as output from subroutine LSCXD/DLSCXD. The nonzeros in column J of the factor are stored from location ILNZ(J) to ILNZ(J + 1) - 1. (Input) The values (RLNZ, ILNZ, NZSUB, INZSUB) give the off-diagonal nonzeros of the factor in a compressed subscript data format.

DIAGNL — Vector of length N containing the diagonals of the Cholesky factor as output from subroutine LNFXD/DLNFXD. (Input)

IPER — Vector of length N containing the ordering as output from subroutine LSCXD/DLSCXD. (Input)
IPER(I) = K indicates that the original row K is the new row I.

B — Vector of length N containing the right-hand side. (Input)

X — Vector of length N containing the solution. (Output)

FORTRAN 90 Interface

Generic: CALL LFSXD (N, MAXSUB, NZSUB, INZSUB, MAXNZ, RLNZ, ILNZ, DIAGNL, IPER, B, X)

Specific: The specific interface names are S_LFSXD and D_LFSXD.

FORTRAN 77 Interface

Single: CALL LFSXD (N, MAXSUB, NZSUB, INZSUB, MAXNZ, RLNZ, ILNZ, DIAGNL, IPER, B, X)

Double: The double precision name is DLFSXD.

Description

Consider the linear equation

Ax = b

where A is sparse, positive definite and symmetric. The sparse coordinate format for the matrix A requires one real and two integer vectors. The real array a contains all the nonzeros in the lower triangle of A including the diagonal. Let the number of nonzeros be nz. The two integer arrays irow and jcol, each of length nz, contain the row and column indices for these entries in A. That is

Airow(i), icol(i) = a(i), i = 1, …, nz

irow(i) ≥ jcol(i) i = 1, …, nz

with all other entries in the lower triangle of A zero.

The routine LFSXD computes the solution of the linear system given its Cholesky factorization. The factorization is performed by calling LSCXD followed by LNFXD. The routine LSCXD computes a minimum degree ordering or uses a user-supplied ordering to set up the sparse data structure for the Cholesky factor, L. Then the routine LNFXD produces the numerical entries in L so that we have

P APT= LLT

Here P is the permutation matrix determined by the ordering.

The numerical computations can be carried out in one of two ways. The first method performs the factorization using a multifrontal technique. This option requires more storage but in certain cases will be faster. The multifrontal method is based on the routines in Liu (1987). For detailed description of this method, see Liu (1990), also Duff and Reid (1983, 1984), Ashcraft (1987), Ashcraft et al. (1987), and Liu (1986, 1989). The second method is fully described in George and Liu (1981). This is just the standard factorization method based on the sparse compressed storage scheme.

Finally, the solution x is obtained by the following calculations:

1) Ly1 = Pb

2) LTy2 = y1

3) x = PTy2

Comments

Informational error

Type	Code	Description
4	1	The input matrix is numerically singular.

Example

As an example, consider the 5 × 5 linear system:

Let

x1T = (1,2,3,4,5)

so that Ax1 = (23, 55, 107, 197, 278)T, and

x2T = (5,4,3,2,1)

so that Ax2 = (55, 83, 103, 97, 82)T. The number of nonzeros in the lower triangle of A is nz = 10. The sparse coordinate form for the lower triangle of A is given by:

or equivalently by

USE LFSXD_INT

USE LNFXD_INT

USE LSCXD_INT

USE WRRRN_INT

INTEGER N, NZ, NRLNZ

PARAMETER (N=5, NZ=10, NRLNZ=10)

INTEGER IJOB, ILNZ(N+1), INVPER(N), INZSUB(N+1), IPER(N),&

IROW(NZ), ISPACE, ITWKSP, JCOL(NZ), MAXNZ, MAXSUB,&

NZSUB(3*NZ)

REAL A(NZ), B1(N), B2(N), DIAGNL(N), RLNZ(NRLNZ), RPARAM(2),&

X(N)

DATA A/10., 20., 1., 30., 4., 40., 2., 3., 5., 50./

DATA B1/23., 55., 107., 197., 278./

DATA B2/55., 83., 103., 97., 82./

DATA IROW/1, 2, 3, 3, 4, 4, 5, 5, 5, 5/

DATA JCOL/1, 2, 1, 3, 3, 4, 1, 2, 4, 5/

! Select minimum degree ordering

! for multifrontal method

IJOB = 3

! Use default workspace

ITWKSP = 0

MAXSUB = 3*NZ

CALL LSCXD (IROW, JCOL, NZSUB, INZSUB, MAXNZ, ILNZ, INVPER, &

MAXSUB=MAXSUB, IPER=IPER, ISPACE=ISPACE)

! Check if NRLNZ is large enough

IF (NRLNZ .GE. MAXNZ) THEN

! Choose multifrontal method

IJOB = 2

CALL LNFXD (A, IROW, JCOL, MAXSUB, NZSUB, INZSUB, MAXNZ, ILNZ,&

IPER, INVPER,ISPACE, DIAGNL, RLNZ, RPARAM, IJOB=IJOB)

! Solve A * X1 = B1

CALL LFSXD (N, MAXSUB, NZSUB, INZSUB, MAXNZ, RLNZ, ILNZ, DIAGNL,&

IPER, B1, X)

! Print X1

CALL WRRRN (’ x1 ’, X, 1, N, 1)

! Solve A * X2 = B2

CALL LFSXD (N, MAXSUB, NZSUB, INZSUB, MAXNZ, RLNZ, ILNZ, &

DIAGNL, IPER, B2, X)

! Print X2

CALL WRRRN (’ x2 ’ X, 1, N, 1)

END IF

END

Output

1 2 3 4 5

1.000 2.000 3.000 4.000 5.000

1 2 3 4 5

5.000 4.000 3.000 2.000 1.000