Table Of Contents

Previous topic

Iterative Solvers

Next topic

Eigenvalue Solver

This Page

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:

  1. natural ordering
  2. MMD applied to the structure of \mathbf{A}^T\mathbf{A}
  3. MMD applied to the structure of \mathbf{A}^T + \mathbf{A}
  4. 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

The pysparseSuperLU Module: A Higher-Level SuperLU Interface

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

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