ToT: conj() on nested-tile contractions; real-plain x complex-ToT products via one real gemm

src/TiledArray/tensor/arena_tensor.h

@@ -439,18 +439,52 @@ void fill(ArenaTensor<T, R>& dst, const U& value) {
   std::fill_n(dst.data(), dst.size(), static_cast<T>(value));
 }
 
-/// `dst *= factor`. No-op on a null view.
+/// `dst *= factor`. No-op on a null view. A ComplexConjugate<...> factor
+/// conjugates (and scales) every arena scalar in place.
 template <typename T, typename R, typename Scalar>
 void scale_to(ArenaTensor<T, R>& dst, Scalar factor) {
   if (!dst) return;
   auto* d = dst.data();
   const auto n = dst.size();
   // operator*= is the permissive one (std::complex<T>::operator*=(const T&)
-  // accepts any arithmetic factor), so no mixed-scalar detail operator is
+  // accepts any arithmetic factor; detail::operator*=(T&, ComplexConjugate)
+  // conjugates a scalar in place), so no mixed-scalar detail operator is
   // needed here
   for (std::size_t i = 0; i < n; ++i) d[i] *= factor;
 }
 
+/// `dst *= conj_op()`, `dst *= conj_op(S)`, `dst *= -conj_op()`: in-place
+/// conjugation (and scale) of every arena scalar. These overloads exist so
+/// that a view is never routed to the generic
+/// detail::operator*=(L&, ComplexConjugate), which needs a value-returning
+/// conj(L) that a non-owning view cannot provide.
+///
+/// N.B. all three are needed; the ComplexConjugate<S> one does NOT subsume
+/// the other two. Against detail's own
+/// operator*=(L&, const ComplexConjugate<void>&) it is ambiguous: this
+/// overload is more specialized in the first parameter, detail's in the
+/// second, so neither wins partial ordering. The <void> and <ComplexNegTag>
+/// overloads are more specialized in BOTH parameters and therefore win.
+template <typename T, typename R>
+ArenaTensor<T, R>& operator*=(ArenaTensor<T, R>& dst,
+                              const detail::ComplexConjugate<void>& factor) {
+  scale_to(dst, factor);
+  return dst;
+}
+template <typename T, typename R>
+ArenaTensor<T, R>& operator*=(
+    ArenaTensor<T, R>& dst,
+    const detail::ComplexConjugate<detail::ComplexNegTag>& factor) {
+  scale_to(dst, factor);
+  return dst;
+}
+template <typename T, typename R, typename S>
+ArenaTensor<T, R>& operator*=(ArenaTensor<T, R>& dst,
+                              const detail::ComplexConjugate<S>& factor) {
+  scale_to(dst, factor);
+  return dst;
+}
+
 /// `dst += src`. Asserts shape compatibility.
 /// A null `src` is an implicit zero and a no-op; a null `dst` with a populated
 /// `src` throws (a view has no storage to write into and cannot allocate, so
