dune-acfem/master/operations_2product_8hh_source.html

#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__

Dune::ACFem::Tensor::ProductTensor
Compute the component-wise product ".[i][j]" over selected indices:
Definition: product.hh:66

Dune::ACFem::Expressions::DontOptimize
OptimizeTag< 0 > DontOptimize
Bottom level is overloaded to do nothing.
Definition: optimizationbase.hh:74

Dune::ACFem::operationName
std::string operationName(F &&f, const std::string &arg)
Verbose print of an operation, helper function to produce noise.
Definition: operationtraits.hh:601

Dune::ACFem::IsTypedValue
BoolConstant< ExpressionTraits< T >::isTypedValue > IsTypedValue
Compile-time true if T is a "typed value", e.g. a std::integral_constant.
Definition: expressiontraits.hh:90

Dune::ACFem::HeadPart
typename GetHeadPartHelper< Cnt, Seq >::Type HeadPart
Extract Cnt many consecutive elements from the front of Seq.
Definition: access.hh:217

Dune::ACFem::size
constexpr std::size_t size()
Gives the number of elements in tuple-likes and std::integer_sequence.
Definition: size.hh:73

Dune::ACFem::CommonHead
IndexConstant< CommonHeadHelper< 0UL, Seq1, Seq2 >::value > CommonHead
Compute the number of identical indices at the head of the sequence.
Definition: filter.hh:173

Dune::ACFem::MakeIndexSequence
MakeSequence< std::size_t, N, Offset, Stride, Repeat > MakeIndexSequence
Make a sequence of std::size_t elements.
Definition: generators.hh:34

Dune::ACFem::OffsetFunctor
MultiplyOffsetFunctor< 1, Offset > OffsetFunctor
Offset functor.
Definition: transform.hh:335

Dune::ACFem::SequenceCat
typename SequenceCatHelper2< S... >::Type SequenceCat
Concatenate the given sequences, in order, to <S0, S1, ... >.
Definition: transform.hh:220

Dune::ACFem::InsertAt
typename InsertAtHelper< Sequence< typename Input::value_type >, Input, Inject, Pos, AssumeSorted >::Type InsertAt
Insert Inject into the sequence Input at the position specified by Pos.
Definition: insertat.hh:87

Dune::ACFem::insertAt
auto insertAt(SrcTuple &&src, DataTuple &&data, BoolConstant< AssumeSorted >=BoolConstant< AssumeSorted >{})
Insert the elements of data at positions given by pos into src in turn.
Definition: subtuple.hh:106

Dune::ACFem::forwardSubTuple
auto forwardSubTuple(const TupleLike &t, IndexSequence< I... >=IndexSequence< I... >{})
Like subTuble() but forward the arguments, possibly in order to be expanded as parameters to another ...
Definition: subtuple.hh:162

Dune::ACFem::Tensor::multiply
constexpr decltype(auto) multiply(T1 &&t1, T2 &&t2)
Greedy tensor contraction for proper tensors, multiplying over all matching dimensions.
Definition: product.hh:322

Dune::ACFem::Tensor::ProductSignature
typename ProductOperationTraits< typename TensorTraits< Left >::TensorType, LPos, typename TensorTraits< Right >::TensorType, RPos >::Signature ProductSignature
Generate the signature of the ProductTensor when the operands have the given signature and the produc...
Definition: product.hh:93

Dune::ACFem::Tensor::ProductTensor< LeftTensor, Seq< LeftIndices... >, RightTensor, Seq< RightIndices... > >::ProductTensor
ProductTensor(Dummy &&...)
Allow default construction if contained types fulfill IsTypedValue.
Definition: product.hh:185

Dune::ACFem::Tensor::operate
constexpr auto operate(Expressions::DontOptimize, OperationTraits< TensorProductOperation< Seq0, Seq1, Dims > >, T0 &&t0, T1 &&t1)
Definition: product.hh:334

Dune::ACFem::Tensor::ProductTensor< LeftTensor, Seq< LeftIndices... >, RightTensor, Seq< RightIndices... > >::operator()
auto operator()(Dims... indices) const
Insert the current view-indices at their proper positions and foward to the underlying "host" tensor.
Definition: product.hh:197

Dune::ACFem::permute
constexpr auto permute(Seq, IndexSequence< Ps... >=IndexSequence< Ps... >{})
Permute the given sequence.
Definition: permutation.hh:116

Dune::ACFem::isSimple
constexpr bool isSimple(Sequence< T, V... > seq)
Definition: compare.hh:334

Dune::ACFem::IsIntegralPack
decltype(isIntegralPack(std::declval< T >()...)) IsIntegralPack
Decide whether the given parameter pack contains only integral types.
Definition: compare.hh:377

Dune::ACFem::PermuteSequence
std::decay_t< decltype(permute(Sequence{}, Perm{}))> PermuteSequence
Apply the given permutation to the positions of the given sequence.
Definition: permutation.hh:137

Dune::ACFem::ModelIntrospection::isZero
constexpr bool isZero
Shortcut identifying a zero model.
Definition: modeltraits.hh:642

std
STL namespace.

Dune::ACFem::TensorProductOperation
Component-wise product over given index-set.
Definition: expressionoperations.hh:389

Dune::ACFem::Tensor::TensorBase
Base class for all tensors.
Definition: tensorbase.hh:144