Copy.hpp 1.58 KB
Newer Older
1
2
#pragma once

3
4
#include <algorithm>

5
6
#include <amdis/linear_algebra/mtl/BlockMTLVector.hpp>
#include <amdis/linear_algebra/mtl/BlockMTLMatrix.hpp>
7

8
9
#include <amdis/linear_algebra/mtl/SystemVector.hpp>
#include <amdis/linear_algebra/mtl/SystemMatrix.hpp>
10

11
#include <amdis/linear_algebra/mtl/Mapper.hpp>
12

13
namespace AMDiS
14
15
{
  // requires MatrixType to have an inserter
16
  template <class Value, std::size_t _N, std::size_t _M, class TargetType>
17
18
19
20
21
22
  void copy(BlockMTLMatrix<Value, _N, _M> const& source, TargetType& target)
  {
    using Mapper = std::conditional_t<(_N == _M), BlockMapper, RectangularMapper>;
    using value_type = typename TargetType::value_type;
    using BaseInserter = mtl::mat::inserter<TargetType, mtl::operations::update_plus<value_type>>;
    using MappedInserter = typename Mapper::template Inserter<BaseInserter>::type;
23
24
25
26
27
28

    std::size_t nnz = 0;
    for (std::size_t rb = 0; rb < _N; ++rb)
      for (std::size_t cb = 0; cb < _M; ++cb)
        nnz += source[rb][cb].nnz();

29
30
31
    {
      target.change_dim(num_rows(source), num_cols(source));
      set_to_zero(target);
32

33
34
      Mapper mapper(source);
      MappedInserter ins(target, mapper, int(1.2 * nnz / target.dim1()));
35
36

      for (std::size_t rb = 0; rb < _N; ++rb) {
37
        mapper.setRow(rb);
38
39
        for (std::size_t cb = 0; cb < _M; ++cb) {
	        mapper.setCol(cb);
40
41
42
43
44
          ins << source[rb][cb];
        }
      }
    }
  }
45

46
47
48
49
50
  template <class FeSpaces, class Value, class TargetType>
  void copy(SystemMatrix<FeSpaces, Value> const& source, TargetType& target)
  {
    copy(source.getMatrix(), target);
  }
51

52
} // end namespace AMDiS