diff --git a/src/TiledArray/tile_op/contract_reduce.h b/src/TiledArray/tile_op/contract_reduce.h index 2cf556c298..df1e7e2727 100644 --- a/src/TiledArray/tile_op/contract_reduce.h +++ b/src/TiledArray/tile_op/contract_reduce.h @@ -470,6 +470,32 @@ class ContractReduce : public ContractReduceBase { }; // class ContractReduce +/// Conjugate (and scale) a finished contraction result. + +/// In place through `conj_to` when the result tile supports it -- every tile +/// in the tree does, and it saves an allocate/copy/free per tile -- otherwise +/// through the value-returning `conj`. The two are independent tile-interface +/// customization points and a tile may implement only the latter, in which +/// case an unpermuted `conj(A*B)` used to fail to instantiate while the +/// permuted form (which already goes through `conj`) compiled. See issue #585. +/// +/// N.B. the choice is made on the ADL call `conj_to(...)`, not on a +/// `conj_to()` member: `btas::Tensor` has no member, only free functions in +/// namespace `btas`, and a member-based test would demote it to the +/// allocating path. +/// \tparam Result the result tile type +/// \tparam Factor the scale, if any (none for a ComplexConjugate) +template +inline Result conj_finalize(Result& temp, const Factor&... factor) { + if constexpr (TiledArray::has_conj_to_v) { + using TiledArray::conj_to; + return conj_to(temp, factor...); + } else { + using TiledArray::conj; + return conj(temp, factor...); + } +} + /// Contract and (sum) reduce operation with a ComplexConjugate factor /// The contraction of \c conj(A*B) (\c Scalar = \c void) or of a scaled @@ -564,22 +590,39 @@ class ContractReduce` unconditionally. + /// - the unpermuted branch takes \c conj_to when the tile has it and falls + /// back to \c conj otherwise (see \c conj_finalize), so \c conj_to is + /// optional -- but preferred, since it saves an allocate/copy/free. + /// + /// A tile missing one of the \c conj overloads fails inside + /// tile_interface.h's default `conj`, reporting only "too many arguments to + /// function call". That is not turned into a named static_assert because the + /// default \c conj CPOs are constrained only on \c Perm being a permutation + /// and have a DEDUCED return type, so `decltype(conj(arg, perm))` has to + /// instantiate the body: detecting them hard-errors instead of yielding + /// false, unlike \c has_conj_to_v (whose CPO is constrained on the member). + /// Making that detectable means constraining the four \c conj overloads the + /// way \c conj_to and \c neg_to already are -- worth doing, but a change to + /// a public header's overload set that wants its own PR. result_type operator()(result_type& temp) const { using TiledArray::empty; TA_ASSERT(!empty(temp)); if constexpr (std::is_void_v) { - if (!ContractReduceBase_::perm()) { - using TiledArray::conj_to; - return conj_to(temp); - } + if (!ContractReduceBase_::perm()) return conj_finalize(temp); using TiledArray::conj; return conj(temp, ContractReduceBase_::perm()); } else { - if (!ContractReduceBase_::perm()) { - using TiledArray::conj_to; - return conj_to(temp, ContractReduceBase_::factor().factor()); - } + if (!ContractReduceBase_::perm()) + return conj_finalize(temp, ContractReduceBase_::factor().factor()); using TiledArray::conj; return conj(temp, ContractReduceBase_::factor().factor(), ContractReduceBase_::perm()); diff --git a/src/TiledArray/tile_op/tile_interface.h b/src/TiledArray/tile_op/tile_interface.h index 942953e869..c7a3e214b2 100644 --- a/src/TiledArray/tile_op/tile_interface.h +++ b/src/TiledArray/tile_op/tile_interface.h @@ -730,6 +730,33 @@ using result_of_conj_t = decltype(conj(std::declval()...)); template using result_of_conj_to_t = decltype(conj_to(std::declval()...)); +namespace detail { +template +struct has_conj_to_helper : public std::false_type {}; +template +struct has_conj_to_helper>, T...> + : public std::true_type {}; +} // namespace detail + +/// Whether the ADL call `conj_to(args...)` is viable, i.e. whether a tile +/// supports IN-PLACE conjugation. + +/// `conj` and `conj_to` are independent customization points: a tile may +/// implement the value-returning `conj` and not the in-place `conj_to`, which +/// is a legal partial implementation of the tile interface. Consumers that can +/// use either should prefer `conj_to` when this is true (it saves an +/// allocate/copy/free per tile) and fall back to `conj` otherwise. +/// +/// N.B. this tests the ADL CALL, not a `conj_to()` member. `btas::Tensor` has +/// no such member -- only free functions in namespace `btas` -- so a +/// member-based test (`detail::has_member_function_conj_to_anyreturn_v`) would +/// report false for it and silently push it onto the allocating path. +/// \tparam T the argument types of the call, e.g. `Result&` or +/// `Result&, const Scalar&` +template +inline constexpr bool has_conj_to_v = + detail::has_conj_to_helper::value; + // Generic element-wise unary operations // --------------------------------------------- diff --git a/tests/tile_op_contract_reduce.cpp b/tests/tile_op_contract_reduce.cpp index b50397097d..7c574f35f3 100644 --- a/tests/tile_op_contract_reduce.cpp +++ b/tests/tile_op_contract_reduce.cpp @@ -23,6 +23,10 @@ * */ +#include +#include + +#include "TiledArray/external/btas.h" #include "TiledArray/tile_op/contract_reduce.h" #include "tiledarray.h" #include "unit_test_config.h" @@ -512,4 +516,138 @@ BOOST_AUTO_TEST_CASE(tensor_contract2) { BOOST_CHECK_EQUAL(result_map, C); } +// --------------------------------------------------------------------------- +// conj finalization on a tile that provides only the value-returning conj. +// +// `conj` and `conj_to` are independent tile-interface customization points, so +// implementing just the former is a legal partial implementation. Before +// issue #585 the ComplexConjugate finalization called `conj_to` +// unconditionally on the unpermuted path, so such a tile compiled for +// `c("k,i") = conj(a*b)` (permuted -> conj) and hard-errored for +// `c("i,k") = conj(a*b)` (unpermuted -> conj_to). +// --------------------------------------------------------------------------- + +namespace { + +/// A tile with the full value-returning `conj` overload set and NO `conj_to`, +/// member or free. Carries only what the finalization touches. +class ConjOnlyTile { + public: + using value_type = std::complex; + + ConjOnlyTile() = default; + explicit ConjOnlyTile(std::vector data) + : data_(std::move(data)) {} + + bool empty() const { return data_.empty(); } + const std::vector& data() const { return data_; } + + ConjOnlyTile conj() const { + std::vector out; + out.reserve(data_.size()); + for (auto const& z : data_) out.push_back(std::conj(z)); + return ConjOnlyTile(std::move(out)); + } + template >> + ConjOnlyTile conj(const Scalar factor) const { + std::vector out; + out.reserve(data_.size()); + for (auto const& z : data_) out.push_back(std::conj(z) * factor); + return ConjOnlyTile(std::move(out)); + } + // the suite never permutes; these exist so the permuted branch of the + // finalization instantiates alongside the unpermuted one under test + template >> + ConjOnlyTile conj(const Perm&) const { + return conj(); + } + template < + typename Scalar, typename Perm, + typename = std::enable_if_t && + TiledArray::detail::is_permutation_v>> + ConjOnlyTile conj(const Scalar factor, const Perm&) const { + return conj(factor); + } + + private: + std::vector data_; +}; + +} // namespace + +BOOST_AUTO_TEST_CASE(conj_to_detection) { + // every tile in the tree supports in-place conjugation, with and without a + // scale -- including btas::Tensor, whose conj_to is a free function in + // namespace btas rather than a member. That is why the trait tests the ADL + // CALL: a member-based test would report false for btas::Tensor and demote + // it to the allocating path. + using TensorZ = Tensor>; + BOOST_CHECK(TiledArray::has_conj_to_v); + BOOST_CHECK((TiledArray::has_conj_to_v)); + + using BtasZ = btas::Tensor, TiledArray::Range>; + BOOST_CHECK(TiledArray::has_conj_to_v); + BOOST_CHECK((TiledArray::has_conj_to_v)); + static_assert( + !TiledArray::detail::has_member_function_conj_to_anyreturn_v, + "btas::Tensor has no conj_to MEMBER; if this ever gains one the ADL " + "call is still the right test, but this guard has lost its point"); + + // the nested tile types too -- Tensor> is MPQC's ToT cell + // type and the one this dispatch must not quietly move off the in-place path + using ArenaToTZ = Tensor>>; + BOOST_CHECK(TiledArray::has_conj_to_v); + BOOST_CHECK((TiledArray::has_conj_to_v)); + using OwnToTZ = Tensor>>; + BOOST_CHECK(TiledArray::has_conj_to_v); + + BOOST_CHECK(!TiledArray::has_conj_to_v); +} + +BOOST_AUTO_TEST_CASE(conj_finalize_in_place_when_supported) { + const std::complex z{1.0, 2.0}; + Tensor> tile(Tensor>::range_type( + std::array{1ul})); + tile.at_ordinal(0) = z; + const auto result = TiledArray::detail::conj_finalize(tile); + BOOST_CHECK(result.at_ordinal(0) == std::conj(z)); + // in place: the argument itself was conjugated, no copy was made + BOOST_CHECK(tile.at_ordinal(0) == std::conj(z)); + + tile.at_ordinal(0) = z; + const auto scaled = TiledArray::detail::conj_finalize(tile, 2.0); + BOOST_CHECK(scaled.at_ordinal(0) == 2.0 * std::conj(z)); + BOOST_CHECK(tile.at_ordinal(0) == 2.0 * std::conj(z)); +} + +// The regression: instantiating and running the ComplexConjugate finalization +// for a tile with no conj_to. On master this is +// "contract_reduce.h: no matching function for call to 'conj_to'". +BOOST_AUTO_TEST_CASE(conj_finalization_without_conj_to) { + const std::complex z0{1.0, 2.0}, z1{-3.0, 0.5}; + + using CR_conj = ContractReduce>; + CR_conj op(math::blas::NoTranspose, math::blas::NoTranspose, + TiledArray::detail::conj_op(), 2u, 2u, 2u); + ConjOnlyTile temp({z0, z1}); + const ConjOnlyTile out = op(temp); + BOOST_REQUIRE_EQUAL(out.data().size(), 2u); + BOOST_CHECK(out.data()[0] == std::conj(z0)); + BOOST_CHECK(out.data()[1] == std::conj(z1)); + + using CR_scaled = + ContractReduce>; + CR_scaled op_scaled(math::blas::NoTranspose, math::blas::NoTranspose, + TiledArray::detail::conj_op(2.0), 2u, 2u, 2u); + ConjOnlyTile temp_scaled({z0, z1}); + const ConjOnlyTile out_scaled = op_scaled(temp_scaled); + BOOST_REQUIRE_EQUAL(out_scaled.data().size(), 2u); + BOOST_CHECK(out_scaled.data()[0] == 2.0 * std::conj(z0)); + BOOST_CHECK(out_scaled.data()[1] == 2.0 * std::conj(z1)); +} + BOOST_AUTO_TEST_SUITE_END()