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.
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.
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.
Optional Arguments
N — Number of equations. (Input)
Default: N = size (IPER,1).
NZ — The number of nonzero coefficients in the linear system. (Input)
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.
ITWKSP — The total workspace needed. (Input)
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
