product.hh

#ifndef __DUNE_ACFEM_TENSORS_OPERATIONS_PRODUCT_HH__
#define __DUNE_ACFEM_TENSORS_OPERATIONS_PRODUCT_HH__
 
#include "../../expressions/storage.hh"
#include "../../expressions/expressionoperations.hh" // FieldPromotion
#include "../../mpl/insertat.hh"
 
#include "../tensorbase.hh"
#include "../modules.hh"
#include "restrictiondetail.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    // forward
    template<class Pos1, class Pos2, class ProductDims>
    struct TensorProductOperation;
 
    namespace Tensor {
 
      template<class LeftTensor, class LeftIndices, class RightTensor, class RightIndices>
      class ProductTensor;
 
      template<class Left, class LPos, class Right, class RPos, class SFINAE = void>
      struct ProductOperationTraits
      {
        using Signature = SequenceCat<SequenceSlice<Left, LPos>,
                                      SequenceSliceComplement<Left, LPos>,
                                      SequenceSliceComplement<Right, RPos> >;
        using Operation = TensorProductOperation<LPos, RPos, SequenceSlice<Left, LPos> >;
        using Functor = OperationTraits<Operation>;
      };
 
      template<class Left, class LPos, class Right, class RPos>
      struct ProductOperationTraits<Left, LPos, Right, RPos,
                                    std::enable_if_t<(IsTensorOperand<Left>::value
                                                      && IsTensorOperand<Right>::value
        )> >
        : ProductOperationTraits<typename TensorTraits<Left>::Signature, LPos,
                                 typename TensorTraits<Right>::Signature, RPos>
      {};
 
      template<class Left, class LPos, class Right, class RPos>
      using ProductSignature = typename ProductOperationTraits<typename TensorTraits<Left>::TensorType, LPos,
                                                               typename TensorTraits<Right>::TensorType, RPos>::Signature;
 
      template<class Left, class LPos, class Right, class RPos>
      using ProductFunctor = typename ProductOperationTraits<typename TensorTraits<Left>::TensorType, LPos,
                                                             typename TensorTraits<Right>::TensorType, RPos>::Functor;
 
      template<class T>
      using TotalProductOperation = TensorProductOperation<
        MakeIndexSequence<TensorTraits<T>::rank>,
        MakeIndexSequence<TensorTraits<T>::rank>,
        typename TensorTraits<T>::Signature
        >;
 
      template<class LeftTensor, std::size_t... LeftIndices, class RightTensor, std::size_t... RightIndices>
      class ProductTensor<LeftTensor, Seq<LeftIndices...>, RightTensor, Seq<RightIndices...> >
        : public TensorBase<typename FieldPromotion<LeftTensor, RightTensor>::Type,
                            ProductSignature<LeftTensor, Seq<LeftIndices...>, RightTensor, Seq<RightIndices...> >,
                            ProductTensor<LeftTensor, Seq<LeftIndices...>, RightTensor, Seq<RightIndices...> > >
        , public Expressions::Storage<ProductFunctor<LeftTensor, Seq<LeftIndices...>, RightTensor, Seq<RightIndices...> >,
                                      LeftTensor, RightTensor>
      // Expressions model rvalues. We are constant if the underlying
      // expression claims to be so or if we store a copy and the
      // underlying expression is "independent" (i.e. determined at
      // runtime only by the contained data)
      {
        using ThisType = ProductTensor;
        using LeftType = LeftTensor;
        using RightType = RightTensor;
        using LeftTensorType = std::decay_t<LeftTensor>;
        using RightTensorType = std::decay_t<RightTensor>;
        using LeftSignature = typename LeftTensorType::Signature;
        using RightSignature = typename RightTensorType::Signature;
        using LeftDefectSignature = SubSequence<LeftSignature, LeftIndices...>;
        using RightDefectSignature = SubSequence<RightSignature, RightIndices...>;
        using LeftSignatureRest = SubSequenceComplement<LeftSignature, LeftIndices...>;
        using RightSignatureRest = SubSequenceComplement<RightSignature, RightIndices...>;
       public:
        using LeftIndexPositions = Seq<LeftIndices...>;
        using RightIndexPositions = Seq<RightIndices...>;
        using Signature = ProductSignature<LeftTensor, LeftIndexPositions, RightTensor, RightIndexPositions>;
        using FieldType = typename FieldPromotion<LeftType, RightType>::Type;
        using FunctorType = OperationTraits<TensorProductOperation<LeftIndexPositions, RightIndexPositions, LeftDefectSignature> >;
       private:
        // Mmmh. Should we really support this?
        using LeftSorted = SortSequence<LeftIndexPositions>;
        using RightSorted = SortSequence<RightIndexPositions>;
        using LeftPermutation = typename LeftSorted::Permutation;
        using RightPermutation = typename RightSorted::Permutation;
        using LeftSortedPos = typename LeftSorted::Result;
        using RightSortedPos = typename RightSorted::Result;
        static constexpr bool leftIsSorted_ = isSimple(LeftPermutation{});
        static constexpr bool rightIsSorted_ = isSimple(RightPermutation{});
        using BaseType = TensorBase<FieldType, Signature, ProductTensor>;
        using StorageType = Expressions::Storage<FunctorType, LeftTensor, RightTensor>;
       public:
        using StorageType::operation;
        using StorageType::operand;
        using BaseType::rank;
        static constexpr std::size_t leftRank_ = LeftSignatureRest::size();
        static constexpr std::size_t rightRank_ = RightSignatureRest::size();
        static constexpr std::size_t frontRank_ = rank - leftRank_ - rightRank_;
        using FrontArgs = MakeIndexSequence<frontRank_>;
        using LeftArgs = MakeIndexSequence<leftRank_, frontRank_>;
        using RightArgs = MakeIndexSequence<rightRank_, frontRank_ + leftRank_>;
       private:
 
        static_assert(sizeof...(LeftIndices) == sizeof...(RightIndices),
                      "Number of contraction indices must coincide.");
        static_assert(LeftTensorType::rank >= sizeof...(LeftIndices),
                      "Left: Number of index-positions must fit into tensor rank.");
        static_assert(RightTensorType::rank >= sizeof...(RightIndices),
                      "Right: Number of index-positions must fit into tensor rank.");
        static_assert(std::is_same<LeftDefectSignature, RightDefectSignature>::value,
                      "Left- and right defect-dimensions must coincide.");
 
       public:
        using DefectSignature = LeftDefectSignature;
        static constexpr std::size_t defectRank_ = DefectSignature::size();
 
        template<class LeftArg, class RightArg,
                 std::enable_if_t<std::is_constructible<LeftType, LeftArg>::value && std::is_constructible<RightType, RightArg>::value, int> = 0>
        ProductTensor(LeftArg&& left, RightArg&& right)
          : StorageType(std::forward<LeftArg>(left), std::forward<RightArg>(right))
        {}
 
        template<
          class... Dummy,
          std::enable_if_t<(sizeof...(Dummy) == 0
                            && IsTypedValue<LeftType>::value
                            && IsTypedValue<RightType>::value), int> = 0>
        ProductTensor(Dummy&&...)
          : StorageType(LeftType{}, RightType{})
        {}
 
        template<class... Dims,
                 std::enable_if_t<(sizeof...(Dims) == rank
                                   &&
                                   IsIntegralPack<Dims...>::value)
                                  , int> = 0>
        auto operator()(Dims... indices) const
        {
          auto frontIndices = forwardSubTuple(std::forward_as_tuple(indices...), FrontArgs{});
          auto leftIndices = forwardSubTuple(std::forward_as_tuple(indices...), LeftArgs{});
          auto rightIndices = forwardSubTuple(std::forward_as_tuple(indices...), RightArgs{});
 
          return
            tensorValue(operand(0_c), insertAt(leftIndices, permute(frontIndices, LeftPermutation{}), LeftSortedPos{}))
            *
            tensorValue(operand(1_c), insertAt(rightIndices, permute(frontIndices, RightPermutation{}), RightSortedPos{}));
        }
 
        template<std::size_t... Indices,
                 std::enable_if_t<(sizeof...(Indices) == rank
                                   && ThisType::template isZero<Indices...>()
          ), int> = 0>
        auto constexpr operator()(Seq<Indices...>) const
        {
          return IntFraction<0>{};
        }
 
        template<std::size_t... Indices,
                 std::enable_if_t<(sizeof...(Indices) == rank
                                   && !ThisType::template isZero<Indices...>()
          ), int> = 0>
        auto constexpr operator()(Seq<Indices...>) const
        {
          using IndexSeq = Seq<Indices...>;
          using ArgFrontIndices = SequenceSlice<IndexSeq, FrontArgs>;
          using ArgLeftIndices = SequenceSlice<IndexSeq, LeftArgs>;
          using ArgRightIndices = SequenceSlice<IndexSeq, RightArgs>;
          using LeftArg = InsertAt<ArgLeftIndices, PermuteSequence<ArgFrontIndices, LeftPermutation>, LeftSortedPos>;
          using RightArg = InsertAt<ArgRightIndices, PermuteSequence<ArgFrontIndices, RightPermutation>, RightSortedPos>;
 
          return operand(0_c)(LeftArg{}) * operand(1_c)(RightArg{});
        }
 
       private:
        template<std::size_t... Indices, class Pos = MakeIndexSequence<sizeof...(Indices)> >
        static bool constexpr isZeroWorker(Seq<Indices...> = Seq<Indices...>{}, Pos = Pos{})
        {
          // This ain't pretty ;)
 
          // Truncate Pos to size of Indices and sort it, otherwise
          // the JoinedDefects stuff will not work properly
          using TruncatedPos = HeadPart<sizeof...(Indices), Pos>;
          using RealPos = typename SortSequence<TruncatedPos>::Result;
          using Permutation = typename SortSequence<TruncatedPos>::Permutation;
          using RealIndices = PermuteSequence<Seq<Indices...>, Permutation>;
 
          // Extract leading product indices
          using FrontPosInd =
            TransformSequence<RealPos, Seq<>, IndexFunctor, AcceptInputInRangeFunctor<std::size_t, 0, frontRank_> >;
          using FrontInd = SequenceSlice<RealIndices, FrontPosInd>;
          using FrontPos = SequenceSlice<RealPos, FrontPosInd>;
 
          using LeftFrontPos = SequenceSlice<LeftIndexPositions, FrontPos>;
          using RightFrontPos = SequenceSlice<RightIndexPositions, FrontPos>;
 
          // Extract the indices referring to the left Tensor
          using LeftPosInd =
            TransformSequence<RealPos, Seq<>, IndexFunctor, AcceptInputInRangeFunctor<std::size_t, frontRank_, frontRank_ + leftRank_> >;
          using LeftPos = TransformSequence<SequenceSlice<RealPos, LeftPosInd>, Seq<>, OffsetFunctor<-(std::ptrdiff_t)frontRank_> >;
          using LeftInd = SequenceSlice<RealIndices, LeftPosInd>;
 
          // The total set of left index positions, contraction and taken from Indices...
          using LeftLookup = SequenceSliceComplement<MakeIndexSequence<frontRank_+leftRank_>, LeftSortedPos>;
          using LeftTotalPos = SequenceCat<LeftFrontPos, TransformedSequence<MapSequenceFunctor<LeftLookup>, LeftPos> >;
          using LeftTotalInd = SequenceCat<FrontInd, LeftInd>;
 
          // Extract the indices referring to the right Tensor
          using RightPosInd =
            TransformSequence<RealPos, Seq<>, IndexFunctor, AcceptInputInRangeFunctor<std::size_t, frontRank_ + leftRank_, frontRank_ + leftRank_ + rightRank_> >;
          using RightPos = TransformSequence<SequenceSlice<RealPos, RightPosInd>, Seq<>, OffsetFunctor<-(std::ptrdiff_t)(frontRank_+leftRank_)> >;
          using RightInd = SequenceSlice<RealIndices, RightPosInd>;
 
          // The total set of right index positions, contraction and taken from Indices...
          using RightLookup = SequenceSliceComplement<MakeIndexSequence<frontRank_+rightRank_>, RightSortedPos>;
          using RightTotalPos = SequenceCat<RightFrontPos, TransformedSequence<MapSequenceFunctor<RightLookup>, RightPos> >;
          using RightTotalInd = SequenceCat<FrontInd, RightInd>;
 
          return
            LeftTensorType::isZero(LeftTotalInd{}, LeftTotalPos{})
            ||
            RightTensorType::isZero(RightTotalInd{}, RightTotalPos{});
        }
 
       public:
        template<std::size_t... Indices, class Pos = MakeIndexSequence<sizeof...(Indices)> >
        static bool constexpr isZero(Seq<Indices...> = Seq<Indices...>{}, Pos = Pos{})
        {
//          if constexpr (ExpressionTraits<LeftTensorType>::isZero || ExpressionTraits<RightTensorType>::isZero) {
//            return true;
//          } else {
            return isZeroWorker(Seq<Indices...>{}, Pos{});
//          }
        }
 
        std::string name() const
        {
          using TL = LeftTensor;
          using TR = RightTensor;
          std::string pfxL = std::is_reference<TL>::value ? (RefersConst<TL>::value ? "cref" : "ref") : "";
          std::string pfxR = std::is_reference<TR>::value ? (RefersConst<TR>::value ? "cref" : "ref") : "";
 
          return operationName(operation(), pfxL+operand(0_c).name(), pfxR+operand(1_c).name());
        }
 
      }; // ProductTensor class
 
      template<class Seq1, class Seq2, class T1, class T2,
               std::enable_if_t<AreProperTensors<T1, T2>::value, int> = 0>
      constexpr decltype(auto) multiply(T1&& t1, T2&& t2)
      {
        return Expressions::finalize(ProductFunctor<T1, Seq1, T2, Seq2>{}, std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      template<class T1, class T2,
               std::enable_if_t<AreProperTensors<T1, T2>::value, int> = 0>
      constexpr decltype(auto) multiply(T1&& t1, T2&& t2)
      {
        constexpr std::size_t numDims = CommonHead<typename TensorTraits<T1>::Signature,
                                                   typename TensorTraits<T2>::Signature>::value;
        using ContractPos = MakeIndexSequence<numDims>;
        return multiply<ContractPos, ContractPos>(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      template<class Seq0, class Seq1, class Dims, class T0, class T1>
      constexpr auto operate(Expressions::DontOptimize, OperationTraits<TensorProductOperation<Seq0, Seq1, Dims> >, T0&& t0, T1&& t1)
      {
        DUNE_ACFEM_RECORD_OPTIMIZATION;
 
        return ProductTensor<T0, Seq0, T1, Seq1>(std::forward<T0>(t0), std::forward<T1>(t1));
      }
 
    } // NS Tensor
 
    // Point here is that the operations defined in
    // ./expressionoperations.hh need the name "multiply".
    using Tensor::multiply;
 
  } // NS ACFem
 
  template<class LeftTensor, class LeftIndices, class RightTensor, class RightIndices>
  struct FieldTraits<ACFem::Tensor::ProductTensor<LeftTensor, LeftIndices, RightTensor, RightIndices> >
    : FieldTraits<typename ACFem::FieldPromotion<LeftTensor, RightTensor>::Type>
  {};
 
 
 
} // NS Dune
 
#endif // __DUNE_ACFEM_TENSORS_OPERATIONS_PRODUCT_HH__