einsum.hh

#ifndef __DUNE_ACFEM_TENSORS_OPERATIONS_EINSUM_HH__
#define __DUNE_ACFEM_TENSORS_OPERATIONS_EINSUM_HH__
 
#include "../../common/types.hh"
//#include "../../common/ostream.hh"
#include "../../mpl/foreach.hh"
#include "../../mpl/insertat.hh"
#include "../../mpl/toarray.hh"
#include "../../expressions/storage.hh"
#include "../../expressions/expressionoperations.hh"
#include "../../expressions/constantoperations.hh"
#include "../tensorbase.hh"
#include "../modules/eye.hh"
#include "../optimization/policy.hh"
#include "restrictiondetail.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    // forward
    template<class Pos1, class Pos2, class ContractDims>
    struct EinsumOperation;
 
    namespace Tensor {
 
      template<class LeftTensor, class LeftIndices, class RightTensor, class RightIndices, bool IdenticalOperands = Expressions::AreRuntimeEqual<LeftTensor, RightTensor>::value>
      class EinsteinSummation;
 
      template<class LDims, class LPos, class RDims, class RPos, class SFINAE = void>
      struct EinsumOperationTraits
        : EinsumOperationTraits<typename TensorTraits<LDims>::Signature, LPos,
                                typename TensorTraits<RDims>::Signature, RPos>
      {};
 
      template<std::size_t... LDims, std::size_t... LPos, std::size_t... RDims, std::size_t... RPos>
      struct EinsumOperationTraits<Seq<LDims...>, Seq<LPos...>, Seq<RDims...>, Seq<RPos...> >
      {
        using Signature = SequenceCat<SubSequenceComplement<Seq<LDims...>, LPos...>,
                                      SubSequenceComplement<Seq<RDims...>, RPos...> >;
        static constexpr std::size_t rank_ = Signature::size();
        using Operation = EinsumOperation<Seq<LPos...>, Seq<RPos...>, SubSequence<Seq<LDims...>, LPos...> >;
        using Functor = OperationTraits<Operation>;
      };
 
      template<class LDims, class LPos, class RDims, class RPos>
      using EinsumSignature = typename EinsumOperationTraits<LDims, LPos, RDims, RPos>::Signature;
 
      template<class LDims, class LPos, class RDims, class RPos>
      using EinsumFunctor = typename EinsumOperationTraits<LDims, LPos, RDims, RPos>::Functor;
 
      using ScalarEinsumOperation = EinsumOperation<Seq<>, Seq<>, Seq<> >;
      using ScalarEinsumFunctor = OperationTraits<ScalarEinsumOperation>;
 
      template<class T0, class T1, bool IdenticalOperands>
      using ScalarEinsumExpression = EinsteinSummation<T0, Seq<>, T1, Seq<>, IdenticalOperands>;
 
      template<class LDims, class LPos, class RDims, class RPos>
      inline constexpr std::size_t EinsumRank = EinsumOperationTraits<LDims, LPos, RDims, RPos>::rank_;
 
      template<class Dims, class Pos>
      using EinsumDimensions = SequenceSlice<Dims, Pos>;
 
      template<class LeftTensor, class LeftIndices, class RightTensor, class RightIndices, bool IdenticalOperands>
      class EinsteinSummation
        : public TensorBase<FloatingPointClosure<typename FieldPromotion<LeftTensor, RightTensor>::Type>,
                            EinsumSignature<LeftTensor, LeftIndices,
                                            RightTensor, RightIndices >,
                            EinsteinSummation<LeftTensor, LeftIndices, RightTensor, RightIndices > >
        , public Expressions::Storage<EinsumFunctor<LeftTensor, LeftIndices, RightTensor, 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 = EinsteinSummation;
        using LeftType = LeftTensor;
        using RightType = RightTensor;
        using LeftTensorType = std::decay_t<LeftTensor>;
        using RightTensorType = std::decay_t<RightTensor>;
       public:
        using LeftIndexPositions = LeftIndices;
        using RightIndexPositions = RightIndices;
       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 LeftSignature = typename LeftTensorType::Signature;
        using RightSignature = typename RightTensorType::Signature;
        using LeftDefectSignature = SequenceSlice<LeftSignature, LeftIndices>;
        using RightDefectSignature = SequenceSlice<RightSignature, RightIndices>;
       public:
        using LeftSignatureRest = SequenceSliceComplement<LeftSignature, LeftIndices>;
        using RightSignatureRest = SequenceSliceComplement<RightSignature, RightIndices>;
        using Signature = EinsumSignature<LeftTensor, LeftIndexPositions, RightTensor, RightIndexPositions>;
        using FieldType = FloatingPointClosure<typename FieldPromotion<LeftTensor, RightTensor>::Type>;
        using FunctorType = OperationTraits<EinsumOperation<LeftIndexPositions, RightIndexPositions, LeftDefectSignature> >;
       private:
        using BaseType = TensorBase<FieldType, Signature, EinsteinSummation>;
        using StorageType = Expressions::Storage<FunctorType, LeftTensor, RightTensor>;
       public:
        static const std::size_t leftRank_ = LeftSignatureRest::size();
        static const std::size_t rightRank_ = RightSignatureRest::size();
        using LeftArgs = MakeIndexSequence<leftRank_>;
        using RightArgs = MakeIndexSequence<rightRank_, leftRank_>; // offset
        using DefectSignature = LeftDefectSignature;
        static constexpr std::size_t defectRank_ = DefectSignature::size();
       private:
        static_assert(LeftIndices::size() == RightIndices::size(),
                      "Number of contraction indices must coincide.");
        static_assert(LeftTensorType::rank >= LeftIndices::size(),
                      "Left: Number of index-positions must fit into tensor rank.");
        static_assert(RightTensorType::rank >= RightIndices::size(),
                      "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.");
        static_assert(LeftIndexPositions{} < LeftTensorType::rank,
                      "Left index positions out of range.");
        static_assert(RightIndexPositions{} < RightTensorType::rank,
                      "Right index positions out of range.");
 
       public:
        using StorageType::operation;
        using StorageType::operand;
        using BaseType::rank;
 
        template<class LeftArg, class RightArg,
                 std::enable_if_t<std::is_constructible<LeftType, LeftArg>::value && std::is_constructible<RightType, RightArg>::value, int> = 0>
        EinsteinSummation(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>
        EinsteinSummation(Dummy&&...)
          : StorageType(LeftType{}, RightType{})
        {}
 
#if 0
        static_assert(!(IsTypedValue<LeftType>::value && std::is_reference<LeftType>::value)
                      &&
                      !(IsTypedValue<RightType>::value && std::is_reference<RightType>::value),
                      "Typed values should not be stored as references.");
#endif
 
#if DUNE_ACFEM_TENSOR_WORKAROUND_GCC(7)
        class GCCBugCompensator1
        {
         public:
          template<class LeftIndexArg, class RightIndexArg>
          GCCBugCompensator1(FieldType& accu,
                             const LeftTensorType& left, const RightTensorType& right,
                             LeftIndexArg&& leftIndices, RightIndexArg&& rightIndices)
            : accu_(accu), left_(left), right_(right), leftIndices_(toArray<std::size_t>(leftIndices)), rightIndices_(toArray<std::size_t>(rightIndices))
          {}
 
          template<class I>
          void operator()(I)
          {
            using DefectIndices = MultiIndex<I::value, DefectSignature>;
            if (!LeftTensorType::isZero(DefectIndices{}, LeftIndexPositions{}) && !RightTensorType::isZero(DefectIndices{}, RightIndexPositions{})) {
              if constexpr (IdenticalOperands) {
                if constexpr (rank == 0 && std::is_same<LeftIndexPositions, RightIndexPositions>::value) {
                  auto val = tensorValue(left_, insertAt(leftIndices_, PermuteSequence<DefectIndices, LeftPermutation>{}, LeftSortedPos{}));
                  accu_ += val * val;
                } else {
                  accu_ +=
                    tensorValue(left_, insertAt(leftIndices_, PermuteSequence<DefectIndices, LeftPermutation>{}, LeftSortedPos{}))
                    *
                    tensorValue(left_, insertAt(rightIndices_, PermuteSequence<DefectIndices, RightPermutation>{}, RightSortedPos{}));
                }
              } else {
                accu_ +=
                  tensorValue(left_, insertAt(leftIndices_, PermuteSequence<DefectIndices, LeftPermutation>{}, LeftSortedPos{}))
                  *
                  tensorValue(right_, insertAt(rightIndices_, PermuteSequence<DefectIndices, RightPermutation>{}, RightSortedPos{}));
              }
            }
          }
         private:
          FieldType& accu_;
          const LeftTensorType& left_;
          const RightTensorType& right_;
          const std::array<std::size_t, leftRank_> leftIndices_;
          const std::array<std::size_t, rightRank_> rightIndices_;
        };
#endif
 
        template<class... Dims,
                 std::enable_if_t<(sizeof...(Dims) == rank
                                   &&
                                   IsIntegralPack<Dims...>::value)
                                  , int> = 0>
        auto operator()(Dims... indices) const
        {
          auto leftIndices = forwardSubTuple(std::forward_as_tuple(indices...), LeftArgs{});
          auto rightIndices = forwardSubTuple(std::forward_as_tuple(indices...), RightArgs{});
          (void)rightIndices;
 
          auto accu = FieldType(0);
#if DUNE_ACFEM_TENSOR_WORKAROUND_GCC(7)
          forLoop<multiDim(DefectSignature{})>(GCCBugCompensator1(accu, operand(0_c), operand(1_c), leftIndices, rightIndices));
#else
          forLoop<multiDim(DefectSignature{})>([&](auto i) {
              using I = decltype(i);
              using DefectIndices = MultiIndex<I::value, DefectSignature>;
              if (!LeftTensorType::isZero(DefectIndices{}, LeftIndexPositions{}) && !RightTensorType::isZero(DefectIndices{}, RightIndexPositions{})) {
                if constexpr (IdenticalOperands) {
                  if constexpr (rank == 0 && std::is_same<LeftIndexPositions, RightIndexPositions>::value) {
                    auto val = tensorValue(operand(0_c), insertAt(leftIndices, PermuteSequence<DefectIndices, LeftPermutation>{}, LeftSortedPos{}));
                    accu += val * val;
                  } else {
                    accu +=
                      tensorValue(operand(0_c), insertAt(leftIndices, PermuteSequence<DefectIndices, LeftPermutation>{}, LeftSortedPos{}))
                      *
                      tensorValue(operand(0_c), insertAt(rightIndices, PermuteSequence<DefectIndices, RightPermutation>{}, RightSortedPos{}));
                  }
                } else {
                  accu +=
                    tensorValue(operand(0_c), insertAt(leftIndices, PermuteSequence<DefectIndices, LeftPermutation>{}, LeftSortedPos{}))
                    *
                    tensorValue(operand(1_c), insertAt(rightIndices, PermuteSequence<DefectIndices, RightPermutation>{}, RightSortedPos{}));
                }
              }
            });
#endif
          return accu;
        }
 
#if DUNE_ACFEM_TENSOR_WORKAROUND_GCC(7)
        template<class LeftIndexArg, class RightIndexArg>
        class GCCBugCompensator2
        {
         public:
          GCCBugCompensator2(FieldType& accu, const LeftTensorType& left, const RightTensorType& right)
            : accu_(accu), left_(left), right_(right)
          {}
 
          template<class I>
          void operator()(I)
          {
            using DefectIndices = MultiIndex<I::value, DefectSignature>;
            using LeftArg = InsertAt<LeftIndexArg, PermuteSequence<DefectIndices, LeftPermutation>, LeftSortedPos>;
            using RightArg = InsertAt<RightIndexArg, PermuteSequence<DefectIndices, RightPermutation>, RightSortedPos>;
            if constexpr (IdenticalOperands) {
              if constexpr (rank == 0 && std::is_same<LeftIndexPositions, RightIndexPositions>::value) {
                if (!LeftTensorType::isZero(LeftArg{})) {
                  auto val = left_(LeftArg{});
                  accu_ += val * val;
                }
              } else {
                if (!LeftTensorType::isZero(LeftArg{}) && !LeftTensorType::isZero(RightArg{})) {
                  accu_ += left_(LeftArg{}) * left_(RightArg{});
                }
              }
            } else {
              if (!LeftTensorType::isZero(LeftArg{}) && !RightTensorType::isZero(RightArg{})) {
                accu_ += left_(LeftArg{}) * right_(RightArg{});
              }
            }
          }
         private:
          FieldType& accu_;
          const LeftTensorType& left_;
          const RightTensorType& right_;
        };
 
        template<class LeftIndexArg, class RightIndexArg>
        class GCCBugCompensator3
        {
         public:
          GCCBugCompensator3(const LeftTensorType& left, const RightTensorType& right)
            : left_(left), right_(right)
          {}
 
          template<class I>
          constexpr auto operator()(I) const
          {
            using DefectIndices = MultiIndex<I::value, DefectSignature>;
            using LeftArg = InsertAt<LeftIndexArg, PermuteSequence<DefectIndices, LeftPermutation>, LeftSortedPos>;
            using RightArg = InsertAt<RightIndexArg, PermuteSequence<DefectIndices, RightPermutation>, RightSortedPos>;
            if constexpr (IdenticalOperands) {
              if constexpr (rank == 0 && std::is_same<LeftIndexPositions, RightIndexPositions>::value) {
                auto val = IntFraction<(int)!LeftTensorType::isZero(LeftArg{})>{} * left_(LeftArg{});
                return val * val;
              } else {
                return
                  (IntFraction<(int)!LeftTensorType::isZero(LeftArg{})>{} * left_(LeftArg{}))
                  *
                  (IntFraction<(int)!RightTensorType::isZero(RightArg{})>{} * left_(RightArg{}));
              }
            } else {
              return
                (IntFraction<(int)!LeftTensorType::isZero(LeftArg{})>{} * left_(LeftArg{}))
                *
                (IntFraction<(int)!RightTensorType::isZero(RightArg{})>{} * right_(RightArg{}));
            }
          }
         private:
          const LeftTensorType& left_;
          const RightTensorType& right_;
        };
#endif
 
        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, int> = 0>
        auto constexpr operator()(Seq<Indices...>) const
        {
          static_assert(sizeof...(Indices) == rank, "This should not happen");
        }
 
        template<std::size_t... Indices,
                 std::enable_if_t<(sizeof...(Indices) == rank
                                   && !ThisType::template isZero<Indices...>()
                                   && multiDim(DefectSignature{}) >= 0
                                   && multiDim(DefectSignature{}) <= Policy::TemplateForEachLimit::value
          ), int> = 0>
        auto constexpr operator()(Seq<Indices...>) const
        {
          using IndexArg = Seq<Indices...>;
          using LeftIndexArg = SequenceSlice<IndexArg, LeftArgs>;
          using RightIndexArg = SequenceSlice<IndexArg, RightArgs>;
          static_assert(multiDim(DefectSignature{}) > 0, "BUG");
#if 1
          return
            addLoop<multiDim(DefectSignature{})>(
#if DUNE_ACFEM_TENSOR_WORKAROUND_GCC(7)
                              GCCBugCompensator3<LeftIndexArg, RightIndexArg>(operand(0_c), operand(1_c)),
#else
                              [&](auto i) {
                                using DefectIndices = MultiIndex<i(), DefectSignature>;
                                using LeftArg = InsertAt<LeftIndexArg, PermuteSequence<DefectIndices, LeftPermutation>, LeftSortedPos>;
                                using RightArg = InsertAt<RightIndexArg, PermuteSequence<DefectIndices, RightPermutation>, RightSortedPos>;
                                if constexpr (IdenticalOperands) {
                                  if constexpr (rank == 0 && std::is_same<LeftIndexPositions, RightIndexPositions>::value) {
                                    auto val = intFraction<!LeftTensorType::isZero(LeftArg{})>() * operand(0_c)(LeftArg{});
                                    return val * val;
                                  } else {
                                    return
                                      (intFraction<!LeftTensorType::isZero(LeftArg{})>() * operand(0_c)(LeftArg{}))
                                      *
                                      (intFraction<!RightTensorType::isZero(RightArg{})>() * operand(0_c)(RightArg{}));
                                  }
                                } else {
                                  return
                                    (intFraction<!LeftTensorType::isZero(LeftArg{})>() * operand(0_c)(LeftArg{}))
                                    *
                                    (intFraction<!RightTensorType::isZero(RightArg{})>() * operand(1_c)(RightArg{}));
                                }
                              },
#endif
                              intFraction<0>());
#else // if 0
          auto accu = FieldType(0);
#if DUNE_ACFEM_TENSOR_WORKAROUND_GCC(7)
          forLoop<multiDim(DefectSignature{})>(GCCBugCompensator2<LeftIndexArg, RightIndexArg>(accu, operand(0_c), operand(1_c)));
#else
          forLoop<multiDim(DefectSignature{})>([&](auto i) {
              using DefectIndices = MultiIndex<i, DefectSignature>;
              using LeftArg = InsertAt<LeftIndexArg, PermuteSequence<DefectIndices, LeftPermutation>, LeftSortedPos>;
              using RightArg = InsertAt<RightIndexArg, PermuteSequence<DefectIndices, RightPermutation>, RightSortedPos>;
 
              assert(!LeftTensorType::isZero(LeftArg{}) || operand(0_c)(LeftArg{}) == 0);
              assert(!RightTensorType::isZero(RightArg{}) || operand(1_c)(RightArg{}) == 0);
 
              if (!LeftTensorType::isZero(LeftArg{}) && !RightTensorType::isZero(RightArg{})) {
                accu += operand(0_c)(LeftArg{}) * operand(1_c)(RightArg{});
              }
            });
#endif
          return accu;
#endif
        }
 
        template<std::size_t... Indices,
                 std::enable_if_t<(sizeof...(Indices) == rank
                                   && !ThisType::template isZero<Indices...>()
                                   && multiDim(DefectSignature{}) > Policy::TemplateForEachLimit::value
          ), int> = 0>
        auto constexpr operator()(Seq<Indices...>) const
        {
#if 0
          // Just forward to the run-time method. Perhaps the compiler
          // still can optimize something ...
          return (*this)(Indices...);
#else
          // Still a compile-time loop, but with pack expansion
          // instead of recursion. Still quite costly ...  This will
          // never return an integral constant, as the pack-expansion
          // loop cannot return values.
          using IndexArg = Seq<Indices...>;
          using LeftIndexArg = SequenceSlice<IndexArg, LeftArgs>;
          using RightIndexArg = SequenceSlice<IndexArg, RightArgs>;
          auto accu = FieldType(0);
#if DUNE_ACFEM_TENSOR_WORKAROUND_GCC(7)
          forLoop<multiDim(DefectSignature{})>(GCCBugCompensator2<LeftIndexArg, RightIndexArg>(accu, operand(0_c), operand(1_c)));
# else
          forLoop<multiDim(DefectSignature{})>([&](auto i) {
              using I = decltype(i);
              using DefectIndices = MultiIndex<I::value, DefectSignature>;
              using LeftArg = InsertAt<LeftIndexArg, PermuteSequence<DefectIndices, LeftPermutation>, LeftSortedPos>;
              using RightArg = InsertAt<RightIndexArg, PermuteSequence<DefectIndices, RightPermutation>, RightSortedPos>;
 
              assert(!LeftTensorType::isZero(LeftArg{}) || operand(0_c)(LeftArg{}) == 0);
              assert(!RightTensorType::isZero(RightArg{}) || operand(1_c)(RightArg{}) == 0);
 
              if (!LeftTensorType::isZero(LeftArg{}) && !RightTensorType::isZero(RightArg{})) {
                accu += operand(0_c)(LeftArg{}) * operand(1_c)(RightArg{});
              }
          });
# endif
          return accu;
#endif
        }
 
       private:
 
        template<class LeftInd, class LeftInj, class LeftPos,
                 class RightInd, class RightInj, class RightPos,
                 std::size_t... N>
        static bool constexpr isZeroExpander(LeftInd, LeftInj, LeftPos,
                                             RightInd, RightInj, RightPos,
                                             Seq<N...>)
        {
          return (... && (LeftTensorType::isZero(InsertAt<PermuteSequence<MultiIndex<N, DefectSignature>, LeftPermutation>, LeftInd, LeftInj>{}, LeftPos{})
                          ||
                          RightTensorType::isZero(InsertAt<PermuteSequence<MultiIndex<N, DefectSignature>, RightPermutation>, RightInd, RightInj>{}, RightPos{})));
        }
 
        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 the indices referring to the left Tensor
          using LeftPosInd =
            TransformSequence<RealPos, Seq<>, IndexFunctor, AcceptInputInRangeFunctor<std::size_t, 0, leftRank_> >;
          using LeftPos = SequenceSlice<RealPos, LeftPosInd>;
          using LeftInd = SequenceSlice<RealIndices, LeftPosInd>;
 
          // The total set of left index positions, contraction and taken from Indices...
          using LeftTotalPos = JoinedDefects<LeftSortedPos, LeftPos>;
          using LeftInj = JoinedInjections<LeftSortedPos, LeftPos>;
 
          // Extract the indices referring to the right Tensor
          using RightPosInd =
            TransformSequence<RealPos, Seq<>, IndexFunctor, AcceptInputInRangeFunctor<std::size_t, leftRank_, leftRank_ + rightRank_> >;
          using RightPos = TransformSequence<SequenceSlice<RealPos, RightPosInd>, Seq<>, OffsetFunctor<-(ssize_t)leftRank_> >;
          using RightInd = SequenceSlice<RealIndices, RightPosInd>;
 
          // The total set of right index positions, contraction and taken from Indices...
          using RightTotalPos = JoinedDefects<RightSortedPos, RightPos>;
          using RightInj = JoinedInjections<RightSortedPos, RightPos>;
 
          return isZeroExpander(LeftInd{}, LeftInj{}, LeftTotalPos{},
                                RightInd{}, RightInj{}, RightTotalPos{},
                                MakeIndexSequence<multiDim(DefectSignature{})>{});
        }
 
       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());
        }
 
      };
 
      template<class Seq1, class Seq2, class T1, class T2,
               std::enable_if_t<AreProperTensors<T1, T2>::value, int> = 0>
      constexpr decltype(auto) einsum(T1&& t1, T2&& t2)
      {
        return finalize(EinsumFunctor<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) einsum(T1&& t1, T2&& t2)
      {
        constexpr std::size_t numDims = CommonHead<typename TensorTraits<T1>::Signature,
                                                   typename TensorTraits<T2>::Signature>::value;
        using ContractPos = MakeIndexSequence<numDims>;
        return einsum<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<EinsumOperation<Seq0, Seq1, Dims> >, T0&& t0, T1&& t1)
      {
        DUNE_ACFEM_RECORD_OPTIMIZATION;
 
        return EinsteinSummation<T0, Seq0, T1, Seq1>(std::forward<T0>(t0), std::forward<T1>(t1));
      }
 
    } // NS Tensor
 
    using Tensor::operate;
    using Tensor::einsum;
 
 
  } // NS ACFem
 
  template<class LeftTensor, class LeftIndices, class RightTensor, class RightIndices>
  struct FieldTraits<ACFem::Tensor::EinsteinSummation<LeftTensor, LeftIndices, RightTensor, RightIndices> >
    : FieldTraits<ACFem::FloatingPointClosure<typename ACFem::FieldPromotion<LeftTensor, RightTensor>::Type> >
  {};
 
 
} // NS Dune
 
#endif // __DUNE_ACFEM_TENSORS_OPERATIONS_EINSUM_HH__