minres.hpp 1.76 KB
Newer Older
1
#pragma once
2
3
4
5
6
7
8

#include <boost/numeric/mtl/concept/collection.hpp>

namespace itl
{
  /// Minimal Residual method
  template <typename Matrix, typename Vector,
9
            typename LeftPreconditioner, typename Iteration>
10
  int minres(const Matrix& A, Vector& x, const Vector& b,
11
             const LeftPreconditioner& P, Iteration& iter)
12
  {
13
    using std::abs; using std::sqrt;
14
15
16
17
18
19
    using math::reciprocal;
    typedef typename mtl::Collection<Vector>::value_type Scalar;

    if (size(b) == 0)
      throw mtl::logic_error("empty rhs vector");

20
21
    Scalar zero= math::zero(b[0]), one= math::one(b[0]);
    Vector v0(size(x), zero), v1(b - A * x), v2(v1), z1(solve(P, v1)), z2(size(x), zero);
22
23
24
    Vector w0(size(x), zero), w1(size(x), zero), w2(size(x), zero);

    Scalar s0(zero), s1(zero), c0(one), c1(one), gamma0(one);
25
    Scalar gamma1(sqrt(abs(dot(z1, v1)))), gamma2(zero), eta(gamma1);
26
27
28
29
30
31
32
33
34
    Scalar sigma1(one), alpha0(zero), alpha1(zero), alpha2(zero), alpha3(zero);

    while (!iter.finished(abs(eta)))
    {
      z1 *= reciprocal(gamma1);
      v2 = A * z1;
      sigma1 = dot(v2, z1);
      v2 += -(sigma1 / gamma1) * v1 - (gamma1 / gamma0) * v0;

35
      z2 = solve(P, v2);
36

37
      gamma2 = sqrt(abs(dot(z2, v2)));
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
      alpha0 = c1 * sigma1 - c0 * s1 * gamma1;
      alpha1 = sqrt(alpha0 * alpha0 + gamma2 * gamma2);
      alpha2 = s1 * sigma1 + c0 * c1 * gamma1;
      alpha3 = s0 * gamma1;

      c0 = c1;
      c1 = alpha0 / alpha1;
      s0 = s1;
      s1 = gamma2 / alpha1;

      w2 = z1 - alpha3 * w0 - alpha2 * w1;
      w2 *=  reciprocal(alpha1);

      x += c1 * eta * w2;
      eta *= -s1;

      w0 = w1;
      w1 = w2;
      v0 = v1;
      v1 = v2;
      z1 = z2;

      gamma0 = gamma1;
      gamma1 = gamma2;

      ++iter;
    }

    return iter;
  }

69
} // end namespace itl