Solves a sparse system of symmetric positive definite linear algebraic equations by Gaussian elimination.
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)
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.
Generic: CALL LSLXD (A, IROW, JCOL, B, X [, ])
Specific: The specific interface names are S_LSLXD and D_LSLXD.
Single: CALL LSLXD (N, NZ, A, IROW, JCOL, B, ITWKSP, X)
Double: The double precision name is DLSLXD.
Consider the linear equation
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
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.
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
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 |
Users 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.
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
x
1
2 3
4 5
1.000 2.000
3.000 4.000 5.000
PHONE: 713.784.3131 FAX:713.781.9260 |