unit_filter_blocked.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/base_header.hpp>
#include <kernel/util/math.hpp>
#include <kernel/lafem/sparse_matrix_bcsr.hpp>
#include <kernel/lafem/dense_vector_blocked.hpp>
#include <kernel/lafem/sparse_vector_blocked.hpp>
#include <kernel/lafem/arch/unit_filter_blocked.hpp>
 
namespace FEAT
{
  namespace LAFEM
  {
    template<
      typename DT_,
      typename IT_,
      int BlockSize_>
    class UnitFilterBlocked
    {
    public:
      typedef DT_ DataType;
      typedef IT_ IndexType;
      static constexpr int BlockSize = BlockSize_;
      typedef Tiny::Vector<DataType, BlockSize> ValueType;
      typedef DenseVectorBlocked<DataType, IndexType, BlockSize> VectorType;
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      using FilterType = UnitFilterBlocked<DT2_, IT2_, BlockSize_>;
 
      template <typename DataType2_, typename IndexType2_>
      using FilterTypeByDI = FilterType<DataType2_, IndexType2_>;
 
      static constexpr bool is_global = false;
      static constexpr bool is_local = true;
 
      static_assert(BlockSize > 1, "BlockSize has to be >= 2 in UnitFilterBlocked!");
 
    private:
      SparseVectorBlocked<DataType, IndexType, BlockSize> _sv;
 
      bool _ignore_nans;
 
    public:
      UnitFilterBlocked() :
        _sv(),
        _ignore_nans(false)
      {
      }
 
      explicit UnitFilterBlocked(Index num_entries, bool ignore_nans = false) :
        _sv(num_entries),
        _ignore_nans(ignore_nans)
      {
      }
 
      explicit UnitFilterBlocked(Index size_in,
                                 DenseVectorBlocked<DT_, IT_, BlockSize_> & values,
                                 DenseVector<IT_, IT_> & indices, bool ignore_nans = false) :
        _sv(size_in, values, indices),
        _ignore_nans(ignore_nans)
      {
        XASSERTM(values.size() == indices.size(), "Vector size mismatch!");
      }
 
      UnitFilterBlocked(UnitFilterBlocked && other) :
        _sv(std::move(other._sv)),
        _ignore_nans(other._ignore_nans)
      {
      }
 
      UnitFilterBlocked & operator=(UnitFilterBlocked && other)
      {
        if(this != &other)
        {
          _sv = std::forward<decltype(other._sv)>(other._sv);
          _ignore_nans = other._ignore_nans;
        }
        return *this;
      }
 
      virtual ~UnitFilterBlocked()
      {
      }
 
      UnitFilterBlocked clone(CloneMode clone_mode = CloneMode::Deep) const
      {
        UnitFilterBlocked other;
        other.clone(*this, clone_mode);
        return other;
      }
 
      void clone(const UnitFilterBlocked & other, CloneMode clone_mode = CloneMode::Deep)
      {
        _sv.clone(other.get_filter_vector(), clone_mode);
        _ignore_nans = other.get_ignore_nans();
      }
 
      template<typename DT2_, typename IT2_, int BS_>
      void convert(const UnitFilterBlocked<DT2_, IT2_, BS_>& other)
      {
        _sv.convert(other.get_filter_vector());
        _ignore_nans = other.get_ignore_nans();
      }
 
      void clear()
      {
        _sv.clear();
      }
 
      std::size_t bytes() const
      {
        return _sv.bytes();
      }
 
      SparseVectorBlocked<DT_, IT_, BlockSize>& get_filter_vector()
      {
        return _sv;
      }
      const SparseVectorBlocked<DT_, IT_, BlockSize>& get_filter_vector() const
      {
        return _sv;
      }
 
      void set_ignore_nans(bool ignore_nans)
      {
        _ignore_nans = ignore_nans;
      }
 
      bool get_ignore_nans() const
      {
        return _ignore_nans;
      }
 
      void add(IndexType idx, ValueType val)
      {
        _sv(idx, val);
      }
 
      Index size() const
      {
        return _sv.size();
      }
 
      Index used_elements() const
      {
        return _sv.used_elements();
      }
 
      IT_* get_indices()
      {
        return _sv.indices();
      }
 
      const IT_* get_indices() const
      {
        return _sv.indices();
      }
 
      ValueType* get_values()
      {
        return _sv.elements();
      }
 
      const ValueType* get_values() const
      {
        return _sv.elements();
      }
 
      void permute(Adjacency::Permutation & perm)
      {
        _sv.permute(perm);
      }
 
#ifdef DOXYGEN
      // The following documentation block is visible to Doxygen only. The actual implementation is matrix type
      // specific and provided below.
 
