.. Description of the itsolvers module
.. module:: itsolvers
.. _itsolvers-page:

=================
Iterative Solvers
=================


The :mod:`itsolvers` Module
===========================

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 user-defined iterative solvers should also follow these conventions, since
other PySparse modules rely on them (e.g. the ``jdsym`` module;
see :ref:`jdsym-page`).

Let's illustrate the calling conventions, using the PCG method.

.. function:: info, iter, relres = pcg(A, b, x, tol, maxit[, K])

   Solve a linear system ``A x = b`` with the preconditioned conjugate gradient
   algorithm. Overwrite ``x`` with the best estimate of the solution found.

   :parameters:

        :A: The coefficient matrix of the linear system of equations. ``A``
             must provide the ``shape`` attribute and the ``matvec`` and
             ``matvec_transp`` methods for multiplying with a vector.
        :b: The right-hand-side of the linear system as a rank-1 NumPy array.
        :x: A rank-1 NumPy array. Upon entry, ``x`` holds the initial guess. On
             exit, ``x`` holds an approximate solution of the linear system.
        :tol: A float value representing the requested error tolerance. The
             exact meaning of this parameter depends on the actual iterative
             solver.
        :maxit: An integer that specifies the maximum number of iterations to
             be executed.
        :K: A preconditioner object that supplies the ``shape`` attribute and
             the ``precon`` method.

   :returns:

        :info: an integer that contains the exit status of the iterative solver.
           ``info >= 0`` indicates, that ``x`` holds an acceptable solution, and
           ``info < 0`` indicates an error condition. ``info`` has one of the
           following values:

           +2. iteration converged, residual is as small as seems reasonable on this machine,

           +1. iteration converged, ``b = 0``, so the exact solution is ``x = 0``.

           +0. iteration converged, relative error appears to be less than ``tol``.

           -1. iteration did not converge, maximum number of iterations was reached.

           -2. iteration did not converge, the system involving the preconditioner was ill-conditioned.

           -3. iteration did not converge, an inner product of the form :math:`\mathbf{x}^T \mathbf{K}^{-1} \mathbf{x}` was not positive, so the preconditioning matrix :math:`\mathbf{K}` does not appear to be positive definite.

           -4. iteration did not converge, the matrix ``A`` appears to be very ill-conditioned.

           -5. iteration did not converge, the method stagnated.

           -6. iteration did not converge, a scalar quantity became too small or too large to continue computing.

        :iter: the number of iterations performed.

        :relres: the relative residual at the approximate solution computed by
           the iterative method. What this actually is depends on the actual
           iterative method used.

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.


``itsolvers`` Module Functions
------------------------------

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.

.. function:: info, iter, relres = pcg(A, b, x, tol, maxit[, K])

   Implements the Preconditioned Conjugate Gradient method.

.. function:: info, iter, relres = minres(A, b, x, tol, maxit[, K])

   Implements the MINRES method.

.. function:: info, iter, relres = qmrs(A, b, x, tol, maxit[, K])

   Implements the QMRS method.

.. function:: info, iter, relres = cgs(A, b, x, tol, maxit[, K])

   Implements the CGS method.


Example: Solving the Poisson System
-----------------------------------

Let's solve the Poisson system

.. math:: \mathbf{L} \mathbf{x} = \mathbf{1},
   :label: eq:python:1
  
using the PCG method. :math:`\mathbf{L}` is the 2D Poisson matrix, introduced in
:ref:`spmatrix-page`, and :math:`\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, 1e-12, 2000)

The code makes use of the Python function ``poisson2d_sym_blk``,
which was defined in :ref:`spmatrix-page`.

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, 1e-12, 2000, Kssor)

The Matlab solution (without preconditioner) may look as follows:

.. code-block:: matlab

   n = 300;
   L = poisson2d_kron(n);
   [x,flag,relres,iter] = pcg(L, ones(n*n,1), 1e-12, 2000, ...
                              [], [], zeros(n*n,1));


Performance comparison with Matlab and native C
-----------------------------------------------

.. 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 :eq:`eq:python: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 matrix-vector 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.

.. _python-vs-matlab-vs-c:
.. table:: Performance comparison of Python, Matlab and native C implementations
           to solve the linear system :eq:`eq:python:1` without
           preconditioning. The execution times are given in seconds. *Assembly*
           is the time for constructing the matrix (or reading it from a file in
           the case of native C).  *Solve* is the time spent in the PCG
           solver. *Total* is the sum of *Assembly* and *Solve*. Matlab version
           6.0 Release 12 was used for these timings.

   +----------+-------+----------+---------+--------+
   | 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 :eq:`eq:python:1`. Matlab is not
only slower when building the matrix, also the matrix-vector 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.
