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/max.hpp>
#include <prc/core/functors/view.hpp>
#include <prc/core/tensor/declarations.hpp>
#include <prc/core/tensor/views/contract.hpp>
 
namespace pRC
{
    template<Index... Is, class X, IsTensorish R = RemoveReference<X>>
        requires(isEven(sizeof...(Is)) && sizeof...(Is) <= R::Dimension) &&
        (max(Is...) < R::Dimension) &&
        (cut<2, 0>(select<Is...>(typename R::Sizes())) ==
            cut<2, 1>(select<Is...>(typename R::Sizes())))
    static inline constexpr auto contract(X &&a)
    {
        if constexpr(IsInvocable<View, X>)
        {
            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)
                {
                    using Sizes =
                        decltype(chip<S1::value(seq)..., S2::value(seq)...>(
                            typename R::Sizes()));
 
                    return TensorViews::Contract<T, Sizes, S1, S2, V>(
                        view(forward<X>(a)));
                });
        }
        else
        {
            return eval(contract<Is...>(a));
        }
    }
 
    template<Index... Is, class XA, class XB,
        IsTensorish RA = RemoveReference<XA>,
        IsTensorish RB = RemoveReference<XB>>
        requires(isEven(sizeof...(Is)) &&
                    sizeof...(Is) / 2 <= min(RA::Dimension, RB::Dimension)) &&
        (reduce<Max>(cut<2, 0>(Sizes<Is...>())) < RA::Dimension &&
            reduce<Max>(cut<2, 1>(Sizes<Is...>())) < RB::Dimension) &&
        (expand(cut<2, 0>(Sizes<Is...>()),
             [](auto const... indices)
             {
                 return select<indices...>(typename RA::Sizes());
             }) ==
            expand(cut<2, 1>(Sizes<Is...>()),
                [](auto const... indices)
                {
                    return select<indices...>(typename RB::Sizes());
                }))
    static inline constexpr auto contract(XA &&a, XB &&b)
    {
        if constexpr(IsInvocable<View, XA> && IsInvocable<View, XB>)
        {
            using SA = decltype(cut<2, 0>(Sequence<Index, Is...>()));
            using SB = decltype(cut<2, 1>(Sequence<Index, Is...>()) +
                Constant<Index, RA::Dimension>());
 
            return expand((SA(), SB()),
                [&a, &b](auto const... indices)
                {
                    return contract<indices...>(
                        tensorProduct(forward<XA>(a), forward<XB>(b)));
                });
        }
        else
        {
            return eval(contract<Is...>(a, b));
        }
    }
}
#endif // pRC_CORE_TENSOR_FUNCTIONS_CONTRACT_H