LNFXD

Computes the numerical Cholesky factorization of a sparse symmetrical matrix A.

Required Arguments

A — Vector of length NZ containing the nonzero coefficients of the lower triangle of the linear system. (Input)

IROW — Vector of length NZ containing the row numbers of the corresponding elements in the lower triangle of A. (Input)

JCOL — Vector of length NZ containing the column numbers of the corresponding elements in the lower triangle of A. (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 nonzeros in the Cholesky factor in compressed format 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 for the nonzeros in column J are stored from location INZSUB (J) to INZSUB(J + 1) - 1.

The row subscripts for the nonzeros in column J are stored from location INZSUB (J) to INZSUB(J + 1) - 1.

MAXNZ — Length of RLNZ as output from subroutine LSCXD/DLSCXD. (Input)

ILNZ — Vector of length N + 1 containing pointers to the Cholesky factor as output from subroutine LSCXD/DLSCXD. (Input)

The row subscripts for the nonzeros in column J of the factor are stored from location ILNZ(J) to ILNZ(J + 1) - 1. (ILNZ, NZSUB, INZSUB) sets up the compressed data structure in column ordered form for the Cholesky factor.

The row subscripts for the nonzeros in column J of the factor are stored from location ILNZ(J) to ILNZ(J + 1) - 1. (ILNZ, NZSUB, INZSUB) sets up the compressed data structure in column ordered form for the Cholesky factor.

IPER — Vector of length N containing the permutation as output from subroutine LSCXD/DLSCXD. (Input)

INVPER — Vector of length N containing the inverse permutation as output from subroutine LSCXD/DLSCXD. (Input)

ISPACE — The storage space needed for the stack of frontal matrices as output from subroutine LSCXD/DLSCXD. (Input)

DIAGNL — Vector of length N containing the diagonal of the factor. (Output)

RLNZ — Vector of length MAXNZ containing the strictly lower triangle nonzeros of the Cholesky factor. (Output)

RPARAM — Parameter vector containing factorization information. (Output)

RPARAM(1) = smallest diagonal element.

RPARAM(2) = largest diagonal element.

RPARAM(1) = smallest diagonal element.

RPARAM(2) = largest diagonal element.

Optional Arguments

N — Number of equations. (Input)

Default: N = size (IPER,1).

Default: N = size (IPER,1).

NZ — The number of nonzero coefficients in the linear system. (Input)

Default: NZ = size (A,1).

Default: NZ = size (A,1).

IJOB — Integer parameter selecting factorization method. (Input)

IJOB = 1 yields factorization in sparse column format.

IJOB = 2 yields factorization using multifrontal method.

Default: IJOB = 1.

IJOB = 1 yields factorization in sparse column format.

IJOB = 2 yields factorization using multifrontal method.

Default: IJOB = 1.

ITWKSP — The total workspace needed. (Input)

If the default is desired, set ITWKSP to zero.

Default: ITWKSP = 0.

If the default is desired, set ITWKSP to zero.

Default: ITWKSP = 0.

FORTRAN 90 Interface

Generic: CALL LNFXD (A, IROW, JCOL, MAXSUB, NZSUB, INZSUB, MAXNZ, ILNZ, IPER, INVPER, ISPACE, DIAGNL, RLNZ, RPARAM [, …])

Specific: The specific interface names are S_LNFXD and D_LNFXD.

FORTRAN 77 Interface

Single: CALL LNFXD (N, NZ, A, IROW, JCOL, IJOB, ITWKSP, MAXSUB, NZSUB, INZSUB, MAXNZ, ILNZ, IPER, INVPER, ISPACE, ITWKSP, DIAGNL, RLNZ, RPARAM)

Double: The double precision name is DLNFXD.

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 LNFXD produces the Cholesky factorization of P APTgiven the symbolic factorization of A which is computed by LSCXD. That is, this routine computes L which satisfies

P APT= LLT

The diagonal of L is stored in DIAGNL and the strictly lower triangular part of L is stored in compressed subscript form in R = RLNZ as follows. The nonzeros in the j-th column of L are stored in locations R(i), …, R(i + k) where i = ILNZ(j) and k = ILNZ(j + 1) − ILNZ(j) − 1. The row subscripts are stored in the vector NZSUB from locations INZSUB(j) to INZSUB(j) + k.

The numerical computations can be carried out in one of two ways. The first method (when IJOB = 2) 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 (when IJOB = 1) is fully described in George and Liu (1981). This is just the standard factorization method based on the sparse compressed storage scheme.

Comments

1. Workspace may be explicitly provided by use of L2FXD/DL2FXD . The reference is:

CALL L2FXD (N, NZ, A, IROW, JCOL, IJOB, MAXSUB, NZSUB, INZSUB, MAXNZ, ILNZ, IPER, INVPER, ISPACE, DIAGNL, RLNZ, RPARAM, WK, LWK, IWK, LIWK)

The additional arguments are as follows:

WK — Real work vector of length LWK.

LWK — The length of WK, LWK should be at least N + 3NZ.

IWK — Integer work vector of length LIWK.

LIWK — The length of IWK, LIWK should be at least 2N.

Note that the parameter ITWKSP is not an argument to this routine.

2. Informational errors

Type | Code | Description |
---|---|---|

4 | 1 | The coefficient matrix is not positive definite. |

4 | 2 | A column without nonzero elements has been found in the coefficient matrix. |

Example

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

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

We first call LSCXD to produce the symbolic information needed to pass on to LNFXD. Then call LNFXD to factor this matrix. The results are displayed below.

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, JCOL(NZ), MAXNZ, MAXSUB,&

NZSUB(3*NZ)

REAL A(NZ), DIAGNL(N), RLNZ(NRLNZ), RPARAM(2) , R(N,N)

!

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

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

MAXSUB = 3*NZ

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

MAXSUB=MAXSUB)

! 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)

! Print results

CALL WRRRN (' diagnl ', DIAGNL, NRA=1, NCA=N, LDA=1)

CALL WRRRN (' rlnz ', RLNZ, NRA= 1, NCA= MAXNZ, LDA= 1)

END IF

!

! Construct L matrix

DO I=1,N

! Diagonal

R(I,I) = DIAG(I)

IF (ILNZ(I) .GT. MAXNZ) GO TO 50

! Find elements of RLNZ for this column

ISTRT = ILNZ(I)

ISTOP = ILNZ(I+1) - 1

! Get starting index for NZSUB

K = INZSUB(I)

DO J=ISTRT, ISTOP

! NZSUB(K) is the row for this element of

RLNZ

R((NZSUB(K)),I) = RLNZ(J)

K = K + 1

END DO

END DO

50 CONTINUE

CALL WRRRN ('L', R, NRA=N, NCA=N)

END

Output

diagnl

1 2 3 4 5

4.472 3.162 7.011 6.284 5.430

rlnz

1 2 3 4 5 6

0.6708 0.6325 0.3162 0.7132 -0.0285 0.6398

L

1 2 3 4 5

1 4.472 0.000 0.000 0.000 0.000

2 0.000 3.162 0.000 0.000 0.000

3 0.671 0.632 7.011 0.000 0.000

4 0.000 0.000 0.713 6.284 0.000

5 0.000 0.316 -0.029 0.640 5.430