vector_mirror.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/lafem/container.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/lafem/dense_vector_blocked.hpp>
#include <kernel/lafem/sparse_vector.hpp>
#include <kernel/lafem/sparse_vector_blocked.hpp>
#include <kernel/lafem/sparse_matrix_csr.hpp>
#include <kernel/lafem/arch/mirror.hpp>
 
namespace FEAT
{
  namespace LAFEM
  {
    template<typename DT_, typename IT_>
    class VectorMirror :
      public Container<DT_, IT_>
    {
    public:
      typedef Container<DT_, IT_> BaseClass;
      typedef DT_ DataType;
      typedef IT_ IndexType;
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      using MirrorType = VectorMirror<DT2_, IT2_>;
 
      template <typename DataType2_, typename IndexType2_>
      using MirrorTypeByDI = MirrorType<DataType2_, IndexType2_>;
 
      typedef IT_* ImageIterator;
 
    public:
      VectorMirror() :
        BaseClass(0)
      {
        this->_scalar_index.push_back(0); // number of indices
      }
 
      explicit VectorMirror(Index size_in, Index num_idx) :
        BaseClass(size_in)
      {
        this->_scalar_index.push_back(num_idx); // number of indices
        if(num_idx > Index(0))
        {
          this->_indices.push_back(MemoryPool::template allocate_memory<IndexType>(num_idx));
          this->_indices_size.push_back(num_idx);
        }
      }
 
      VectorMirror(VectorMirror&& other) :
        BaseClass(std::forward<BaseClass>(other))
      {
      }
 
      VectorMirror& operator=(VectorMirror&& other)
      {
        if(this == &other)
          return *this;
 
        this->move(std::forward<BaseClass>(other));
 
        return *this;
      }
 
      static VectorMirror make_identity(Index size_in)
      {
        VectorMirror mir(size_in, size_in);
        IndexType* idx = mir.indices();
        for(Index i(0); i < size_in; ++i)
          idx[i] = IndexType(i);
        return mir;
      }
 
      static VectorMirror make_empty(const DenseVector<DT_, IT_>& tmpl_vec)
      {
        return VectorMirror(tmpl_vec.size(), 0u);
      }
 
      template<int block_size_>
      static VectorMirror make_empty(const DenseVectorBlocked<DT_, IT_, block_size_>& tmpl_vec)
      {
        return VectorMirror(tmpl_vec.template size<LAFEM::Perspective::native>(), 0u);
      }
 
      VectorMirror clone(CloneMode clone_mode = CloneMode::Weak) const
      {
        VectorMirror t;
        t.clone(*this, clone_mode);
        return t;
      }
 
      template<typename DT2_, typename IT2_>
      void clone(const VectorMirror<DT2_, IT2_> & other, CloneMode clone_mode = CloneMode::Weak)
      {
        Container<DT_, IT_>::clone(other, clone_mode);
      }
 
      template<typename DT2_, typename IT2_>
      void convert(const VectorMirror<DT2_, IT2_>& other)
      {
        this->assign(other);
      }
 
      Index num_indices() const
      {
        return this->_scalar_index.empty() ? Index(0) : this->_scalar_index.at(1);
      }
 
      bool empty() const
      {
        return (this->num_indices() == Index(0));
      }
 
      IT_* indices()
      {
        return this->_indices.empty() ? nullptr : this->_indices.at(0);
      }
 
      const IT_* indices() const
      {
        return this->_indices.empty() ? nullptr : this->_indices.at(0);
      }
 
      template<typename DT2_, typename IT2_>
      Index buffer_size(const DenseVector<DT2_, IT2_>& DOXY(vector)) const
      {
        return num_indices();
      }
 
      template<typename DT2_, typename IT2_, int block_size_>
      Index buffer_size(const DenseVectorBlocked<DT2_, IT2_, block_size_>& DOXY(vector)) const
      {
        return num_indices()*Index(block_size_);
      }
 
      template<typename DT2_, typename IT2_>
      Index buffer_size(const SparseVector<DT2_, IT2_>& DOXY(vector)) const
      {
        return num_indices();
      }
 
      template<typename DT2_, typename IT2_, int block_size_>
      Index buffer_size(const SparseVectorBlocked<DT2_, IT2_, block_size_>& DOXY(vector)) const
      {
        return num_indices()*Index(block_size_);
      }
 
      template<typename Vector_>
      DenseVector<DataType, IndexType> create_buffer(const Vector_& vector) const
      {
        return DenseVector<DataType, IndexType>(buffer_size(vector));
      }
 
      void gather(
        LAFEM::DenseVector<DataType, IndexType>& buffer,
        const LAFEM::DenseVector<DataType, IndexType>& vector,
        const Index buffer_offset = Index(0)) const
      {
        XASSERT(buffer_offset + this->num_indices() <= buffer.size());
        XASSERTM(this->size() == vector.size(), "size mismatch between mirror and vector");
 
        if(this->empty())
          return;
 
        LAFEM::Arch::Mirror::gather_dv(
          buffer_offset, this->num_indices(), this->indices(), buffer.elements(), vector.elements());
      }
 
      void scatter_axpy(
        LAFEM::DenseVector<DataType, IndexType>& vector,
        const LAFEM::DenseVector<DataType, IndexType>& buffer,
        const DataType alpha = DataType(1),
        const Index buffer_offset = Index(0)) const
      {
        XASSERT(buffer_offset + this->num_indices() <= buffer.size());
        XASSERTM(this->size() == vector.size(), "size mismatch between mirror and vector");
 
        if(this->empty())
          return;
 
        LAFEM::Arch::Mirror::scatter_dv(
          buffer_offset, this->num_indices(), this->indices(), buffer.elements(), vector.elements(), alpha);
      }
 
      template<int block_size_>
      void gather(
        LAFEM::DenseVector<DataType, IndexType>& buffer,
        const LAFEM::DenseVectorBlocked<DataType, IndexType, block_size_>& vector,
        const Index buffer_offset = Index(0)) const
      {
        XASSERT(buffer_offset + Index(block_size_)*this->num_indices() <= buffer.size());
        XASSERTM(this->size() == vector.size(), "size mismatch between mirror and vector");
 
        if(this->empty())
          return;
 
        LAFEM::Arch::Mirror::gather_dvb(
          Index(block_size_), buffer_offset, this->num_indices(), this->indices(),
          buffer.elements(), vector.template elements<Perspective::pod>());
      }
 
      template<int block_size_>
      void scatter_axpy(
        LAFEM::DenseVectorBlocked<DataType, IndexType, block_size_>& vector,
        const LAFEM::DenseVector< DataType, IndexType>& buffer,
        const DataType alpha = DataType(1),
        const Index buffer_offset = Index(0)) const
      {
        XASSERT(buffer_offset + Index(block_size_)*this->num_indices() <= buffer.size());
        XASSERTM(this->size() == vector.size(), "size mismatch between mirror and vector");
 
        if(this->empty())
          return;
 
        LAFEM::Arch::Mirror::scatter_dvb(
          Index(block_size_), buffer_offset, this->num_indices(), this->indices(),
          buffer.elements(), vector.template elements<Perspective::pod>(), alpha);
      }
 
      void gather(
        LAFEM::DenseVector<DataType, IndexType>& buffer,
        const LAFEM::SparseVector<DataType, IndexType>& vector,
        const Index buffer_offset = Index(0)) const
      {
        XASSERT(buffer_offset + this->num_indices() <= buffer.size());
        XASSERTM(this->size() == vector.size(), "size mismatch between mirror and vector");
 
        if(this->empty())
          return;
 
        LAFEM::Arch::Mirror::gather_sv(
          buffer_offset, this->num_indices(), this->indices(), buffer.elements(),
          vector.used_elements(), vector.elements(), vector.indices());
      }
 
      void scatter_axpy(
        LAFEM::SparseVector<DataType, IndexType>& vector,
        const LAFEM::DenseVector<DataType, IndexType>& buffer,
        const DataType alpha = DataType(1),
        const Index buffer_offset = Index(0)) const
      {
        XASSERT(buffer_offset + this->num_indices() <= buffer.size());
        XASSERTM(this->size() == vector.size(), "size mismatch between mirror and vector");
 
        if(this->empty())
          return;
 
        LAFEM::Arch::Mirror::scatter_sv(
          buffer_offset, this->num_indices(), this->indices(), buffer.elements(),
          vector.used_elements(), vector.elements(), vector.indices(), alpha);
      }
 
      template<int block_size_>
      void gather(
        LAFEM::DenseVector<DataType, IndexType>& buffer,
        const LAFEM::SparseVectorBlocked<DataType, IndexType, block_size_>& vector,
        const Index buffer_offset = Index(0)) const
      {
        XASSERT(buffer_offset + Index(block_size_)*this->num_indices() <= buffer.size());
        XASSERTM(this->size() == vector.size(), "size mismatch between mirror and vector");
 
        if(this->empty())
          return;
 
        LAFEM::Arch::Mirror::gather_svb(
          Index(block_size_), buffer_offset, this->num_indices(), this->indices(),buffer.elements(),
          vector.used_elements(), vector.template elements<Perspective::pod>(), vector.indices());
      }
 
      template<int block_size_>
      void scatter_axpy(
        LAFEM::SparseVectorBlocked<DataType, IndexType, block_size_>& vector,
        const LAFEM::DenseVector<DataType, IndexType>& buffer,
        const DataType alpha = DataType(1),
        const Index buffer_offset = Index(0)) const
      {
        XASSERT(buffer_offset + Index(block_size_)*this->num_indices() <= buffer.size());
        XASSERTM(this->size() == vector.size(), "size mismatch between mirror and vector");
 
        if(this->empty())
          return;
 
        LAFEM::Arch::Mirror::scatter_svb(
          Index(block_size_), buffer_offset, this->num_indices(), this->indices(), buffer.elements(),
          vector.used_elements(), vector.template elements<Perspective::pod>(), vector.indices(), alpha);
      }
 
      template<Perspective perspective_, typename DT2_, typename IT2_>
      Index mask_scatter(const DenseVector<DT2_, IT2_>& vector, std::vector<int>& mask,
        const int value, const Index offset = Index(0)) const
      {
        XASSERT(Index(mask.size()) >= vector.template size<perspective_>());
        XASSERT(Index(mask.size()) >= this->size() + offset);
        const Index n = this->num_indices();
        const IT_* idx = this->indices();
 
        for(Index i(0); i < n; ++i)
          mask[offset + idx[i]] = value;
 
        return vector.template size<perspective_>();
      }
 
      template<Perspective perspective_, typename DT2_, typename IT2_, int block_size_>
      Index mask_scatter(const DenseVectorBlocked<DT2_, IT2_, block_size_>& vector, std::vector<int>& mask,
        const int value, const Index offset = Index(0)) const
      {
        XASSERT(Index(mask.size()) >= vector.template size<perspective_>());
        const Index n = this->num_indices();
        const IT_* idx = this->indices();
 
        if(perspective_ == LAFEM::Perspective::native)
        {
          XASSERT(Index(mask.size()) >= this->size() + offset);
          for(Index i(0); i < n; ++i)
            mask[offset + idx[i]] = value;
        }
        else // POD
        {
          XASSERT(Index(mask.size()) >= Index(block_size_)*this->size() + offset);
          for(Index i(0); i < n; ++i)
          {
            const Index ibs = idx[i] * Index(block_size_);
            for(int k(0); k < block_size_; ++k)
              mask[offset + ibs + Index(k)] = value;
          }
        }
        return vector.template size<perspective_>();
      }
 
 
      friend std::ostream & operator<< (std::ostream & lhs, const VectorMirror & b)
      {
        Index n = b.num_indices();
        const IT_* idx = b.indices();
 
        lhs << "[";
        for (Index i(0) ; i < n ; ++i)
        {
          lhs << "  " << idx[i];
        }
        lhs << "]";
 
        return lhs;
      }
 
      /* ******************************************************************* */
      /*  A D J A C T O R   I N T E R F A C E   I M P L E M E N T A T I O N  */
      /* ******************************************************************* */
    public:
      Index get_num_nodes_domain() const
      {
        return this->num_indices();
      }
 
      Index get_num_nodes_image() const
      {
        return this->size();
      }
 
      ImageIterator image_begin(Index domain_node) const
      {
        XASSERTM(domain_node < num_indices(), "Domain node index out of range");
        return &this->_indices.at(0)[domain_node];
      }
 
      ImageIterator image_end(Index domain_node) const
      {
        XASSERTM(domain_node < num_indices(), "Domain node index out of range");
        return &this->_indices.at(0)[domain_node+1];
      }
    }; // class VectorMirror<...>
  } // namespace LAFEM
} // namespace FEAT