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
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
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |