transpose.hh

#ifndef __DUNE_ACFEM_TENSORS_OPTIMIZATION_TRANSPOSE_HH__
#define __DUNE_ACFEM_TENSORS_OPTIMIZATION_TRANSPOSE_HH__
 
#include "../../expressions/optimization.hh"
#include "../../mpl/sequencesetoperations.hh"
#include "../operations/transpose.hh"
#include "../operations/transposetraits.hh"
#include "../modules.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    namespace Tensor::Optimization {
 
      namespace Transpose {
 
        using SubExprTag = Policy::SubExpressionOptimizationTag;
        using DefaultTag = Policy::DefaultOptimizationTag;
 
        template<class Permutation, class T, std::enable_if_t<IsTransposeExpression<T>::value, int> = 0>
        constexpr decltype(auto) operate(DefaultTag, F<TransposeOperation<Permutation> >, T&& t)
        {
          DUNE_ACFEM_RECORD_OPTIMIZATION;
 
          // permute the outer permutation with the inner
          using TType = std::decay_t<T>;
          using Perm = PermuteSequence<Permutation, typename TType::Permutation>;
 
          DUNE_ACFEM_EXPRESSION_RESULT(
            operate<TransposeOperation<Perm> >(std::forward<T>(t).operand(0_c))
            , "nested transpositions"
            );
        }
 
        template<class Permutation, class T,
                 std::enable_if_t<IsSelfTransposed<Permutation, T>::value, int> = 0>
        constexpr decltype(auto) operate(OptimizeTerminal1, F<TransposeOperation<Permutation> >, T&& t)
        {
          DUNE_ACFEM_RECORD_OPTIMIZATION;
 
          return forwardReturnValue<T>(t);
        }
 
        template<class Permutation, class T,
                 std::enable_if_t<(IsUnaryMinusExpression<T>::value
                                   && Permutation::size() == TensorTraits<T>::rank
          ), int> = 0>
        constexpr decltype(auto) operate(SubExprTag, F<TransposeOperation<Permutation> >, T&& t)
        {
          DUNE_ACFEM_RECORD_OPTIMIZATION;
 
          DUNE_ACFEM_EXPRESSION_RESULT(
            operate<MinusOperation>(
              operate<TransposeOperation<Permutation> >(std::forward<T>(t).operand(0_c))
              )
            , "self transpose"
            );
        }
 
        template<class Permutation, class T, class SFINAE = void>
        struct IsTransposeOfLeftProductWithScalar
          : BoolConstant<TensorTraits<Operand<0, T> >::rank == 0>
        {};
 
        template<class Permutation, class T>
        struct IsTransposeOfLeftProductWithScalar<
          Permutation, T,
          std::enable_if_t<!IsProductExpression<T>::value> >
          : FalseType
        {};
 
        template<class Permutation, class T,
                 std::enable_if_t<IsTransposeOfLeftProductWithScalar<Permutation, T>::value, int> = 0>
        constexpr decltype(auto) operate(SubExprTag, F<TransposeOperation<Permutation> >, T&& t)
        {
          DUNE_ACFEM_RECORD_OPTIMIZATION;
 
          DUNE_ACFEM_EXPRESSION_RESULT(
            operate(
              std::forward<T>(t).operation(),
              std::forward<T>(t).operand(0_c),
              operate<TransposeOperation<Permutation> >(
                std::forward<T>(t).operand(1_c)
                )
              )
            , "transpose inside product"
            );
        }
 
        template<class Permutation, class T,
                 std::enable_if_t<(IsConstKroneckerDelta<T>::value
                                   && Permutation::size() <= TensorTraits<T>::rank
          ), int> = 0>
        constexpr decltype(auto) operate(OptimizeTerminal0, F<TransposeOperation<Permutation> >, T&& t)
        {
          DUNE_ACFEM_RECORD_OPTIMIZATION;
 
          using TType = std::decay_t<T>;
          using Inverse = InversePermutation<Permutation>;
          using Signature = PermuteSequence<typename TType::Signature, Inverse>;
          using Pivots = PermuteSequence<typename std::decay_t<T>::PivotSequence, Inverse>;
          using ValueType = typename TType::ValueType;
 
          return KroneckerDelta<Signature, Pivots, ValueType>{};
        }
 
        template<class T, class Perm, class SFINAE = void>
        struct EinsumCommutationHelper
        {
          using Einsum = std::decay_t<T>;
          static constexpr std::size_t leftRank_ = TensorTraits<Operand<0, Einsum> >::rank - Einsum::defectRank_;
          static constexpr std::size_t rightRank_ = TensorTraits<Operand<1, Einsum> >::rank - Einsum::defectRank_;
          using LeftPos = MakeIndexSequence<leftRank_>;
          using RightPos = MakeIndexSequence<rightRank_, leftRank_>; // offset
 
          using LeftCommuted = MakeIndexSequence<leftRank_, rightRank_>;
          using RightCommuted = MakeIndexSequence<rightRank_>;
 
          using CommutationPermutation = SequenceCat<LeftCommuted, RightCommuted>;
          using InversePermutation = SequenceCat<RightPos, LeftPos>;
 
          static constexpr bool isCommutation = std::is_same<Perm, CommutationPermutation>::value;
 
          // we try to normalize permutations if the operands are
          // invariant w.r.t. to the permutation which results in
          // ascending permutation indices.
          using LeftPermutation = SequenceSlice<Perm, LeftPos>;
          using RightPermutation = SequenceSlice<Perm, RightPos>;
          using LeftNormalization = typename SortSequence<LeftPermutation>::Permutation;
          using RightNormalization = typename SortSequence<RightPermutation>::Permutation;
 
          static constexpr bool leftNormalizable =
            IsSelfTransposed<LeftNormalization, Operand<0, Einsum> >::value
            &&
            !isSimple(LeftNormalization{});
          static constexpr bool rightNormalizable =
            IsSelfTransposed<RightNormalization, Operand<1, Einsum> >::value
            &&
            !isSimple(RightNormalization{});
 
          static constexpr bool isNormalizable = (leftNormalizable || rightNormalizable);
          using Normalization =
            SequenceCat<ConditionalType<leftNormalizable, LeftNormalization, LeftPos>,
                        ConditionalType<rightNormalizable, OffsetSequence<leftRank_, RightNormalization>, RightPos> >;
 
#ifndef NDEBUG
          using LeftTransposed = PermuteSequenceValues<LeftPos, Perm>;
          using RightTransposed = PermuteSequenceValues<RightPos, Perm>;
 
          static_assert(isCommutation
                        ==
                        (std::is_same<LeftTransposed, LeftCommuted>::value
                         &&
                         std::is_same<RightTransposed, RightCommuted>::value),
                        "Inconsistent commutation permutation.");
#endif
        };
 
        template<class A, class B>
        struct EinsumCommutationHelper<A, B, std::enable_if_t<!IsEinsumExpression<A>::value> >
        {};
 
        template<class Permutation, class T, class SFINAE = void>
        struct IsCommutationOfEinsumTranspose
          : BoolConstant<(!IsSelfTransposed<Permutation, T>::value
                          && (IsEye<Operand<1, T> >::value // don't undo move eye left
                              ||
                              !IsEye<Operand<0, T> >::value)
                          && EinsumCommutationHelper<T, Permutation>::isCommutation)>
        {};
 
        template<class Permutation, class T>
        struct IsCommutationOfEinsumTranspose<
          Permutation, T,
          std::enable_if_t<!IsEinsumExpression<T>::value> >
          : FalseType
        {};
 
        template<class Permutation, class T,
                 std::enable_if_t<IsCommutationOfEinsumTranspose<Permutation, T>::value, int> = 0>
        constexpr decltype(auto) operate(SubExprTag, F<TransposeOperation<Permutation> >, T&& t)
        {
          DUNE_ACFEM_RECORD_OPTIMIZATION;
 
          DUNE_ACFEM_EXPRESSION_RESULT(
            operate<Operation<T> >(
              std::forward<T>(t).operand(1_c),
              std::forward<T>(t).operand(0_c)
              )
            , "transpose commutation"
            );
        }
 
        template<class Permutation, class T, class SFINAE = void>
        struct IsTransposeOfEinsumWithRightEye
          : BoolConstant<(!IsSelfTransposed<Permutation, T>::value
                          && !IsEye<Operand<0, T> >::value // no point in doing so
                          && IsEye<Operand<1, T> >::value
                          && !EinsumCommutationHelper<T, Permutation>::isCommutation)>
        {};
 
        template<class Permutation, class T>
        struct IsTransposeOfEinsumWithRightEye<
          Permutation, T,
          std::enable_if_t<!IsEinsumExpression<T>::value> >
          : FalseType
        {};
 
        template<class Permutation, class T,
                 std::enable_if_t<IsTransposeOfEinsumWithRightEye<Permutation, T>::value, int> = 0>
        constexpr decltype(auto) operate(SubExprTag, F<TransposeOperation<Permutation> >, T&& t)
        {
          DUNE_ACFEM_RECORD_OPTIMIZATION;
 
          using Traits = EinsumCommutationHelper<T, Permutation>;
 
          DUNE_ACFEM_EXPRESSION_RESULT(
            operate<TransposeOperation<Permutation> >(
              operate(
                DontOptimize{},
                F<TransposeOperation<typename Traits::InversePermutation> >{},
                operate(
                  std::forward<T>(t).operation(),
                  std::forward<T>(t).operand(1_c),
                  std::forward<T>(t).operand(0_c)
                  )))
            , "transpose normalize eyes"
            );
        }
 
        template<class Permutation, class T, class SFINAE = void>
        struct IsNormalizableTranspose
          : BoolConstant<EinsumCommutationHelper<T, Permutation>::isNormalizable>
        {};
 
        template<class Permutation, class T>
        struct IsNormalizableTranspose<
          Permutation, T,
          std::enable_if_t<(IsSelfTransposed<Permutation, T>::value
                            || !IsEinsumExpression<T>::value
          )> >
          : FalseType
        {};
 
        template<class Permutation, class T,
                 std::enable_if_t<IsNormalizableTranspose<Permutation, T>::value, int> = 0>
        constexpr decltype(auto) operate(DefaultTag, F<TransposeOperation<Permutation> >, T&& t)
        {
          DUNE_ACFEM_RECORD_OPTIMIZATION;
 
          using Traits = EinsumCommutationHelper<T, Permutation>;
 
          DUNE_ACFEM_EXPRESSION_RESULT(
            operate<TransposeOperation<Permutation> >(
              operate(DontOptimize{},
                      F<TransposeOperation<typename Traits::Normalization> >{},
                      std::forward<T>(t)
                ))
            , "sort transpositions"
            );
        }
 
 
        template<class F, class T0, class T1, class SFINAE = void>
        constexpr inline bool IsLeftSelfTransposedContractionV = HasInvariantValuesV<
          MPL::SequenceSetMinus<MakeIndexSequence<TensorTraits<T0>::rank>, typename EinsumTraits<F>::LeftIndexPositions>,
          typename TransposeTraits<T0>::Permutation>;
 
        template<class F, class T0, class T1>
        constexpr inline bool IsLeftSelfTransposedContractionV<
          F, T0, T1,
          std::enable_if_t<(!FunctorHas<IsEinsumOperation, F>::value
                            || std::is_same<F, ScalarEinsumFunctor>::value
                            || !IsTransposeExpression<T0>::value
          )> > = false;
 
        template<class F, class T0, class T1,
                 std::enable_if_t<IsLeftSelfTransposedContractionV<F, T0, T1>, int> = 0>
        constexpr decltype(auto) operate(DefaultTag, F, T0&& t0, T1&& t1)
        {
          DUNE_ACFEM_RECORD_OPTIMIZATION;
 
          using Inverse = InversePermutation<typename TransposeTraits<T0>::Permutation>;
          using Traits = EinsumTraits<F>;
          using Operation = EinsumOperation<
            PermuteSequenceValues<typename Traits::LeftIndexPositions, Inverse>,
            typename Traits::RightIndexPositions,
            typename Traits::Dimensions>;
 
          DUNE_ACFEM_EXPRESSION_RESULT(
            operate<Operation>(
              std::forward<T0>(t0).operand(0_c),
              std::forward<T1>(t1))
            , "left self transposed contraction"
            );
        }
 
 
        template<class F, class T0, class T1, class SFINAE = void>
        constexpr inline bool IsRightSelfTransposedContractionV = HasInvariantValuesV<
          MPL::SequenceSetMinus<MakeIndexSequence<TensorTraits<T1>::rank>, typename EinsumTraits<F>::LeftIndexPositions>,
          typename TransposeTraits<T1>::Permutation>;
 
        template<class F, class T0, class T1>
        constexpr inline bool IsRightSelfTransposedContractionV<
          F, T0, T1,
          std::enable_if_t<(!FunctorHas<IsEinsumOperation, F>::value
                            || std::is_same<F, ScalarEinsumFunctor>::value
                            || !IsTransposeExpression<T1>::value
          )> > = false;
 
        template<class F, class T0, class T1,
                 std::enable_if_t<IsRightSelfTransposedContractionV<F, T0, T1>, int> = 0>
        constexpr decltype(auto) operate(OptimizeNext<DefaultTag>, F, T0&& t0, T1&& t1)
        {
          DUNE_ACFEM_RECORD_OPTIMIZATION;
 
          using Inverse = InversePermutation<typename TransposeTraits<T1>::Permutation>;
          using Traits = EinsumTraits<F>;
          using Operation = EinsumOperation<
            typename Traits::LeftIndexPositions,
            PermuteSequenceValues<typename Traits::RightIndexPositions, Inverse>,
            typename Traits::Dimensions>;
 
          DUNE_ACFEM_EXPRESSION_RESULT(
            operate<Operation>(
              std::forward<T0>(t0),
              std::forward<T1>(t1).operand(0_c)
              )
            , "right self transposed contraction"
            );
        }
 
 
      } // Transpose::
 
    } // Tensor::Optimization::
 
    namespace Expressions {
 
      using Tensor::Optimization::Transpose::operate;
 
    } // Expressions
 
  } // ACFem
 
} // Dune
 
#endif // __DUNE_ACFEM_TENSORS_OPTIMIZATION_TRANSPOSE_HH__