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:

0. natural ordering
1. MMD applied to the structure of \(\mathbf{A}^T\mathbf{A}\)
2. MMD applied to the structure of \(\mathbf{A}^T + \mathbf{A}\)
3. COLAMD, approximate minimum degree column ordering

(default: 2).

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.

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

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 of A (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 to A.nnz. To obtain the number of nonzero element in the factors, see lunz().

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 the pysparseUmfpackSolver 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.

class pysparse.direct.directSolver.PysparseDirectSolver(A, **kwargs)

PysparseDirectSolver is a generic class and should be subclassed.

solve(b, **kwargs)

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 and diag_pivot_thresh must be small, e.g., 0.0 or 0.1. Since the value of diag_pivot_thresh is up to the user, setting symmetric to True does not automatically set permc_spec and diag_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.

0. natural ordering,
1. MMD applied to the structure of \(\mathbf{A}^T \mathbf{A}\)
2. MMD applied to the structure of \(\mathbf{A}^T + \mathbf{A}\)
3. COLAMD.

Default: 2.

LU

A superlu_context object encapsulating the factorization.

sol

The solution of the linear system after a call to solve().

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.

fetch_lunz()

Retrieve the number of nonzeros in the factors L and U together. The result is stored in the member lunz of the class instance.

solve(rhs, transpose=False)

Solve the linear system A x = rhs, where A is the input matrix and rhs is a Numpy vector of appropriate dimension. The result is placed in the sol member of the class instance.

If the optional argument transpose is True, the transpose system A^T x = rhs is solved.

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.

sol

The solution of the linear system after a call to solve().

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.

fetch_lunz()

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(rhs, **kwargs)

Solve the linear system A x = rhs. The result is placed in the sol 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