LSLZD
Solves a complex sparse Hermitian positive definite system of linear equations by Gaussian elimination.
Required Arguments
A — Complex 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 — Complex vector of length N containing the right-hand side of the linear system. (Input)
X — Complex 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 LSLZD (A, IROW, JCOL, B, X [, …])
Specific: The specific interface names are S_LSLZD and D_LSLZD.
FORTRAN 77 Interface
Single: CALL LSLZD (N, NZ, A, IROW, JCOL, B, ITWKSP, X)
Double: The double precision name is DLSLZD.
Description
Consider the linear equation
Ax = b
where A is sparse, positive definite and Hermitian. The sparse coordinate format for the matrix A requires one complex and two integer vectors. The complex 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
LSLZD solves a system of linear algebraic equations having a complex, sparse, Hermitian 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
LNFZD to perform the numerical factorization. The solution of the linear system is then found using
LFSZD.
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 LNFZD produces the numerical entries in L so that we have
P APT = LLH
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) LH y2 = y1
3) x = PT y2
The routine LFSZD 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 L2LZD/DL2LZD. The reference is:
CALL L2LZD (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 — Complex 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 for 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 L2LZD, then set IPARAM(1) to zero and call the routine L2LZD. Otherwise, if any nondefault parameters are desired for IPARAM or RPARAM, then the following steps should be taken before calling L2LZD.
CALL L4LZD (IPARAM, RPARAM)
Set nondefault values for desired IPARAM, RPARAM elements.
Note that the call to L4LZD 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 | Sparse column |
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 absolute value of the largest diagonal element in the Cholesky factorization.
RPARAM(2) = The absolute value of the smallest diagonal element in the Cholesky factorization.
If double precision is required, then DL4LZD is called and RPARAM is declared double precision.
Example
As an example, consider the 3 × 3 linear system:
Let xT = (1 + i, 2 + 2i, 3 + 3i) so that Ax = (−2 + 2i, 5 + 15i, 36 + 28i)T. The number of nonzeros in the lower triangle of A is nz = 5. The sparse coordinate form for the lower triangle of A is given by:
or equivalently by
USE LSLZD_INT
USE WRCRN_INT
INTEGER N, NZ
PARAMETER (N=3, NZ=5)
!
INTEGER IROW(NZ), JCOL(NZ)
COMPLEX A(NZ), B(N), X(N)
!
DATA A/(2.0,0.0), (4.0,0.0), (10.0,0.0), (-1.0,-1.0), (1.0,-2.0)/
DATA B/(-2.0,2.0), (5.0,15.0), (36.0,28.0)/
DATA IROW/1, 2, 3, 2, 3/
DATA JCOL/1, 2, 3, 1, 2/
! Solve A * X = B
CALL LSLZD (A, IROW, JCOL, B, X)
! Print results
CALL WRCRN (’ x ’, X, 1, N, 1)
END
Output
x
1 2 3
( 1.000, 1.000) ( 2.000, 2.000) ( 3.000, 3.000)
