preonly.hpp 720 Bytes
Newer Older
1
#pragma once
2

3
#include <cassert>
4
5
6
7
8
9
10
11
#include <boost/numeric/mtl/concept/collection.hpp>

namespace itl
{
  /// Solver that simply applies a preconditioner to the rhs vector
  template <typename Matrix, typename Vector, typename Preconditioner, typename Iteration>
  int preonly(const Matrix& A, Vector& x, const Vector& b, const Preconditioner& P, Iteration& iter)
  {
12
    assert(size(b) != 0);
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

    typedef typename mtl::Collection<Vector>::value_type Scalar;

    // simple richardson iteration
    Vector r(b - A*x);
    Scalar res = two_norm(r);
    for (; !iter.finished(res); ++iter)
    {
      x += Vector(solve(P, r));
      r = b - A*x;
      res = two_norm(r);
    }
    return iter;
  }

28
} // end namespace itl