contract.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_CORE_TENSOR_FUNCTIONS_CONTRACT_H
#define pRC_CORE_TENSOR_FUNCTIONS_CONTRACT_H
 
#include <prc/core/basic/functions/is_even.hpp>
#include <prc/core/basic/functions/min.hpp>
#include <prc/core/functors/view.hpp>
#include <prc/core/tensor/views/contract.hpp>
 
namespace pRC
{
    template<Index... Is>
    struct Contract;
 
    template<Index... Is, class X, class R = RemoveReference<X>,
        If<IsTensorish<R>> = 0, If<IsInvocable<View, X>> = 0,
        If<IsSatisfied<(isEven(sizeof...(Is)))>> = 0,
        If<IsSatisfied<(sizeof...(Is) <= typename R::Dimension())>> = 0>
    static inline constexpr auto contract(X &&a)
    {
        using V = RemoveReference<decltype(view(forward<X>(a)))>;
 
        using E = ResultOf<X, typename R::Subscripts>;
        using T = RemoveConstReference<ResultOf<Add, E, E>>;
 
        using S1 = decltype(cut<2, 0>(Sequence<Index, Is...>()));
        using S2 = decltype(cut<2, 1>(Sequence<Index, Is...>()));
 
        return expand(makeSeries<Index, sizeof...(Is) / 2>(),
            [&a](auto const... seq)
            {
                static_assert(select<S1::value(seq)...>(typename R::Sizes()) ==
                        select<S2::value(seq)...>(typename R::Sizes()),
                    "Sizes of dimensions to be contracted differ.");
 
                using Sizes =
                    decltype(chip<S1::value(seq)..., S2::value(seq)...>(
                        typename R::Sizes()));
 
                return TensorViews::ContractUnary<T, Sizes, S1, S2, V>(
                    view(forward<X>(a)));
            });
    }
 
    template<Index... Is, class XA, class XB, class RA = RemoveReference<XA>,
        class RB = RemoveReference<XB>, If<IsTensorish<RA>> = 0,
        If<IsInvocable<View, XA>> = 0, If<IsTensorish<RB>> = 0,
        If<IsInvocable<View, XB>> = 0,
        If<IsSatisfied<(isEven(sizeof...(Is)))>> = 0,
        If<IsSatisfied<(sizeof...(Is) / 2 <=
            min(typename RA::Dimension(), typename RB::Dimension()))>> = 0>
    static inline constexpr auto contract(XA &&a, XB &&b)
    {
        using SA = decltype(cut<2, 0>(Sequence<Index, Is...>()));
        using SB = decltype(cut<2, 1>(Sequence<Index, Is...>()) +
            Constant<Index, typename RA::Dimension{}>());
 
        return expand((SA(), SB()),
            [&a, &b](auto const... indices)
            {
                return contract<indices...>(
                    tensorProduct(forward<XA>(a), forward<XB>(b)));
            });
    }
 
    template<Index... Is, class X, class R = RemoveReference<X>,
        If<IsTensorish<R>> = 0, If<Not<IsInvocable<View, X>>> = 0,
        If<IsInvocable<Contract<Is...>, X &>> = 0>
    static inline constexpr auto contract(X &&a)
    {
        return eval(contract<Is...>(a));
    }
 
    template<Index... Is, class XA, class XB, class RA = RemoveReference<XA>,
        class RB = RemoveReference<XB>, If<IsTensorish<RA>> = 0,
        If<IsTensorish<RB>> = 0,
        If<Not<All<IsInvocable<View, XA>, IsInvocable<View, XB>>>> = 0,
        If<IsInvocable<Contract<Is...>, XA &, XB &>> = 0>
    static inline constexpr auto contract(XA &&a, XB &&b)
    {
        return eval(contract<Is...>(a, b));
    }
}
#endif // pRC_CORE_TENSOR_FUNCTIONS_CONTRACT_H