Logo MTL4

Using Predefined Linear Solvers

The following program illustrates how to solve a linear system:

#include <iostream>
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/itl/itl.hpp>

using namespace mtl;
using namespace itl;

int main(int, char**)
{
    const int size = 40, N = size * size; 
    typedef compressed2D<double>  matrix_type;

    // Set up a matrix 1,600 x 1,600 with 5-point-stencil
    matrix_type                   A(N, N);
    matrix::laplacian_setup(A, size, size);

    // Create an ILU(0) preconditioner
    pc::ilu_0<matrix_type>        P(A);
    
    // Set b such that x == 1 is solution; start with x == 0
    dense_vector<double>          x(N, 1.0), b(N);
    b= A * x; x= 0;
    
    // Termination criterion: r < 1e-6 * b or N iterations
    noisy_iteration<double>       iter(b, 500, 1.e-6);
    
    // Solve Ax == b with left preconditioner P
    bicgstab(A, x, b, P, iter);

    return 0;
}

Currently the folling solvers (in alphabetical order) are available:

All Krylov sub-space methods solve the linear system Ax = b as in the example above. A left preconditioner L is used in all methods and some methods also incorporate a right preconditioner R. The iteration object controls the termination of the iteration, see below. Some algorithms take an additional argument (ell, kmax_in, restart, s) specifying the dimension of the vector space regarding to which new search directions are orthogonalized. The difference between gmres and gmres_full is that continues until one criterion in iter holds. gmres_full computes at most kmax_in iterations (or size(x) depending on what is smaller) regardless on whether the termination criterion is reached or not. Needless to say that gmres is implemented by means of gmres_full.

As preconditioners we provide at the moment:

The iteration object can be chosen between: Mandatory arguments for the iteration objects' constructors are the initial residuum r0 (which is of course b if one starts with x = 0), the maximal number of iteration and the relative error reduction (more precisely residuum reduction). Optionally, the absolute residuum a (fourth argument) can be given as termination criterion. Thus the iterative methods are terminated when either: In the cyclic iteration, the user can specify after how many iterations the residual informations are printed, default is 100. For cyclic and noisy iterations, one can declare on which ostream the information is printed. This enables printing it into log files or for parallel computing printing only on one processor. By default the output is printed into std::out.

General assumptions on solver iterations:

Return to Introduction Krylov-Subspace Methods                                Table of Content                                Proceed to Iteration


Using Predefined Linear Solvers -- MTL 4 -- Peter Gottschling and Andrew Lumsdaine -- Gen. with rev. 7542 on 7 Apr 2011 by doxygen 1.5.9 -- © 2010 by SimuNova UG.