@@ -599,6 +633,9 @@ auto gemm(ArenaTensor<T, R>& result, const ArenaTensor<T, R>& left,
       (gemm_helper.right_op() == math::blas::NoTranspose) ? N : K;
   const integer ldc = N;
 
+  static_assert(!detail::is_complex_conjugate_v<Scalar>,
+                "ArenaTensor gemm: a ComplexConjugate<...> factor cannot be "
+                "applied inside the gemm; conjugate the finished result");
   math::blas::gemm(gemm_helper.left_op(), gemm_helper.right_op(), M, N, K,
                    static_cast<T>(factor), left.data(), lda, right.data(), ldb,
                    T(1), result.data(), ldc);

src/TiledArray/tensor/complex.h

@@ -220,6 +220,37 @@ inline ComplexConjugate<ComplexNegTag> ComplexConjugate<void>::operator-()
 template <typename S>
 struct is_numeric<ComplexConjugate<S>> : public std::true_type {};
 
+// Mixed-scalar products (an integral or floating scalar with a std::complex of
+// a different scalar type). Declared BEFORE the ComplexConjugate operators
+// below: those evaluate `conj(value) * op.factor()`, and an integer-literal
+// scale (`2 * conj(a*b)` -> ComplexConjugate<int>) resolves only if these are
+// visible at the template definition (they are not found by ADL).
+template <typename L, typename R,
+          typename = std::enable_if_t<std::is_integral_v<L>>>
+TILEDARRAY_FORCE_INLINE auto operator*(const L l, const std::complex<R> r) {
+  return static_cast<R>(l) * r;
+}
+
+template <typename L, typename R,
+          typename = std::enable_if_t<std::is_integral_v<R>>>
+TILEDARRAY_FORCE_INLINE auto operator*(const std::complex<L> l, const R r) {
+  return l * static_cast<L>(r);
+}
+
+template <typename L, typename R>
+TILEDARRAY_FORCE_INLINE
+    std::enable_if_t<std::is_floating_point_v<L>, std::complex<R>>
+    operator*(const L l, const std::complex<R> r) {
+  return std::complex<R>(l, 0.) * r;
+}
+
+template <typename L, typename R>
+TILEDARRAY_FORCE_INLINE
+    std::enable_if_t<std::is_floating_point_v<R>, std::complex<L>>
+    operator*(const std::complex<L> l, const R r) {
+  return l * std::complex<L>(r, 0.);
+}
+
 /// ComplexConjugate operator factory function
 
 /// \tparam S The scalar type
@@ -291,39 +322,48 @@ TILEDARRAY_FORCE_INLINE L& operator*=(L& value,
   return value;
 }
 
+template <typename T>
+struct is_complex_conjugate : std::false_type {};
+template <typename T>
+struct is_complex_conjugate<ComplexConjugate<T>> : std::true_type {};
+/// true if \c T is a ComplexConjugate<...> contraction factor
+template <typename T>
+inline constexpr bool is_complex_conjugate_v =
+    is_complex_conjugate<std::remove_cv_t<T>>::value;
+
+/// Whether inner cells of scalar type \c Inner can be viewed as arrays of
+/// \c Real for the interleaved real GEMM on complex cells: \c Inner is
+/// std::complex<Real> laid out as two contiguous \c Real, the standard's
+/// array-oriented access guarantee ([complex.numbers.general]/4). The size
+/// and alignment terms are part of the predicate rather than an assertion, so
+/// an ABI that breaks the layout falls back to the per-cell loop instead of a
+/// reinterpret_cast into undefined behavior.
+template <typename Inner, typename Real>
+inline constexpr bool is_interleaved_real_view_v =
+    std::is_same_v<std::complex<Real>, Inner> &&
+    sizeof(Inner) == 2 * sizeof(Real) && alignof(Inner) >= alignof(Real);
+
+/// The numeric multiplier to bake into a per-element (per-cell) multiply-add
+/// op for a contraction with factor \c factor: the factor itself (converted
+/// to \c Numeric) for a numeric factor, and \c Numeric(1) for a
+/// ComplexConjugate<...> factor -- the conjugation (and, for
+/// ComplexConjugate<Scalar>, the scale) of such a factor is applied to the
+/// finished result by ContractReduce's finalization step, not per element.
+template <typename Numeric, typename Scalar>
+TILEDARRAY_FORCE_INLINE Numeric elem_factor(const Scalar& factor) {
+  if constexpr (is_complex_conjugate_v<Scalar>)
+    return Numeric(1);
+  else
+    return static_cast<Numeric>(factor);
+}
+
 template <typename T>
 inline auto abs(const ComplexConjugate<T>& a) {
   return std::abs(a.factor());
 }
 
 inline int abs(const ComplexConjugate<void>& a) { return 1; }
 
-template <typename L, typename R,
-          typename = std::enable_if_t<std::is_integral_v<L>>>
-TILEDARRAY_FORCE_INLINE auto operator*(const L l, const std::complex<R> r) {
-  return static_cast<R>(l) * r;
-}
-
-template <typename L, typename R,
-          typename = std::enable_if_t<std::is_integral_v<R>>>
-TILEDARRAY_FORCE_INLINE auto operator*(const std::complex<L> l, const R r) {
-  return l * static_cast<L>(r);
-}
-
-template <typename L, typename R>
-TILEDARRAY_FORCE_INLINE
-    std::enable_if_t<std::is_floating_point_v<L>, std::complex<R>>
-    operator*(const L l, const std::complex<R> r) {
-  return std::complex<R>(l, 0.) * r;
-}
-
-template <typename L, typename R>
-TILEDARRAY_FORCE_INLINE
-    std::enable_if_t<std::is_floating_point_v<R>, std::complex<L>>
-    operator*(const std::complex<L> l, const R r) {
-  return l * std::complex<L>(r, 0.);
-}
-
 }  // namespace detail
 
 namespace conversions {

src/TiledArray/tensor/tensor.h

@@ -2062,14 +2062,41 @@ class Tensor {
     // early exit for empty this
     if (empty()) return *this;
 
-    if constexpr (is_arena_tensor_v<value_type>) {
-      // Arena inner cells: route through each cell's own in-place scale_to (the
-      // free arena kernel), which handles a ComplexConjugate factor by
-      // conjugating each arena scalar in place. Going through `cell *= factor`
-      // would instead select the generic operator*=(.., ComplexConjugate) ->
-      // detail::conj(cell), which has no value-returning conj for ArenaTensor.
-      return inplace_unary(
-          [factor](value_type& MADNESS_RESTRICT c) { c.scale_to(factor); });
+    if constexpr (detail::is_complex_conjugate_v<Scalar> &&
+                  !detail::is_complex_v<numeric_type>) {
+      // real elements: conj is the identity, only the scale (if any) remains
+      if constexpr (std::is_same_v<Scalar, detail::ComplexConjugate<void>>)
+        return *this;
+      else if constexpr (std::is_same_v<Scalar, detail::ComplexConjugate<
+                                                    detail::ComplexNegTag>>)
+        return neg_to();
+      else
+        return scale_to(factor.factor());
+    } else if constexpr (is_arena_tensor_v<value_type>) {
+      // Arena inner cells: the free arena kernel scales -- or, for a
+      // ComplexConjugate factor, conjugates and scales -- every arena scalar
+      // in place. (The member ArenaTensor::scale_to takes numeric factors
+      // only; `cell *= factor` is the same kernel via the ArenaTensor
+      // operator*= overloads.)
+      return inplace_unary([factor](value_type& MADNESS_RESTRICT c) {
+        ::TiledArray::scale_to(c, factor);
+      });
+    } else if constexpr (detail::is_complex_conjugate_v<Scalar> &&
+                         detail::is_ta_tensor_v<value_type>) {
+      // Owning nested cells with a ComplexConjugate factor: conjugate (and
+      // scale) each cell IN PLACE. `cell *= factor` would resolve to
+      // `cell = conj(cell) * S`, i.e. a fresh allocation, copy and free of
+      // every cell.
+      return inplace_unary([factor](value_type& MADNESS_RESTRICT res) {
+        if constexpr (std::is_same_v<Scalar, detail::ComplexConjugate<void>>)
+          res.conj_to();
+        else if constexpr (std::is_same_v<Scalar, detail::ComplexConjugate<
+                                                      detail::ComplexNegTag>>) {
+          res.conj_to();
+          res.neg_to();
+        } else
+          res.conj_to(factor.factor());
+      });
     } else {
       return inplace_unary(
           [factor](value_type& MADNESS_RESTRICT res) { res *= factor; });
@@ -3387,6 +3414,14 @@ class Tensor {
   Tensor& gemm(const Tensor<As...>& A, const Tensor<Bs...>& B, const W alpha,
                const math::GemmHelper& gemm_helper) {
     numeric_type beta = 1;
+    // A ComplexConjugate<...> factor is not a gemm alpha: conj does not
+    // distribute into the sum of products, so the expression layer applies it
+    // to the finished result (ContractReduce's finalization; see
+    // ContEngine::outer_factor) and hands the gemm a numeric alpha.
+    static_assert(!detail::is_complex_conjugate_v<W>,
+                  "Tensor::gemm: a ComplexConjugate<...> alpha cannot be "
+                  "applied inside the gemm; conjugate the finished result");
+    const numeric_type alpha_n = static_cast<numeric_type>(alpha);
     if (this->empty()) {
       *this =
           Tensor(gemm_helper.make_result_range<range_type>(A.range_, B.range()),
@@ -3411,7 +3446,7 @@ class Tensor {
       }
       for (size_t i = 0; i < this->nbatch(); ++i) {
         auto Ci = this->batch(i);
-        TiledArray::gemm(alpha, A.batch(i), B.batch(i),
+        TiledArray::gemm(alpha_n, A.batch(i), B.batch(i),
                          twostep ? numeric_type(0) : numeric_type(1), Ci,
                          gemm_helper);
       }
@@ -3465,7 +3500,7 @@ class Tensor {
 #else   // TA_ENABLE_TILE_OPS_LOGGING
     for (size_t i = 0; i < this->nbatch(); ++i) {
       auto Ci = this->batch(i);
-      TiledArray::detail::gemm(alpha, A.batch(i), B.batch(i), beta, Ci,
+      TiledArray::detail::gemm(alpha_n, A.batch(i), B.batch(i), beta, Ci,
                                gemm_helper);
     }
 #endif  // TA_ENABLE_TILE_OPS_LOGGING
@@ -3572,7 +3607,21 @@ class Tensor {
     if constexpr (detail::is_numeric_v<V> && is_tensor_view_v<U> &&
                   is_tensor_view_v<value_type>) {
       using Real = std::remove_cv_t<typename value_type::value_type>;
-      if constexpr (std::is_same_v<std::remove_cv_t<V>, Real>) {
+      using Vr = std::remove_cv_t<V>;
+      // Both slabs are viewed through Vr* below, so the left ToT's inner
+      // scalar must be the result's. Same element type: one gemm in that
+      // type. Real plain scalars (Vr) against complex<Vr> inner cells: the
+      // complex slabs are viewed as real matrices with the inner extent
+      // doubled (re,im interleaved), so one real gemm with alpha = beta = 1
+      // accumulates both parts exactly (detail::is_interleaved_real_view_v
+      // carries the layout terms; a mismatch falls back to the per-cell loop).
+      constexpr bool left_matches =
+          std::is_same_v<std::remove_cv_t<typename U::value_type>, Real>;
+      constexpr bool same_type = left_matches && std::is_same_v<Vr, Real>;
+      constexpr bool interleaved = left_matches && !std::is_same_v<Vr, Real> &&
+                                   detail::is_interleaved_real_view_v<Real, Vr>;
+      constexpr integer cw = interleaved ? 2 : 1;  // reals per inner element
+      if constexpr (same_type || interleaved) {
         // Both operand orientations are served: the plain (right) matrix's
         // transpose is passed to BLAS as its op, and a transposed nested
         // (left) tile only changes which cells form row m (cells (k,m) are
@@ -3685,6 +3734,9 @@ class Tensor {
                 if (detail::scale_gemm_timing_enabled()) {
                   detail::g_scale[0].gemm_runs.fetch_add(
                       1, std::memory_order_relaxed);
+                  // element multiply-adds (A in inner elements), the same
+                  // unit as fb_flop and the same-type path, so the
+                  // [scale-timing] coverage ratios compare like with like
                   detail::g_scale[0].gemm_flop.fetch_add(
                       2ull * static_cast<std::uint64_t>(K) *
                           static_cast<std::uint64_t>(N) *
@@ -3699,10 +3751,12 @@ class Tensor {
                     right_no_trans ? TiledArray::math::blas::Transpose
                                    : TiledArray::math::blas::NoTranspose,
                     TiledArray::math::blas::NoTranspose,
-                    /*M=*/N, /*N=*/Ai, /*K=*/K, Real(1),
+                    /*M=*/N, /*N=*/Ai * cw, /*K=*/K, Vr(1),
                     /*A=*/right_data, /*lda=*/right_no_trans ? N : K,
-                    /*B=*/lc(0).data(), /*ldb=*/ldb, Real(1),
-                    /*C=*/rc0[0].data(), /*ldc=*/ldc);
+                    /*B=*/reinterpret_cast<const Vr*>(lc(0).data()),
+                    /*ldb=*/ldb * cw, Vr(1),
+                    /*C=*/reinterpret_cast<Vr*>(rc0[0].data()),
+                    /*ldc=*/ldc * cw);
               } else {  // per-cell AXPY fallback for this row
                 if (detail::scale_gemm_timing_enabled()) {
                   // classify fallback reason (re-scan; observation only, does
@@ -3773,7 +3827,17 @@ class Tensor {
     if constexpr (detail::is_numeric_v<U> && is_tensor_view_v<V> &&
                   is_tensor_view_v<value_type>) {
       using Real = std::remove_cv_t<typename value_type::value_type>;
-      if constexpr (std::is_same_v<std::remove_cv_t<U>, Real>) {
+      using Ur = std::remove_cv_t<U>;
+      // see the tot_x_t block: the right ToT's inner scalar must be the
+      // result's; then same type, or real plain x complex<Ur> inner cells via
+      // the re,im-interleaved real gemm
+      constexpr bool right_matches =
+          std::is_same_v<std::remove_cv_t<typename V::value_type>, Real>;
+      constexpr bool same_type = right_matches && std::is_same_v<Ur, Real>;
+      constexpr bool interleaved = right_matches && !std::is_same_v<Ur, Real> &&
+                                   detail::is_interleaved_real_view_v<Real, Ur>;
+      constexpr integer cw = interleaved ? 2 : 1;  // reals per inner element
+      if constexpr (same_type || interleaved) {
         // Both orientations, as in the tot_x_t block: the plain (left)
         // matrix's transpose goes to BLAS as its op; a transposed nested
         // (right) tile stores cell (k,n) at n*K + k, so column n's cells are
@@ -3866,6 +3930,7 @@ class Tensor {
                 if (detail::scale_gemm_timing_enabled()) {
                   detail::g_scale[1].gemm_runs.fetch_add(
                       1, std::memory_order_relaxed);
+                  // element multiply-adds, the unit of fb_flop (see above)
                   detail::g_scale[1].gemm_flop.fetch_add(
                       2ull * static_cast<std::uint64_t>(M) *
                           static_cast<std::uint64_t>(K) *
@@ -3880,10 +3945,13 @@ class Tensor {
                     left_no_trans ? TiledArray::math::blas::NoTranspose
                                   : TiledArray::math::blas::Transpose,
                     TiledArray::math::blas::NoTranspose,
-                    /*M=*/M, /*N=*/Ai, /*K=*/K, Real(1),
+                    /*M=*/M, /*N=*/Ai * cw, /*K=*/K, Ur(1),
                     /*A=*/left_data, /*lda=*/left_no_trans ? K : M,
-                    /*B=*/right_data[rcell(0, n)].data(), /*ldb=*/sbc, Real(1),
-                    /*C=*/this_data[n].data(), /*ldc=*/ldc);
+                    /*B=*/
+                    reinterpret_cast<const Ur*>(right_data[rcell(0, n)].data()),
+                    /*ldb=*/sbc * cw, Ur(1),
+                    /*C=*/reinterpret_cast<Ur*>(this_data[n].data()),
+                    /*ldc=*/ldc * cw);
               } else {  // per-cell AXPY fallback for this column
                 if (detail::scale_gemm_timing_enabled()) {
                   // classify fallback reason (re-scan; observation only) +

src/TiledArray/tile_op/batched_contract_reduce.h

@@ -119,11 +119,17 @@ class BatchedContractReduce {
   /// Create a new, empty result object
   result_type operator()() const { return result_type(); }
 
-  /// Post processing step (no result permutation supported)
-  result_type operator()(const result_type& temp) const {
+  /// Post processing step (no result permutation supported): the wrapped
+  /// op's own finalization -- a no-op for a numeric factor, conjugation (and
+  /// scale) of the finished tile for a ComplexConjugate<...> one
+  result_type operator()(result_type& temp) const {
     using TiledArray::empty;
     TA_ASSERT(!empty(temp));
-    return temp;
+    // the wrapped op's finalization also applies ITS permutation; the batched
+    // op is built perm-free (ContEngine), and a perm leaking in here would
+    // permute the fused result silently
+    TA_ASSERT(!op_.perm());
+    return op_(temp);
   }
 
   /// Reduce two result objects (both carry the full fused range)