      void filter_mat(MatrixType & matrix) const
      {
      }
 
      void filter_offdiag_row_mat(MatrixType & matrix) const
      {
      }
 
      void filter_offdiag_col_mat(MatrixType & matrix) const
      {
      }
 
#endif
      template<int BlockWidth_>
      void filter_mat(SparseMatrixBCSR<DT_, IT_, BlockSize_, BlockWidth_> & matrix) const
      {
        if(_sv.used_elements() == Index(0))
          return;
 
        XASSERTM(_sv.size() == matrix.rows(), "Matrix size does not match!");
 
        Arch::UnitFilterBlocked::filter_unit_mat
          (matrix.template val<LAFEM::Perspective::pod>(), matrix.row_ptr(), matrix.col_ind(), BlockSize_, BlockWidth_,
          _sv.template elements<LAFEM::Perspective::pod>(), _sv.indices(), _sv.used_elements(), _ignore_nans);
      }
 
      template<int BlockWidth_>
      void filter_offdiag_row_mat(SparseMatrixBCSR<DT_, IT_, BlockSize_, BlockWidth_> & matrix) const
      {
        if(_sv.used_elements() == Index(0))
          return;
 
        XASSERTM(_sv.size() == matrix.rows(), "Matrix size does not match!");
 
        Arch::UnitFilterBlocked::filter_offdiag_row_mat
          (matrix.template val<LAFEM::Perspective::pod>(), matrix.row_ptr(), BlockSize_, BlockWidth_,
          _sv.template elements<LAFEM::Perspective::pod>(), _sv.indices(), _sv.used_elements(), _ignore_nans);
      }
 
      template<int BlockWidth_>
      void filter_offdiag_col_mat(SparseMatrixBCSR<DT_, IT_, BlockSize_, BlockWidth_> &) const
      {
        // nothing to do here
      }
 
 
      template<int BlockWidth_>
      void filter_weak_matrix_rows(SparseMatrixBCSR<DT_, IT_, BlockSize_, BlockWidth_> & matrix_a,
        const SparseMatrixBCSR<DT_, IT_, BlockSize_, BlockWidth_>& matrix_m) const
      {
        if(_sv.used_elements() == Index(0))
          return;
 
        if(matrix_a.val() == matrix_m.val())
        {
          XABORTM("Matrices are not allowed to hold the same data");
        }
 
        XASSERTM(_sv.size() == matrix_a.rows(), "Matrix size does not match!");
        XASSERTM(_sv.size() == matrix_m.rows(), "Matrix size does not match!");
 
        const IndexType* row_ptr(matrix_a.row_ptr());
        const IndexType* col_idx(matrix_m.col_ind());
        XASSERTM(row_ptr == matrix_m.row_ptr(), "matrix A and M must share their layout");
        XASSERTM(col_idx == matrix_m.col_ind(), "matrix A and M must share their layout");
 
        Arch::UnitFilterBlocked::filter_weak_matrix_rows
          (matrix_a.template val<LAFEM::Perspective::pod>(), matrix_m.template val<LAFEM::Perspective::pod>(), row_ptr, BlockSize_, BlockWidth_,
          _sv.template elements<LAFEM::Perspective::pod>(), _sv.indices(), _sv.used_elements());
      }
 
 
      void filter_rhs(VectorType& vector) const
      {
        if(_sv.used_elements() == Index(0))
          return;
 
        XASSERTM(_sv.size() == vector.size(), "Vector size does not match!");
        Arch::UnitFilterBlocked::filter_rhs
          (vector.template elements<Perspective::pod>(), BlockSize_, _sv.template elements<Perspective::pod>(),
          _sv.indices(), _sv.used_elements(), _ignore_nans);
      }
 
      void filter_sol(VectorType& vector) const
      {
        // same as rhs
        filter_rhs(vector);
      }
 
      void filter_def(VectorType& vector) const
      {
        if(_sv.used_elements() == Index(0))
          return;
 
        XASSERTM(_sv.size() == vector.size(), "Vector size does not match!");
        Arch::UnitFilterBlocked::filter_def
          (vector.template elements<Perspective::pod>(), BlockSize_, _sv.template elements<Perspective::pod>(),
          _sv.indices(), _sv.used_elements(), _ignore_nans);
      }
 
      void filter_cor(VectorType& vector) const
      {
        // same as def
        filter_def(vector);
      }
    }; // class UnitFilterBlocked<...>
  } // namespace LAFEM
} // namespace FEAT