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 LUfactorisation of a general nonsymmetric sparse matrix with partial pivoting.
The superlu module exports a single function, called factorize.
The factorize function computes an LUfactorisation of the matrix A.
Parameters : 



Keywords : 


Return type:  an object of type superlu_context. This object encapsulates the L and U factors of A (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.
An abstract encapsulation of the LU factorization of a matrix by SuperLU.
A 2tuple describing the dimension of the matrix factorized. It is equal to A.shape.
The nnz attribute holds the total number of nonzero entries stored in both the L and U factors.
The solve method accepts two rank1 NumPy arrays b and x of appropriate size and assigns the solution of the linear system \(\mathbf{A}\mathbf{x} = \mathbf{b}\) to x. If the optional parameter trans is set to the string 'T', the transposed system \(\mathbf{A}^T\mathbf{x} = \mathbf{b}\) is solved instead.
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 LowLevel 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 LowLevel Sparse Matrix Types.
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.0e3):
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, 1e12, 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 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 nonsymmetric 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.
The factorize function computes an LUfactorisation of the matrix A.
Parameters : 



Keywords : 


Return type:  an object of type umfpack_context. This object encapsulates the L and U factors of A (see below). 
An abstract encapsulation of the LU factorization of a matrix by UMFPACK.
A 2tuple describing the dimension of the matrix factorized. It is equal to A.shape.
The nnz attribute holds the total number of nonzero entries stored in the input matrix. It is equal to A.nnz. To obtain the number of nonzero element in the factors, see lunz().
Parameters : 


Return the factors, permutation and scaling information. See the documentation of the pysparseUmfpackSolver below for more details.
Return the number of nonzeros in factors, i.e., in \(\mathbf{L} + \mathbf{U}\).
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 LowLevel 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 LowLevel Sparse Matrix Types.
This section anticipates on HigherLevel Sparse Matrix Classes and shows usage of higherlevel interfaces to the LU factorization packages.
A framework for solving sparse linear systems of equations using a direct factorization.
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. 720755, 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. 915952, 1999. 
[LD03]  X. S. Li and J. W. Demmel, SuperLU_DIST: A Scalable DistributedMemory Sparse Direct Solver for Unsymmetric Linear Systems, ACM Transactions on Mathematical Software 29(2), pp. 110140, 2003. 
[SLU]  http://crd.lbl.gov/~xiaoye/SuperLU 
Bases: pysparse.direct.directSolver.PysparseDirectSolver
PysparseSuperLUSolver is a wrapper class around the SuperLu library for the factorization of fullrank nbym matrices. Only matrices with real coefficients are currently supported.
Parameters : 



Keywords : 

A superlu_context object encapsulating the factorization.
The CPU time to perform the factorization.
The CPU time to perform the forward and backward sweeps.
The number of nonzero elements in the factors L and U together after a call to fetch_lunz().
Not yet available.
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.099685128150953e14
Note that this example uses the vectorized Poisson constructors of LowLevel Sparse Matrix Types.
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 preordering strategy for the unsymmetricpattern multifrontal method, ACM Transactions on Mathematical Software, 30(2), pp. 165195, 2004. 
[D04b]  T. A. Davis, Algorithm 832: UMFPACK, an unsymmetricpattern multifrontal method, ACM Transactions on Mathematical Software, 30(2), pp. 196199, 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. 119, 1999. 
[DD97]  T. A. Davis and I. S. Duff, An unsymmetricpattern multifrontal method for sparse LU factorization, SIAM Journal on Matrix Analysis and Applications, 18(1), pp. 140158, 1997. 
[UMF]  http://www.cise.ufl.edu/research/sparse/umfpack 
Bases: pysparse.direct.directSolver.PysparseDirectSolver
PysparseUmfpackSolver is a wrapper class around the UMFPACK library for the factorization of fullrank nbym matrices. Only matrices with real coefficients are currently supported.
Parameters : 



Keywords : 

An umfpack_context object encapsulating the factorization.
The L factor of the input matrix.
The U factor of the input matrix.
The row permutation used for the factorization.
The column permutation used for the factorization.
The row scaling used during the factorization. See the documentation of fetch_factors().
The CPU time to perform the factorization.
The CPU time to perform the forward and backward sweeps.
Nature of the row scaling. See fetch_factors().
The number of nonzero elements in the factor L.
The number of nonzero elements in the factor U from which the diagonal was removed.
The number of nonzero elements on the diagonal of the factor U.
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
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 rowscaling diagonal matrix such that
Retrieve the number of nonzeros in the factors. The results are stored in the members lnz, unz and nz_udiag of the class instance.
Solve the linear system A x = rhs. The result is placed in the sol member of the class instance.
Parameters : 



Keywords : 

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.10998989668e15