block_diagonal.hpp 2.33 KB
Newer Older
1
2
3
4
#pragma once

#include <boost/numeric/itl/pc/diagonal.hpp>

5
#include <amdis/linear_algebra/mtl/BlockMTLMatrix.hpp>
6
7
8
9
10
11

namespace itl
{
  namespace pc
  {
    /// Diagonal Preconditioner
12
    template <class MTLMatrix, std::size_t _N, std::size_t _M, class Value>
13
14
15
16
    class diagonal<AMDiS::BlockMTLMatrix<MTLMatrix, _N, _M>, Value>
    {
      using Self       = diagonal;
      using MatrixType = AMDiS::BlockMTLMatrix<MTLMatrix, _N, _M>;
17

18
19
20
21
22
23
24
25
26
    public:
      using value_type = Value;
      using size_type  = typename mtl::Collection<MatrixType>::size_type;

      /// Constructor takes matrix reference
      explicit diagonal(const MatrixType& A) : inv_diag(num_rows(A))
      {
        MTL_THROW_IF(num_rows(A) != num_cols(A), mtl::matrix_not_square());
        using math::reciprocal;
27

28
29
30
        std::array<mtl::irange, _N> r_rows;
        A.getRowRanges(r_rows);

31
        for (std::size_t i = 0; i < _N; i++)
32
          inv_diag[r_rows[i]] = mtl::mat::diagonal(A[i][i]);
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

        for (size_type i = 0; i < num_rows(A); ++i)
          inv_diag[i]= reciprocal(inv_diag[i]);
      }

      /// 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
      {
        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:
76
      mtl::vec::dense_vector<value_type>    inv_diag;
77
78
79
80
81
    };


  }
} // namespace itl::pc