ContractReduce: fall back to value-returning conj when a tile has no conj_to

src/TiledArray/tile_op/contract_reduce.h

@@ -470,6 +470,32 @@ class ContractReduce : public ContractReduceBase<Result, Left, Right, Scalar> {
 
 };  // 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<void>)
+template <typename Result, typename... Factor>
+inline Result conj_finalize(Result& temp, const Factor&... factor) {
+  if constexpr (TiledArray::has_conj_to_v<Result&, const Factor&...>) {
+    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<Result, Left, Right,
   result_type operator()() const { return result_type(); }
 
   /// Post processing step
+
+  /// Requirements on \c Result. `perm()` is a runtime value, so BOTH branches
+  /// below are instantiated for every tile used here:
+  /// - the permuted branch always needs the value-returning \c conj --
+  ///   `conj(result, perm)`, or `conj(result, factor, perm)` when the factor
+  ///   carries a scale -- whether or not a permutation is ever applied. This
+  ///   mirrors the primary template, whose finalization likewise instantiates
+  ///   `Permute<Result, Result>` 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<Scalar>) {
-      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());

src/TiledArray/tile_op/tile_interface.h

@@ -730,6 +730,33 @@ using result_of_conj_t = decltype(conj(std::declval<T>()...));
 template <typename... T>
 using result_of_conj_to_t = decltype(conj_to(std::declval<T>()...));
 
+namespace detail {
+template <typename, typename... T>
+struct has_conj_to_helper : public std::false_type {};
+template <typename... T>
+struct has_conj_to_helper<std::void_t<result_of_conj_to_t<T...>>, 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 <typename... T>
+inline constexpr bool has_conj_to_v =
+    detail::has_conj_to_helper<void, T...>::value;
+
 // Generic element-wise unary operations
 // ---------------------------------------------
 

tests/tile_op_contract_reduce.cpp

@@ -23,6 +23,10 @@
  *
  */
 
+#include <complex>
+#include <vector>
+
+#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<double>;
+
+  ConjOnlyTile() = default;
+  explicit ConjOnlyTile(std::vector<value_type> data)
+      : data_(std::move(data)) {}
+
+  bool empty() const { return data_.empty(); }
+  const std::vector<value_type>& data() const { return data_; }
+
+  ConjOnlyTile conj() const {
+    std::vector<value_type> out;
+    out.reserve(data_.size());
+    for (auto const& z : data_) out.push_back(std::conj(z));
+    return ConjOnlyTile(std::move(out));
+  }
+  template <typename Scalar,
+            typename = std::enable_if_t<std::is_arithmetic_v<Scalar>>>
+  ConjOnlyTile conj(const Scalar factor) const {
+    std::vector<value_type> 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 <typename Perm, typename = std::enable_if_t<
+                               TiledArray::detail::is_permutation_v<Perm>>>
+  ConjOnlyTile conj(const Perm&) const {
+    return conj();
+  }
+  template <
+      typename Scalar, typename Perm,
+      typename = std::enable_if_t<std::is_arithmetic_v<Scalar> &&
+                                  TiledArray::detail::is_permutation_v<Perm>>>
+  ConjOnlyTile conj(const Scalar factor, const Perm&) const {
+    return conj(factor);
+  }
+
+ private:
+  std::vector<value_type> 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<std::complex<double>>;
+  BOOST_CHECK(TiledArray::has_conj_to_v<TensorZ&>);
+  BOOST_CHECK((TiledArray::has_conj_to_v<TensorZ&, const double&>));
+
+  using BtasZ = btas::Tensor<std::complex<double>, TiledArray::Range>;
+  BOOST_CHECK(TiledArray::has_conj_to_v<BtasZ&>);
+  BOOST_CHECK((TiledArray::has_conj_to_v<BtasZ&, const double&>));
+  static_assert(
+      !TiledArray::detail::has_member_function_conj_to_anyreturn_v<BtasZ&>,
+      "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<ArenaTensor<...>> is MPQC's ToT cell
+  // type and the one this dispatch must not quietly move off the in-place path
+  using ArenaToTZ = Tensor<TiledArray::ArenaTensor<std::complex<double>>>;
+  BOOST_CHECK(TiledArray::has_conj_to_v<ArenaToTZ&>);
+  BOOST_CHECK((TiledArray::has_conj_to_v<ArenaToTZ&, const double&>));
+  using OwnToTZ = Tensor<Tensor<std::complex<double>>>;
+  BOOST_CHECK(TiledArray::has_conj_to_v<OwnToTZ&>);
+
+  BOOST_CHECK(!TiledArray::has_conj_to_v<ConjOnlyTile&>);
+}
+
+BOOST_AUTO_TEST_CASE(conj_finalize_in_place_when_supported) {
+  const std::complex<double> z{1.0, 2.0};
+  Tensor<std::complex<double>> tile(Tensor<std::complex<double>>::range_type(
+      std::array<std::size_t, 1>{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<double> z0{1.0, 2.0}, z1{-3.0, 0.5};
+
+  using CR_conj = ContractReduce<ConjOnlyTile, ConjOnlyTile, ConjOnlyTile,
+                                 TiledArray::detail::ComplexConjugate<void>>;
+  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<ConjOnlyTile, ConjOnlyTile, ConjOnlyTile,
+                     TiledArray::detail::ComplexConjugate<double>>;
+  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()