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;
matrix_type A(N, N);
matrix::laplacian_setup(A, size, size);
pc::ilu_0<matrix_type> P(A);
dense_vector<double> x(N, 1.0), b(N);
b= A * x; x= 0;
noisy_iteration<double> iter(b, 500, 1.e-6);
bicgstab(A, x, b, P, iter);
return 0;
}
Currently the folling solvers (in alphabetical order) are available:
- Bi-Conjugate Gradient: itl::bicg(A, x, b, L, iter);
- Bi-Conjugate Gradient Stabilized: itl::bicgstab(A, x, b, L, iter);
- Bi-Conjugate Gradient Stabilized(2): itl::bicgstab_2(A, x, b, L, iter);
- Bi-Conjugate Gradient Stabilized(ell): itl::bicgstab_ell(A, x, b, L, R, iter, ell);
- Conjugate Gradient: itl::cg(A, x, b, L, iter);
- Conjugate Gradient Squared: itl::cgs(A, x, b, L, iter);
- Generalized Minimal Residual method (without restart): itl::gmres_full(A, x, b, L, R, iter, kmax_in);
- Generalized Minimal Residual method with restart: itl::gmres(A, x, b, L, R, iter, restart);
- Induced Dimension Reduction on s dimensions (IDR(s)): itl::idr_s(A, x, b, L, R, iter, s);
- Quasi-minimal residual: itl::qmr(A, x, b, L, R, iter); and
- Transposed-free Quasi-minimal residual: itl::tfqmr(A, x, b, L, R, iter).
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:
- Diagonal inversion: itl::pc::diagonal<Matrix>;
- Incomplete LU factorization without fill-in: itl::pc::ilu_0<Matrix>; and
- Incomplete Cholesky factorization without fill-in: itl::pc::ic_0<Matrix>;
The iteration object can be chosen between:
- Basic iteration does not generate output: basic_iteration(r0, m, r, a= 0);
- Cyclic iteration prints residual information every c iteration: cyclic_iteration(r0, m, r, a= 0, c= 100, out= std::cout); and
- Noisy iteration prints residual in each iteration: noisy_iteration(r0, m, r, a= 0, out= std::cout).
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:
- The maximum number of iterations is reached (failure);
- The relative residuum reduction was achieved; or
- The absolute residuum is below a if specified.
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:
- 0th iteration is the starting residue.
- Once the input value (x) is changed you have made at least one iteration.
- Fractions of iterations are counted as whole iterations.
Return to Introduction Krylov-Subspace Methods Table of Content Proceed to Iteration