The itsolvers module provides a set of iterative methods for solving linear systems of equations.
The iterative methods are callable like ordinary Python functions. All these functions expect the same parameter list, and all function return values also follow a common standard.
Any userdefined iterative solvers should also follow these conventions, since other PySparse modules rely on them (e.g. the jdsym module; see Eigenvalue Solver).
Let’s illustrate the calling conventions, using the PCG method.
Solve a linear system A x = b with the preconditioned conjugate gradient algorithm. Overwrite x with the best estimate of the solution found.
Parameters : 



Returns: 

The iterative solvers may accept additional parameters, which are passed as keyword arguments.
Note that not all iterative solvers check for all above error conditions.
The module functions defined in the precon module implement various iterative methods (PCG, MINRES, QMRS and CGS). The parameters and return values conform to the conventions described above.
Implements the Preconditioned Conjugate Gradient method.
Implements the MINRES method.
Implements the QMRS method.
Implements the CGS method.
Let’s solve the Poisson system
using the PCG method. \(\mathbf{L}\) 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.precon import precon
from pysparse.itsolvers import krylov
import numpy
n = 300
L = poisson2d_sym_blk(n)
b = numpy.ones(n*n)
x = numpy.empty(n*n)
info, iter, relres = krylov.pcg(L.to_sss(), b, x, 1e12, 2000)
The code makes use of the Python function poisson2d_sym_blk, which was defined in LowLevel Sparse Matrix Types.
Incorporating e.g. a SSOR preconditioner is straightforward:
from pysparse.sparse import spmatrix
from pysparse.precon import precon
from pysparse.itsolver import krylov
import numpy
n = 300
L = poisson2d_sym_blk(n)
b = numpy.ones(n*n)
x = numpy.empty(n*n)
S = L.to_sss()
Kssor = precon.ssor(S)
info, iter, relres = krylov.pcg(S, b, x, 1e12, 2000, Kssor)
The Matlab solution (without preconditioner) may look as follows:
n = 300;
L = poisson2d_kron(n);
[x,flag,relres,iter] = pcg(L, ones(n*n,1), 1e12, 2000, ...
[], [], zeros(n*n,1));
Todo
Update the timings below.
Warning
The timings below are Roman’s old benchmarks. I don’t know on which machine they were run. We should update them.
To evaluate the performance of the Python implementation we solve the 2D Poisson system (1) using the PCG method. The Python timings are compared with results of a Matlab and a native C implementation.
The native C and the Python implementation use the same core algorithms for PCG method and the matrixvector multiplication. On the other hand, C reads the matrix from an external file instead of building it on the fly. In contrast to the Python implementation, the native C version does not suffer from the overhead generated by the runtime argument parsing and calling overhead.
Function  Size  Assembly  Solve  Total 
Python  n=100  0.03  1.12  1.15 
n=300  0.21  49.65  49.86  
n=500  0.62  299.39  300.01  
Native C  n=100  0.30  0.96  1.26 
n=300  3.14  48.38  51.52  
n=500  10.86  288.67  299.53  
Matlab  n=100  0.21  8.85  9.06 
n=300  2.05  387.26  389.31  
n=500  6.23  1905.67  1911.8 
This table shows the execution times for the Python, the Matlab and the native C implementation for solving the linear system (1). Matlab is not only slower when building the matrix, also the matrixvector multiplication seems to be implemented inefficiently. Considering Solve, the performance of Python and native C is comparable. The Python overhead is under a factor of 4.