VectorBackend.hpp

#pragma once
 
#include <boost/numeric/mtl/vector/dense_vector.hpp>
 
#include <dune/common/ftraits.hh>
 
#include <amdis/Output.hpp>
#include <amdis/common/FakeContainer.hpp>
#include <amdis/linearalgebra/VectorBase.hpp>
#include <amdis/typetree/MultiIndex.hpp>
 
namespace AMDiS
{
  class DefaultIndexDistribution;
 
  template <class T>
  class MTLVector
      : public VectorBase<MTLVector<T>>
  {
    using Self = MTLVector;
 
  public:
    using BaseVector = mtl::dense_vector<T>;
 
    using value_type = typename BaseVector::value_type;
 
    using size_type  = typename BaseVector::size_type;
 
  private:
    using block_type = value_type;
 
    using field_type = typename Dune::FieldTraits<value_type>::field_type;
 
  public:
    explicit MTLVector(DefaultIndexDistribution const&) {}
 
    BaseVector const& vector() const
    {
      return vector_;
    }
 
    BaseVector& vector()
    {
      return vector_;
    }
 
    std::size_t size() const
    {
      return mtl::vec::size(vector_);
    }
 
 
    template <class Basis>
    void init(Basis const& basis, bool clear)
    {
      vector_.change_dim(basis.dimension());
      if (clear)
        set_to_zero(vector_);
    }
 
    template <class MultiIndex>
    value_type& at(MultiIndex const& idx)
    {
      const size_type i = flatMultiIndex(idx);
      test_exit_dbg(i < size(),"Index {} out of range [0,{})", i, size());
      return vector_[i];
    }
 
    template <class MultiIndex>
    value_type const& at(MultiIndex const& idx) const
    {
      const size_type i = flatMultiIndex(idx);
      test_exit_dbg(i < size(), "Index {} out of range [0,{})", i, size());
      return vector_[i];
    }
 
 
    Self& axpy (const value_type& alpha, const Self& x)
    {
      vector_ += alpha * x.vector_;
      return *this;
    }
 
    Self& aypx (const value_type& alpha, const Self& x)
    {
      vector_ = alpha * vector_ + x.vector_;
      return *this;
    }
 
    void set (const value_type& alpha)
    {
      vector_ = alpha;
    }
 
    void scale (const value_type& alpha)
    {
      vector_ *= alpha;
    }
 
 
  private:
    BaseVector vector_;
  };
 
 
  namespace Recursive
  {
    template <class T>
    struct ForEach<MTLVector<T>> : ForEach<VectorBase<MTLVector<T>>> {};
 
    template <class T>
    struct Transform<MTLVector<T>> : Transform<VectorBase<MTLVector<T>>> {};
 
    template <class T>
    struct InnerProduct<MTLVector<T>> : InnerProduct<VectorBase<MTLVector<T>>> {};
 
  } // end namespace Recursive
} // end namespace AMDiS