masslumping.hpp 3.87 KB
Newer Older
1
#pragma once
2

3
#include <boost/numeric/linear_algebra/inverse.hpp>
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95

#include <boost/numeric/mtl/vector/dense_vector.hpp>
#include <boost/numeric/mtl/utility/exception.hpp>
#include <boost/numeric/mtl/utility/tag.hpp>
#include <boost/numeric/mtl/detail/base_cursor.hpp>
#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/mtl/operation/resource.hpp>
#include <boost/numeric/mtl/interface/vpt.hpp>
#include <boost/numeric/itl/pc/solver.hpp>

namespace itl
{
  namespace pc
  {
    /// Diagonal Preconditioner
    template <typename Matrix, typename Value= typename mtl::Collection<Matrix>::value_type>
    class masslumping
    {
    public:
      typedef Value                                         value_type;
      typedef typename mtl::Collection<Matrix>::size_type   size_type;
      typedef masslumping                                   self;

      /// Constructor takes matrix reference
      explicit masslumping(const Matrix& A) : inv_diag(num_rows(A))
      {
        MTL_THROW_IF(num_rows(A) != num_cols(A), mtl::matrix_not_square());
        using namespace mtl;
        using namespace mtl::tag;
        using mtl::traits::range_generator;
        using math::reciprocal;

        typedef typename range_generator<tag::row, Matrix>::type c_type;
        typedef typename range_generator<tag::nz, c_type>::type  ic_type;

        typename mtl::traits::row<Matrix>::type row(A);
        typename mtl::traits::const_value<Matrix>::type value(A);

        value_type row_sum;
        c_type cursor(begin<tag::row>(A));
        for (c_type cend(end<tag::row>(A)); cursor != cend; ++cursor)
        {
          row_sum = math::zero(value_type());
          ic_type icursor(begin<tag::nz>(cursor));
          size_type r = row(icursor);

          for (ic_type icend(end<tag::nz>(cursor)); icursor != icend; ++icursor)
            row_sum += value(*icursor);

          inv_diag[r]= reciprocal(row_sum);
        }
      }

      /// Member function solve, better use free function solve
      template <typename Vector>
      Vector solve(const Vector& x) const
      {
        Vector y(resource(x));
        solve(x, y);
        return y;
      }

      template <typename VectorIn, typename VectorOut>
      void solve(const VectorIn& x, VectorOut& y) const
      {
        mtl::vampir_trace<5051> tracer;
        y.checked_change_resource(x);
        MTL_THROW_IF(size(x) != size(inv_diag), mtl::incompatible_size());
        for (size_type i= 0; i < size(inv_diag); ++i)
          y[i]= inv_diag[i] * x[i];
      }

      /// Member function for solving adjoint problem, better use free function adjoint_solve
      template <typename Vector>
      Vector adjoint_solve(const Vector& x) const
      {
        Vector y(resource(x));
        adjoint_solve(x, y);
        return y;
      }

      template <typename VectorIn, typename VectorOut>
      void adjoint_solve(const VectorIn& x, VectorOut& y) const
      {
        using mtl::conj;
        y.checked_change_resource(x);
        MTL_THROW_IF(size(x) != size(inv_diag), mtl::incompatible_size());
        for (size_type i= 0; i < size(inv_diag); ++i)
          y[i]= conj(inv_diag[i]) * x[i];
      }

    protected:
96
      mtl::vec::dense_vector<value_type>    inv_diag;
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
    };

    /// Solve approximately a sparse system in terms of inverse lumped diagonal
    template <typename Matrix, typename Vector>
    solver<masslumping<Matrix>, Vector, false>
    inline solve(const masslumping<Matrix>& P, const Vector& x)
    {
      return solver<masslumping<Matrix>, Vector, false>(P, x);
    }

    /// Solve approximately the adjoint of a sparse system in terms of inverse lumped diagonal
    template <typename Matrix, typename Vector>
    solver<masslumping<Matrix>, Vector, true>
    inline adjoint_solve(const masslumping<Matrix>& P, const Vector& x)
    {
      return solver<masslumping<Matrix>, Vector, true>(P, x);
    }

115
116
  } // end namespace pc
} // end namespace itl