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 = PTy2Comments
Informational error
Type Code
4 1 The input matrix is numerically singular.
Example
As an example, consider the 5 × 5 linear system:
Let
so that Ax1 = (23, 55, 107, 197, 278)T, and
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
x1
1
2 3
4 5
1.000 2.000
3.000 4.000 5.000
x2
1 2 3 4 5
5.000 4.000 3.000 2.000 1.000
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |