SOLVER(3) MathKeisan SOLVER SOLVER(3)
NAME
Direct Sparse Solver for Symmetric Positive Definite Linear Systems -
ldlt_preprocess, ldlt_fact, ldlt_solve, ldlt_getdiag, ldlt_refine,
ldlt_free
SYNOPSIS
INTEGER token, n, colptr(n+1), rowind(nnz), Lnnz, order
CALL LDLT_PREPROCESS ( token, n, colptr, rowind, Lnnz, order )
INTEGER token, colptr(n+1), rowind(nnz)
REAL*8 nzvals(nnz), flops
CALL LDLT_FACT ( token, colptr, rowind, nzvals, flops )
INTEGER token
REAL*8 x(n), rhs(n), flops
CALL LDLT_SOLVE ( token, x, rhs, flops )
INTEGER token
REAL*8 diag(n)
CALL LDLT_GETDIAG ( token, diag )
INTEGER token, colptr(n+1), rowind(nnz)
REAL*8 nzvals(nnz), rhs(n), x(n), newx(n), flops
CALL LDLT_REFINE ( token, colptr, rowind, nzvals, rhs, x, newx, flops )
INTEGER token
CALL LDLT_FREE ( token )
DESCRIPTION
ldlt is a sparse direct solver for symmetric positive definite linear
systems. There are subroutines for setup and memory allocation, fac-
torization, and forward and back substitution. Additional routines
allow for the freeing of all allocated memory, for access to the diago-
nal of the factorization, and for performing one step of iterative
refinement.
The interface allows for using one factorization for multiple right-
hand side vectors, or one reordering for multiple factorizations. In
the first case, factor once and solve multiple times. In the second
case, reorder once and factorize multiple times.
Under certain circumstances the error in the calculated solution may be
too large, such as when the matrix has a large condition number. In
these cases, the accuracy of the initially computed solution may be
improved by one step of iterative refinement, based on the original
matrix and the factorization.
PARALLEL FACTORIZATION
ldlt allows the numerical factorization to be done in parallel. The
algorithm used for parallel numerical factorization is slightly differ-
ent from that for the serial factorization. The user does not need to
change the code to select the parallel factorization routine. All that
needs to be done is to set the environment variable OPM_NUM_THREADS to
the number of parallel threads desired. If OPM_NUM_THREADS is not set,
the serial numerical factorization will be invoked.
SPARSE MATRIX FORMAT
In order to use ldlt, one must store the sparse matrix in compressed
column format (This is the format used in the Harwell-Boeing matrix
collection). The matrix must be symmetric, and only the lower triangu-
lar part of the matrix is required.
The nonzero entries are stored column by column contiguously in a dou-
ble precision array nzvals. Row indicies for entries in nzvals are
stored in an integer array rowind. Both nzvals and rowind contain nnz
entries, where nnz is the number of nonzeros in the lower triangular
part of the sparse matrix, including the diagonal. An integer array
colptr of length n+1 is used to keep the starting point of each column.
More precisely, colptr(j) gives the starting location of the j-th col-
umn in nzvals and rowind.
Below is a simple example to illustrate the compressed column storage
scheme.
| 10.0 |
| 40.0 |
A = | 2.0 60.0 |
| 3.0 5.0 80.0 |
| 7.0 90.0 |
This matrix is stored as follows:
+-----+-----+-----+-----+-----+-----+-----+-----+-----+
nzvals = | 10.0| 2.0 | 3.0 | 40.0| 5.0 | 60.0| 7.0 | 80.0| 90.0|
+-----+-----+-----+-----+-----+-----+-----+-----+-----+
+---+---+---+---+---+---+---+---+---+
rowind = | 1 | 3 | 4 | 2 | 4 | 3 | 5 | 4 | 5 |
+---+---+---+---+---+---+---+---+---+
+----+----+----+----+----+----+
colptr = | 1 | 4 | 6 | 8 | 9 | 10 |
+----+----+----+----+----+----+
ARGUMENTS
Input
INTEGER token
Matrix identifier for the linear system to be solved.
INTEGER n
Dimension of the linear system.
INTEGER order
Method used to order the original sparse matrix.
order = 0: C-version MMD (provides supernode info)
order = 1: FORTRAN-version MMD
order = 2: Node-based Nested Dissection
order = 3: Edge-based Nested Dissection
Node-based Nested Dissection frequently gives the best ordering,
so it is suggested that you use order=2.
INTEGER colptr(n+1)
Array containing the index of the starting point of all columns
in nzvals and rowind. For example, colptr(j) indicates the
index in nzvals and in rowind corresponding to the first entry
in column j.
INTEGER rowind(nnz)
Array containing the row indices of the nonzero entries in the
sparse matrix. nnz is the number of non-zero entries.
REAL*8 nzvals(nnz)
Array containing the nonzero entries of the sparse matrix.
REAL*8 rhs(n)
Array containing the right-hand side vector.
REAL*8 x(n)
Array containing the previously computed solution to the linear
system. For ldlt_refine only.
Output
INTEGER Lnnz
Number of nonzero entries in the LDLT factorization of the given
sparse matrix.
REAL*8 x(n)
Array containing the solution to the linear system.
REAL*8 flops
Number of floating point operations performed.
REAL*8 diag(n)
Array containing the diagonal entries of the Cholesky factoriza-
tion of the given sparse matrix.
REAL*8 newx(n)
Array containing the new solution after one step of iterative
refinement.
EXAMPLE CODE
The example below shows how to use the ldlt routines to solve a single
sparse linear system.
program example
c
integer maxn, maxnnz
parameter (maxn = 10, maxnnz = 100)
c
integer n, Lnnz, token, i
double precision fflops, sflops, rflops
integer colptr(maxn+1), rowind(maxnnz)
double precision nzvals(maxnnz),
& rhs(maxn), x(maxn), newx(maxn)
c
token = 0
n = 5
order = 2
c
c === original sparse matrix ===
c
colptr(1) = 1
colptr(2) = 4
colptr(3) = 6
colptr(4) = 8
colptr(5) = 9
colptr(6) = 10
c
rowind(1) = 1
rowind(2) = 3
rowind(3) = 4
rowind(4) = 2
rowind(5) = 4
rowind(6) = 3
rowind(7) = 5
rowind(8) = 4
rowind(9) = 5
c
nzvals(1) = 10.0
nzvals(2) = 1.0
nzvals(3) = 1.0
nzvals(4) = 10.0
nzvals(5) = 1.0
nzvals(6) = 10.0
nzvals(7) = 1.0
nzvals(8) = 10.0
nzvals(9) = 10.0
c
c === preprocess ===
c
call ldlt_preprocess(token, n, colptr, rowind, Lnnz, order)
c
c === numerical factorization ===
c
call ldlt_fact(token, colptr, rowind, nzvals, fflops)
c
c === construct the right-hand side ===
c
do i = 1, n
rhs(i) = 1.0
end do
call ldlt_solve(token, x, rhs, sflops)
c
c === refine the solution ===
c
call ldlt_refine (token, colptr, rowind, nzvals, rhs,
& x, newx, rflops)
c
c === free working storage ===
c
call ldlt_free(token)
end
If you have any problems with this library, please contact us at:
http://www.mathkeisan.com/
SEE ALSO
mathkeisan(3), MathKeisan User's Guide at:
http://www.mathkeisan.com/
The command:
man mathkeisan
gives more information on linking.
MathKeisan SOLVER(3)