LSLXD
Solves a sparse system of symmetric positive definite linear algebraic equations by Gaussian elimination.
Required Arguments
A — Vector of length NZ containing the nonzero coefficients in the lower triangle of the linear system. (Input)
The sparse matrix has nonzeroes only in entries (IROW (i), JCOL(i)) for i = 1 to NZ, and at this location the sparse matrix has value A(i).
IROW — Vector of length NZ containing the row numbers of the corresponding elements in the lower triangle of A. (Input)
Note IROW(i) ≥ JCOL(i), since we are only indexing the lower triangle.
JCOL — Vector of length NZ containing the column numbers of the corresponding elements in the lower triangle of A. (Input)
B — Vector of length N containing the right-hand side of the linear system. (Input)
X — Vector of length N containing the solution to the linear system. (Output)
Optional Arguments
N — Number of equations. (Input)
Default: N = size (B,1).
NZ — The number of nonzero coefficients in the lower triangle of the linear system. (Input)
Default: NZ = size (A,1).
ITWKSP — The total workspace needed. (Input)
If the default is desired, set ITWKSP to zero.
Default: ITWKSP = 0.
FORTRAN 90 Interface
Generic: CALL LSLXD (A, IROW, JCOL, B, X [, …])
Specific: The specific interface names are S_LSLXD and D_LSLXD.
FORTRAN 77 Interface
Single: CALL LSLXD (N, NZ, A, IROW, JCOL, B, ITWKSP, X)
Double: The double precision name is DLSLXD.
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
LSLXD solves a system of linear algebraic equations having a real, sparse and positive definite coefficient matrix. It first uses the routine
LSCXD to compute a symbolic factorization of a permutation of the coefficient matrix. It then calls
LNFXD to perform the numerical factorization. The solution of the linear system is then found using
LFSXD.
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
The routine LFSXD accepts b and the permutation vector which determines P. It then returns x.
Comments
1. Workspace may be explicitly provided, if desired, by use of L2LXD/DL2LXD. The reference is:
CALL L2LXD (N, NZ, A, IROW, JCOL, B, X, IPER, IPARAM, RPARAM, WK, LWK, IWK, LIWK)
The additional arguments are as follows:
IPER — Vector of length N containing the ordering.
IPARAM — Integer vector of length 4. See Comment 3.
RPARAM — Real vector of length 2. See Comment 3.
WK — Real work vector of length LWK.
LWK — The length of WK, LWK should be at least 2N + 6NZ.
IWK — Integer work vector of length LIWK.
LIWK — The length of IWK, LIWK should be at least 15N + 15NZ + 9.
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. |
3. If the default parameters are desired for L2LXD, then set IPARAM(1) to zero and call the routine L2LXD. Otherwise, if any nondefault parameters are desired for IPARAM or RPARAM, then the following steps should be taken before calling L2LXD.
CALL L4LXD (IPARAM, RPARAM)
Set nondefault values for desired IPARAM, RPARAM elements.
Note that the call to L4LXD will set IPARAM and RPARAM to their default values, so only nondefault values need to be set above. The arguments are as follows: |
IPARAM — Integer vector of length 4.
IPARAM(1) = Initialization flag.
IPARAM(2) = The numerical factorization method.
IPARAM(2) | Action |
0 | Multifrontal |
1 | Markowitz column search |
Default: 0.
IPARAM(3) = The ordering option.
IPARAM(3) | Action |
0 | Minimum degree ordering |
1 | User’s ordering specified in IPER |
Default: 0.
IPARAM(4) = The total number of nonzeros in the factorization matrix.
RPARAM — Real vector of length 2.
RPARAM(1) = The value of the largest diagonal element in the Cholesky factorization.
RPARAM(2) = The value of the smallest diagonal element in the Cholesky factorization.
If double precision is required, then DL4LXD is called and RPARAM is declared double precision.
Example
As an example consider the 5× 5 linear system:
Let xT = (1, 2, 3, 4, 5) so that Ax = (23, 55, 107, 197, 278)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 LSLXD_INT
USE WRRRN_INT
INTEGER N, NZ
PARAMETER (N=5, NZ=10)
!
INTEGER IROW(NZ), JCOL(NZ)
REAL A(NZ), B(N), X(N)
!
DATA A/10., 20., 1., 30., 4., 40., 2., 3., 5., 50./
DATA B/23., 55., 107., 197., 278./
DATA IROW/1, 2, 3, 3, 4, 4, 5, 5, 5, 5/
DATA JCOL/1, 2, 1, 3, 3, 4, 1, 2, 4, 5/
! Solve A * X = B
CALL LSLXD (A, IROW, JCOL, B, X)
! Print results
CALL WRRRN (’ x ’, X, 1, N, 1)
END
Output
x
1 2 3 4 5
1.000 2.000 3.000 4.000 5.000
