transpose.hh

#ifndef __DUNE_ACFEM_TENSORS_OPERATIONS_TRANSPOSE_HH__
#define __DUNE_ACFEM_TENSORS_OPERATIONS_TRANSPOSE_HH__
 
#include "../../expressions/storage.hh"
#include "../../mpl/permutation.hh"
#include "../../mpl/sort.hh"
#include "../tensorbase.hh"
#include "../modules.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    // forward
    template<class Perm>
    struct TransposeOperation;
 
    namespace Tensor {
 
      template<class Tensor, class Perm>
      class Transposition;
 
      template<class Tensor, std::size_t... Perm>
      class Transposition<Tensor, Seq<Perm...> >
        : public TensorBase<typename TensorTraits<Tensor>::FieldType,
                            Seq<Get<Perm, typename TensorTraits<Tensor>::Signature>::value...>,
                            Transposition<Tensor, Seq<Perm...> > >
        , public Expressions::Storage<OperationTraits<TransposeOperation<Seq<Perm...> > >, Tensor>
      {
        static_assert(sizeof...(Perm) == TensorTraits<Tensor>::rank,
                      "A transposition must be provided for all tensor dimensions.");
 
        using ArgType = Tensor;
        using HostType = std::decay_t<Tensor>;
        using BaseType = TensorBase<typename HostType::FieldType,
                                    Seq<Get<Perm, typename TensorTraits<Tensor>::Signature>::value...>,
                                    Transposition<Tensor, Seq<Perm...> > >;
        using StorageType = Expressions::Storage<OperationTraits<TransposeOperation<Seq<Perm...> > >, Tensor>;
      public:
        using Permutation = Seq<Perm...>;
        using StorageType::operation;
        using StorageType::operand;
        using BaseType::rank;
        using typename BaseType::Signature;
        using typename BaseType::FieldType;
 
       public:
        template<class Arg, std::enable_if_t<std::is_constructible<ArgType, Arg>::value, int> = 0>
        Transposition(Arg&& host)
          : StorageType(std::forward<Arg>(host))
        {}
 
        template<
          class... Dummy,
          std::enable_if_t<(sizeof...(Dummy) == 0
                            && IsTypedValue<ArgType>::value), int> = 0>
        Transposition(Dummy&&...)
          : StorageType(ArgType{})
        {}
 
        template<class... Dims,
                 std::enable_if_t<(sizeof...(Dims) == rank
                                   &&
                                   IsIntegralPack<Dims...>::value)
                                  , int> = 0>
        decltype(auto) operator()(Dims... indices)
        {
          return operand(0_c)(std::get<Perm>(std::forward_as_tuple(indices...))...);
        }
 
        template<class... Dims,
                 std::enable_if_t<(sizeof...(Dims) == rank
                                   &&
                                   IsIntegralPack<Dims...>::value)
                                  , int> = 0>
        decltype(auto) operator()(Dims... indices) const
        {
          return operand(0_c)(std::get<Perm>(std::forward_as_tuple(indices...))...);
        }
 
        template<std::size_t... Indices, std::enable_if_t<sizeof...(Indices) == rank, int> = 0>
        decltype(auto) constexpr operator()(Seq<Indices...>) const
        {
          return operand(0_c)(Seq<Get<Perm, Seq<Indices...> >::value...>{});
        }
 
        template<std::size_t... Indices, class Pos = MakeIndexSequence<sizeof...(Indices)> >
        static bool constexpr isZero(Seq<Indices...> = Seq<Indices...>{}, Pos = Pos{})
        {
          // Truncate Pos to size of Indices...
          using RealPos = HeadPart<sizeof...(Indices), Pos>;
 
          // Map to the transposed positions.
          using TransposedPos = PermuteSequenceValues<RealPos, Seq<Perm...> >;
 
          using SortedPos = typename SortSequence<TransposedPos>::Result;
          using Permutation = typename SortSequence<TransposedPos>::Permutation;
          using SortedIndices = PermuteSequence<Seq<Indices...>, Permutation>;
 
          return HostType::isZero(SortedIndices{}, SortedPos{});
        }
 
        std::string name() const
        {
          using T = Tensor;
          std::string pfx = std::is_reference<T>::value ? (RefersConst<T>::value ? "cref" : "ref") : "";
          return operationName(operation(), pfx+operand(0_c).name());
        }
 
      };
 
      template<class T>
      struct IsTransposition
        : FalseType
      {};
 
      template<class T, std::size_t... Perm>
      struct IsTransposition<Transposition<T, Seq<Perm...> > >
        : TrueType
      {};
 
      template<std::size_t... Perm, class T,
               std::enable_if_t<(IsProperTensor<T>::value
                                 && isPermutation(ResizedPermutation<Seq<Perm...>, TensorTraits<T>::rank>{})
                                 && sizeof...(Perm) <= TensorTraits<T>::rank
        ), int> = 0>
      constexpr decltype(auto) transpose(T&& t, Seq<Perm...> = Seq<Perm...>{})
      {
        using Operation = TransposeOperation<ResizedPermutation<Seq<Perm...>, TensorTraits<T>::rank> >;
 
        return Expressions::finalize<Operation>(std::forward<T>(t));
      }
 
      template<class Perm, class T, std::enable_if_t<Perm::size() == TensorTraits<T>::rank, int> = 0>
      constexpr auto operate(Expressions::DontOptimize, OperationTraits<TransposeOperation<Perm> >, T&& t)
      {
        DUNE_ACFEM_RECORD_OPTIMIZATION;
 
        return Transposition<T, Perm>(std::forward<T>(t));
      }
 
    } // NS Tensor
 
    //using Tensor::transpose;
 
  } // NS ACFem
 
  template<class T, class P>
  struct FieldTraits<ACFem::Tensor::Transposition<T, P> >
    : FieldTraits<std::decay_t<T> >
  {};
 
 
 
} // NS Dune
 
#endif // __DUNE_ACFEM_TENSORS_OPERATIONS_TRANSPOSE_HH__