bicgstab.hpp 3.72 KB
Newer Older
1
// Software License for MTL
2
//
3
4
5
6
7
// Copyright (c) 2007 The Trustees of Indiana University.
//               2008 Dresden University of Technology and the Trustees of Indiana University.
//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
8
//
9
// This file is part of the Matrix Template Library
10
//
11
12
13
14
15
16
17
18
19
20
21
22
23
24
// See also license.mtl.txt in the distribution.

#ifndef ITL_BICGSTAB_INCLUDE
#define ITL_BICGSTAB_INCLUDE

#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/mtl/utility/exception.hpp>
#include <boost/numeric/itl/utility/exception.hpp>
#include <boost/numeric/mtl/operation/resource.hpp>
#include <boost/numeric/mtl/interface/vpt.hpp>

namespace itl {

///  Bi-Conjugate Gradient Stabilized
25
template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB,
26
	   class Preconditioner, class Iteration >
27
int bicgstab(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
28
29
30
31
32
33
34
	     const Preconditioner& M, Iteration& iter)
{
  typedef typename mtl::Collection<HilbertSpaceX>::value_type Scalar;
  typedef HilbertSpaceX                                       Vector;
  mtl::vampir_trace<7004> tracer;

  Scalar     rho_1(0), rho_2(0), alpha(0), beta(0), gamma, omega(0);
35
  Vector     p(resource(x)), phat(resource(x)), s(resource(x)), shat(resource(x)),
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
             t(resource(x)), v(resource(x)), r(resource(x)), rtilde(resource(x));

  r = b - A * x;
  rtilde = r;

  while (! iter.finished(r)) {
    ++iter;
    rho_1 = dot(rtilde, r);
    MTL_THROW_IF(rho_1 == 0.0, unexpected_orthogonality());

    if (iter.first())
      p = r;
    else {
      MTL_THROW_IF(omega == 0.0, unexpected_orthogonality());
      beta = (rho_1 / rho_2) * (alpha / omega);
      p = r + beta * (p - omega * v);
    }
    phat = solve(M, p);
    v = A * phat;

    gamma = dot(rtilde, v);
    MTL_THROW_IF(gamma == 0.0, unexpected_orthogonality());

    alpha = rho_1 / gamma;
    s = r - alpha * v;
61

62
63
64
65
66
67
68
    if (iter.finished(s)) {
      x += alpha * phat;
      break;
    }
    shat = solve(M, s);
    t = A * shat;
    omega = dot(t, s) / dot(t, t);
69

70
71
    x += omega * shat + alpha * phat;
    r = s - omega * t;
72
73

    rho_2 = rho_1;
74
75
76
77
78
  }
  return iter;
}

/// Solver class for BiCGStab method; right preconditioner ignored (prints warning if not identity)
79
template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
80
81
82
83
84
	   typename RightPreconditioner= pc::identity<LinearOperator> >
class bicgstab_solver
{
  public:
    /// Construct solver from a linear operator; generate (left) preconditioner from it
85
    explicit bicgstab_solver(const LinearOperator& A) : A(A), L(A)
86
87
88
89
90
91
    {
	if (!pc::static_is_identity<RightPreconditioner>::value)
	    std::cerr << "Right Preconditioner ignored!" << std::endl;
    }

    /// Construct solver from a linear operator and (left) preconditioner
92
    bicgstab_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L)
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
    {
	if (!pc::static_is_identity<RightPreconditioner>::value)
	    std::cerr << "Right Preconditioner ignored!" << std::endl;
    }

    /// Solve linear system approximately as specified by \p iter
    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
    {
	return bicgstab(A, x, b, L, iter);
    }

    /// Perform one BiCGStab iteration on linear system
    template < typename HilbertSpaceB, typename HilbertSpaceX >
    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
    {
	itl::basic_iteration<double> iter(x, 1, 0, 0);
	return solve(b, x, iter);
    }
112

113
114
115
116
117
118
119
120
  private:
    const LinearOperator& A;
    Preconditioner        L;
};

} // namespace itl

#endif // ITL_BICGSTAB_INCLUDE