Direct Solvers¶
The Low-Level C Modules¶
The superlu
Module¶
The superlu
module interfaces the SuperLU library to make it usable by
Python code. SuperLU is a software package written in C, that is able to compute
an LU-factorisation of a general non-symmetric sparse matrix with
partial pivoting.
The superlu
module exports a single function, called factorize
.
-
superlu.
factorize
(A, **kwargs)¶ The factorize function computes an LU-factorisation of the matrix
A
.Parameters: A: A csr_mat
object that represents the matrix to be factorized.Keywords: diag_pivot_thresh: the partial pivoting threshold, in the interval [0,1].
diag_pivot_thresh=0
corresponds to no pivoting.diag_pivot_thresh=1
corresponds to partial pivoting (default: 1.0).drop_tol: the drop tolerance, in the interval [0,1].
drop_tol=0
corresponds to the exact factorization (default: 0.0).relax: the degree of relaxing supernodes (default: 1).
panel_size: the maximum number of columns that form a panel (default: 10).
permc_spec: the matrix ordering used to control sparsity of the factors:
- natural ordering
- MMD applied to the structure of \mathbf{A}^T\mathbf{A}
- MMD applied to the structure of \mathbf{A}^T + \mathbf{A}
- COLAMD, approximate minimum degree column ordering
(default: 2).
Return type: an object of type
superlu_context
. This object encapsulates the L and U factors ofA
(see below).
Note
The drop_tol
has no effect in SuperLU version 2.0 and below. In SuperLU
version 3.0 and above, the default value of permc_spec
is 3.
superlu_context
Object Attributes and Methods¶
-
class
superlu.
superlu_context
¶ An abstract encapsulation of the LU factorization of a matrix by SuperLU.
-
shape
¶ A 2-tuple describing the dimension of the matrix factorized. It is equal to
A.shape
.
-
nnz
¶ The
nnz
attribute holds the total number of nonzero entries stored in both the L and U factors.
-
solve
(b, x, trans)¶ The
solve
method accepts two rank-1 NumPy arraysb
andx
of appropriate size and assigns the solution of the linear system \mathbf{A}\mathbf{x} = \mathbf{b} tox
. If the optional parametertrans
is set to the string'T'
, the transposed system \mathbf{A}^T\mathbf{x} = \mathbf{b} is solved instead.
-
Example: Solving a 2D Poisson System with SuperLU¶
Let’s now solve the 2D Poisson system \mathbf{A} \mathbf{x} = \mathbf{1} using an LU factorization. Here, \mathbf{A} is the 2D Poisson matrix, introduced in Low-Level Sparse Matrix Types and \mathbf{1} is a vector with all entries equal to one.
The Python solution for this task looks as follows:
from pysparse.sparse import spmatrix
from pysparse.direct import superlu
import numpy
n = 100
A = poisson2d_sym_blk(n)
b = numpy.ones(n*n)
x = numpy.empty(n*n)
LU = superlu.factorize(A.to_csr(), diag_pivot_thresh=0.0)
LU.solve(b, x)
The code makes use of the Python function poisson2d_sym_blk()
, which
was defined in Low-Level Sparse Matrix Types.
Example: An Incomplete LU Factorization Preconditioner¶
Warning
SuperLU 3.0 and above accept a drop_tol
argument although the source
files mention that incomplete factorization is not implemented. Therefore,
changing drop_tol
has no effect on the factorization at the moment and we
must wait for it to be implemented. In the meantime, we can still demonstrate
in this section how to implement an incomplete factorization preconditioner
in Pysparse, even though in the present situation, it will be a complete
factorization preconditioner!
Versions of SuperLU above 3.0 accept the drop_tol
argument that allows the
computation of incomplete factors, realizing a tradeoff between computational
cost and factor density. The following example show how to use an incomplete LU
factorization as a preconditioner in any of the iterative methods of the
itsolvers
module:
from pysparse.tools import poisson
from pysparse.direct import superlu
from pysparse.itsolvers import krylov
import numpy
class ILU_Precon:
"""
A preconditioner based on an
incomplete LU factorization.
Input: A matrix in CSR format.
Keyword argument: Drop tolerance.
"""
def __init__(self, A, drop=1.0e-3):
self.LU = superlu.factorize(A, drop_tol=drop)
self.shape = self.LU.shape
def precon(self, x, y):
self.LU.solve(x,y)
n = 300
A = poisson.poisson2d_sym_blk(n).to_csr() # Convert right away
b = numpy.ones(n*n)
x = numpy.empty(n*n)
K = ILU_Precon(A)
info, niter, relres = krylov.pcg(A, b, x, 1e-12, 2000, K)
Note
Note that the 2D Poisson matrix is symmetric and positive definite, although barely. Indeed its smallest eigenvalue is 2 (1 - \cos(\pi/(n+1))) \approx (\pi/(n+1))^2. Therefore, a Cholesky factorization would be more appropriate. In the future, we intend to interface the Cholmod library.
The umfpack
Module¶
The umfpack
module interfaces the UMFPACK library to make it usable by
Python code. UMFPACK is a software package written in C, that is able to compute
an LU factorization of a general non-symmetric sparse matrix with
partial pivoting.
Note
The major difference with the superlu
modules is
that umfpack
receives a matrix in ll_mat
format instead of
csr_mat format.
The umfpack
module exports a single function, called factorize
.
-
superlu.
factorize
(A, **kwargs) The factorize function computes an LU-factorisation of the matrix
A
.Parameters: A: A ll_mat
object that represents the matrix to be factorized.Keywords: strategy: Pivoting strategy. Possible values are:
- “UMFPACK_STRATEGY_AUTO”
- “UMFPACK_STRATEGY_UNSYMMETRIC”
- “UMFPACK_STRATEGY_SYMMETRIC”
- “UMFPACK_STRATEGY_2BY2”
tol2by2: Tolerance for the 2-by-2 strategy.
scale: Scaling used during factorization. Possible values are:
- “UMFPACK_SCALE_NONE”
- “UMFPACK_SCALE_SUM”
- “UMFPACK_SCALE_MAX”
tolpivot: Relative pivot tolerance for threshold partial pivoting with row interchanges.
tolsympivot: If diagonal pivoting is attempted, this parameter controls when the diagonal is selected in a given pivot column.
Return type: an object of type
umfpack_context
. This object encapsulates the L and U factors ofA
(see below).
umfpack_context
Object Attributes and Methods¶
-
class
superlu.
umfpack_context
¶ An abstract encapsulation of the LU factorization of a matrix by UMFPACK.
-
shape
¶ A 2-tuple describing the dimension of the matrix factorized. It is equal to
A.shape
.
-
nnz
¶ The
nnz
attribute holds the total number of nonzero entries stored in the input matrix. It is equal toA.nnz
. To obtain the number of nonzero element in the factors, seelunz()
.
-
solve
(b, x, method, irsteps)¶ Parameters: b: The right-hand side of the system \mathbf{A} x = b as a Numpy array. x: A Numpy array to hold the solution of \mathbf{A} x = b. method: (optional) Different systems may be solved by setting the method
argument appropriately. See the documentation of thepysparseUmfpackSolver
below for more details.irsteps: (optional) The number of iterative refinement steps to perform.
-
lu
()¶ Return the factors, permutation and scaling information. See the documentation of the
pysparseUmfpackSolver
below for more details.
-
lunz
()¶ Return the number of nonzeros in factors, i.e., in \mathbf{L} + \mathbf{U}.
-
Example: Solving a 2D Poisson System with UMFPACK¶
We now solve again the 2D Poisson system \mathbf{A} \mathbf{x} = \mathbf{1} using an LU factorization. Here, \mathbf{A} is the 2D Poisson matrix, introduced in Low-Level Sparse Matrix Types and \mathbf{1} is a vector with all entries equal to one.
The Python solution using UMFPACK looks as follows:
from pysparse.sparse import spmatrix
from pysparse.direct import umfpack
import numpy
n = 100
A = poisson2d_sym_blk(n)
b = numpy.ones(n*n)
x = numpy.empty(n*n)
LU = umfpack.factorize(A, strategy="UMFPACK_STRATEGY_SYMMETRIC")
LU.solve(b, x)
The code makes use of the Python function poisson2d_sym_blk()
, which
was defined in Low-Level Sparse Matrix Types.
Higher-Level Python Interfaces¶
This section anticipates on Higher-Level Sparse Matrix Classes and shows usage of higher-level interfaces to the LU factorization packages.
The Abstract directSolver
Module¶
A framework for solving sparse linear systems of equations using a direct factorization.
The pysparseSuperLU
Module: A Higher-Level SuperLU Interface¶
A framework for solving sparse linear systems of equations using an LU factorization, by means of the supernodal sparse LU factorization package SuperLU ([DEGLL99], [DGL99], [LD03]).
This package is appropriate for factorizing sparse square unsymmetric or rectangular matrices.
See [SLU] for more information.
References:
[DEGLL99] | J. W. Demmel, S. C. Eisenstat, J. R. Gilbert, X. S. Li and J. W. H. Liu, A supernodal approach to sparse partial pivoting, SIAM Journal on Matrix Analysis and Applications 20(3), pp. 720-755, 1999. |
[DGL99] | J. W. Demmel, J. R. Gilbert and X. S. Li, An Asynchronous Parallel Supernodal Algorithm for Sparse Gaussian Elimination, SIAM Journal on Matrix Analysis and Applications 20(4), pp. 915-952, 1999. |
[LD03] | X. S. Li and J. W. Demmel, SuperLU_DIST: A Scalable Distributed-Memory Sparse Direct Solver for Unsymmetric Linear Systems, ACM Transactions on Mathematical Software 29(2), pp. 110-140, 2003. |
[SLU] | http://crd.lbl.gov/~xiaoye/SuperLU |
-
class
pysparse.direct.pysparseSuperLU.
PysparseSuperLUSolver
(A, **kwargs)¶ Bases:
pysparse.direct.directSolver.PysparseDirectSolver
PysparseSuperLUSolver is a wrapper class around the SuperLu library for the factorization of full-rank n-by-m matrices. Only matrices with real coefficients are currently supported.
Parameters: A: The matrix to be factorized, supplied as a PysparseMatrix instance. Keywords: symmetric: a boolean indicating that the user wishes to use symmetric mode. In symmetric mode,
permc_spec=2
must be chosen anddiag_pivot_thresh
must be small, e.g., 0.0 or 0.1. Since the value ofdiag_pivot_thresh
is up to the user, settingsymmetric
toTrue
does not automatically setpermc_spec
anddiag_pivot_thresh
to appropriate values.diag_pivot_thresh: a float value between 0 and 1 representing the threshold for partial pivoting (0 = no pivoting, 1 = always perform partial pivoting). Default: 1.0.
drop_tol: the value of a drop tolerance, between 0 and 1, if an incomplete factorization is desired (0 = exact factorization). This keyword does not exist if using SuperLU version 2.0 and below. In more recent version of SuperLU, the keyword is accepted but has no effect. Default: 0.0
relax: an integer controling the degree of relaxing supernodes. Default: 1.
panel_size: an integer specifying the maximum number of columns to form a panel. Default: 10.
permc_spec: an integer specifying the ordering strategy used during the factorization.
- natural ordering,
- MMD applied to the structure of \mathbf{A}^T \mathbf{A}
- MMD applied to the structure of \mathbf{A}^T + \mathbf{A}
- COLAMD.
Default: 2.
-
LU
¶ A
superlu_context
object encapsulating the factorization.
-
factorizationTime
¶ The CPU time to perform the factorization.
-
solutionTime
¶ The CPU time to perform the forward and backward sweeps.
-
lunz
¶ The number of nonzero elements in the factors L and U together after a call to
fetch_lunz()
.
-
fetch_factors
()¶ Not yet available.
Example: The 2D Poisson System with SuperLU¶
The solution of a 2D Poisson system with PysparseSuperLUSolver
may look
like this:
from pysparse.sparse.pysparseMatrix import PysparseMatrix
from pysparse.direct.pysparseSuperLU import PysparseSuperLUSolver
from pysparse.tools.poisson_vec import poisson2d_sym_blk_vec
from numpy import ones
from numpy.linalg import norm
n = 200
A = PysparseMatrix( matrix=poisson2d_sym_blk_vec(n) )
x_exact = ones(n*n)/n
b = A * x_exact
LU = PysparseSuperLUSolver(A)
LU.solve(b)
print 'Factorization time: ', LU.factorizationTime
print 'Solution time: ', LU.solutionTime
print 'Error: ', norm(LU.sol - x_exact)/norm(x_exact)
The above script produces the output:
Factorization time: 0.494116
Solution time: 0.017096
Error: 2.099685128150953e-14
Note that this example uses the vectorized Poisson constructors of Low-Level Sparse Matrix Types.
The pysparseUmfpack
Module: A Higher-Level UMFPACK Interface¶
A framework for solving sparse linear systems of equations using an LU factorization, by means of the unsymmetric multifrontal sparse LU factorization package UMFPACK ([D04a], [D04b], [DD99], [DD97]).
This package is appropriate for factorizing sparse square unsymmetric or rectangular matrices.
See [UMF] for more information.
References:
[D04a] | T. A. Davis, A column pre-ordering strategy for the unsymmetric-pattern multifrontal method, ACM Transactions on Mathematical Software, 30(2), pp. 165-195, 2004. |
[D04b] | T. A. Davis, Algorithm 832: UMFPACK, an unsymmetric-pattern multifrontal method, ACM Transactions on Mathematical Software, 30(2), pp. 196-199, 2004. |
[DD99] | T. A. Davis and I. S. Duff, A combined unifrontal/multifrontal method for unsymmetric sparse matrices, ACM Transactions on Mathematical Software, 25(1), pp. 1-19, 1999. |
[DD97] | T. A. Davis and I. S. Duff, An unsymmetric-pattern multifrontal method for sparse LU factorization, SIAM Journal on Matrix Analysis and Applications, 18(1), pp. 140-158, 1997. |
[UMF] | http://www.cise.ufl.edu/research/sparse/umfpack |
-
class
pysparse.direct.pysparseUmfpack.
PysparseUmfpackSolver
(A, **kwargs)¶ Bases:
pysparse.direct.directSolver.PysparseDirectSolver
PysparseUmfpackSolver is a wrapper class around the UMFPACK library for the factorization of full-rank n-by-m matrices. Only matrices with real coefficients are currently supported.
Parameters: A: A PysparseMatrix instance representing the matrix to be factorized. Keywords: strategy: string that specifies what kind of ordering and pivoting strategy UMFPACK should use. Valid values are ‘auto’, ‘unsymmetric’, ‘symmetric’ and ‘2by2’. Default: ‘auto’ tol2by2: tolerance for the 2 by 2 strategy. Default: 0.1 scale: string that specifies the scaling UMFPACK should use. Valid values are ‘none’, ‘sum’, and ‘max’. Default: ‘sum’. tolpivot: relative pivot tolerance for threshold partial pivoting with row interchanges. Default: 0.1 tolsympivot: if diagonal pivoting is attempted, this parameter is used to control when the diagonal is selected in a given pivot column. Default: 0.0 -
LU
¶ An
umfpack_context
object encapsulating the factorization.
-
L
¶ The L factor of the input matrix.
-
U
¶ The U factor of the input matrix.
-
P
¶ The row permutation used for the factorization.
-
Q
¶ The column permutation used for the factorization.
-
R
¶ The row scaling used during the factorization. See the documentation of
fetch_factors()
.
-
factorizationTime
¶ The CPU time to perform the factorization.
-
solutionTime
¶ The CPU time to perform the forward and backward sweeps.
-
do_recip
¶ Nature of the row scaling. See
fetch_factors()
.
-
lnz
¶ The number of nonzero elements in the factor L.
-
unz
¶ The number of nonzero elements in the factor U from which the diagonal was removed.
-
nz_udiag
¶ The number of nonzero elements on the diagonal of the factor U.
-
fetch_factors
()¶ Retrieve the L and U factors of the input matrix along with the permutation matrices P and Q and the row scaling matrix R such that
\mathbf{P R A Q} = \mathbf{L U}.The matrices P, R and Q are stored as Numpy arrays. L and U are stored as PysparseMatrix instances and are lower triangular and upper triangular, respectively.
R is a row-scaling diagonal matrix such that
- the i-th row of A has been divided by R[i] if
do_recip = True
, - the i-th row of A has been multiplied by R[i] if
do_recip = False
.
- the i-th row of A has been divided by R[i] if
-
fetch_lunz
()¶ Retrieve the number of nonzeros in the factors. The results are stored in the members
lnz
,unz
andnz_udiag
of the class instance.
-
solve
(rhs, **kwargs)¶ Solve the linear system
A x = rhs
. The result is placed in thesol
member of the class instance.Parameters: rhs: a Numpy vector of appropriate dimension. Keywords: method: specifies the type of system being solved:
"UMFPACK_A"
\mathbf{A} x = b (default) "UMFPACK_At"
\mathbf{A}^T x = b "UMFPACK_Pt_L"
\mathbf{P}^T \mathbf{L} x = b "UMFPACK_L"
\mathbf{L} x = b "UMFPACK_Lt_P"
\mathbf{L}^T \mathbf{P} x = b "UMFPACK_Lt"
\mathbf{L}^T x = b "UMFPACK_U_Qt"
\mathbf{U} \mathbf{Q}^T x = b "UMFPACK_U"
\mathbf{U} x = b "UMFPACK_Q_Ut"
\mathbf{Q} \mathbf{U}^T x = b "UMFPACK_Ut"
\mathbf{U}^T x = b irsteps: number of iterative refinement steps to attempt. Default: 2
-
Example: The 2D Poisson System with UMFPACK¶
The solution of a 2D Poisson system with PysparseUmfpackSolver
may look
like this:
from pysparse.tools.poisson_vec import poisson2d_sym_blk_vec
from numpy import ones
from numpy.linalg import norm
n = 200
A = PysparseMatrix( matrix=poisson2d_sym_blk_vec(n) )
x_exact = ones(n*n)/n
b = A * x_exact
LU = PysparseUmfpackSolver(A)
LU.solve(b)
print 'Factorization time: ', LU.factorizationTime
print 'Solution time: ', LU.solutionTime
print 'Error: ', norm(LU.sol - x_exact)/norm(x_exact)
This script produces the output:
Factorization time: 0.520043
Solution time: 0.031086
Error: 1.10998989668e-15