Tensor: an empty-operand fast path must own its copy, not alias the source

`detail::clone_or_cast<Tensor>(right)` deep-copies only when `Right` *is*
`Tensor`, where it calls `clone()`. For any other nested-tensor `Right` -- most
importantly the `TensorInterface<..., BlockRange>` that `Tensor::block()`
yields -- it falls into a generic branch that `std::copy`s the *cells*, and a
cell is a shallow handle: an `ArenaTensor` is a non-owning view, and
`TA::Tensor` is itself reference-counted. The copy therefore shared storage
with `right`, and the empty-`this` fast paths mutate it immediately afterwards:

    t.subt_to(p.block(...));        // neg_to()   -> p's cells negated
    t.add_to(p.block(...), 3.0);    // scale_to() -> p's cells scaled

Both corrupt a const source operand. The value-returning `add`/`subt`
empty-left exits have the same defect one step removed: they hand the caller a
result aliasing `right`, so the caller's first mutation writes through.

Route all eight of these through a new `owning_copy_of`, which for a nested
tensor builds the copy with `arena_trivial_unary` -- a fresh slab, elements
copied, no shared cells -- and keeps `clone_or_cast` for the flat and
same-type cases where it is already a deep copy. Being range-addressed, that
kernel also handles a strided source correctly.

This is reachable only because the null-cell guard earlier in this branch made
`add_to`/`subt_to` compile for a non-contiguous `Right` at all; previously they
failed to instantiate. The unscaled `add_to`'s aliasing predates that and is
fixed here too.

Also in this commit:

- `arena_trivial_scaled` joins its unary and binary siblings in addressing
  operands through `arena_cell_offset` and taking counts from the `total_size`
  CPO. It still used `nbatch()` and linear `data()[ord]`, so the mixed
  `Tensor<ArenaTensor>` +/-/* `Tensor<scalar>` overloads failed to compile
  against a `block()` operand under TA_ASSERT_THROW and silently read the wrong
  scalars under TA_ASSERT_IGNORE -- a TA_ASSERT_POLICY-dependent difference.

- einsum's view-inner ToT x T diagnostic now fires only once a source inner
  cell is actually populated. It threw on entry to the loop, which rejected the
  fully screened case that used to work: `axpy_to` early-returns on a null
  source, so an all-null operand never needed the destination to have storage
  and correctly produced a null result tile.

- A test pins intersection sparsity for the plain owning ToT. `Tensor::mult`'s
  intersection rule applies in the `is_ta_tensor_v` branch as well, so it is a
  behavior change for `TA::Tensor<TA::Tensor<T>>`, not only for arena tiles;
  the pre-existing mult_mismatched_null_inners is written permissively and
  passes under either rule.

Two of the six new corruption regressions were confirmed failing before the
fix (`-1 != 1` and `3 != 1`); the rest exercise the same path.

Author: evaleev

src/TiledArray/einsum/tiledarray.h

@@ -948,17 +948,35 @@ auto einsum(expressions::TsrExpr<ArrayA_> A, expressions::TsrExpr<ArrayB_> B,
               if constexpr (inner_is_view) {
                 // `tile` is built with default-constructed (null) inner cells
                 // and a view cell cannot allocate, so it has no storage to
-                // accept the first contribution: the free `axpy_to` refuses a
-                // null destination rather than silently drop the source. Only
+                // accept a contribution: the free `axpy_to` refuses a null
+                // destination rather than silently drop the source. Only
                 // run_regime_a_arena grows result cells, so reaching this
                 // legacy path means the arena plan was inactive (e.g. under
-                // detail::arena_disabled()). Report that here rather than let
-                // the per-cell kernel raise a lower-level diagnostic --
-                // mirrors the ToT x ToT branch above.
-                TA_EXCEPTION(
-                    "TA::einsum: ToT x T product with view inner cells (e.g. "
-                    "ArenaTensor) is supported only via the regime-A arena "
-                    "fast path, which was inactive for this expression");
+                // detail::arena_disabled()).
+                //
+                // Report that here rather than let the per-cell kernel raise a
+                // lower-level diagnostic -- but only once a source cell is
+                // actually populated. A fully screened operand contributes
+                // nothing, leaves the null result tile correct, and must keep
+                // working: `axpy_to` early-returns on a null source, so that
+                // case never needed the destination to have storage.
+                auto aik = ai.batch(k);
+                auto bik = bi.batch(k);
+                const auto vol = aik.total_size();
+                TA_ASSERT(vol == bik.total_size());
+                for (auto i = 0; i < vol; ++i) {
+                  bool src_populated;
+                  if constexpr (IsArrayToT<ArrayA>)
+                    src_populated = !aik.data()[i].empty();
+                  else
+                    src_populated = !bik.data()[i].empty();
+                  if (src_populated)
+                    TA_EXCEPTION(
+                        "TA::einsum: ToT x T product with view inner cells "
+                        "(e.g. ArenaTensor) is supported only via the regime-A "
+                        "arena fast path, which was inactive for this "
+                        "expression");
+                }
               } else {
                 auto aik = ai.batch(k);
                 auto bik = bi.batch(k);

src/TiledArray/tensor/arena_kernels.h

@@ -489,21 +489,26 @@ OuterTensor arena_trivial_scaled(const ToTSide& tot_outer,
                                  FillOp&& fill_op) {
   using elem_t = typename OuterTensor::value_type::value_type;
   using inner_range_t = typename OuterTensor::value_type::range_type;
+  // Either side may be a strided view; address both through `arena_cell_offset`
+  // and take the cell/batch counts from the `total_size` CPO, as the unary and
+  // binary kernels do. `nbatch()` is a TA::Tensor member a view does not have.
   TA_ASSERT(tot_outer.range().volume() == scalar_outer.range().volume());
-  TA_ASSERT(tot_outer.nbatch() == scalar_outer.nbatch());
+  TA_ASSERT(TiledArray::total_size(tot_outer) ==
+            TiledArray::total_size(scalar_outer));
   auto range_fn = [&tot_outer](std::size_t ord) -> inner_range_t {
-    const auto& t = tot_outer.data()[ord];
+    const auto& t = tot_outer.data()[arena_cell_offset(tot_outer, ord)];
     return t.empty() ? inner_range_t{} : t.range();
   };
-  OuterTensor result = arena_outer_init<OuterTensor>(
-      tot_outer.range(), tot_outer.nbatch(), range_fn, alignof(elem_t),
-      /*zero_init=*/false);
-  const std::size_t N_cells = tot_outer.range().volume() * tot_outer.nbatch();
+  const auto [N_cells, tot_nbatch] = arena_cells_and_batches(tot_outer);
+  OuterTensor result =
+      arena_outer_init<OuterTensor>(tot_outer.range(), tot_nbatch, range_fn,
+                                    alignof(elem_t), /*zero_init=*/false);
   for (std::size_t ord = 0; ord < N_cells; ++ord) {
     auto& dst = result.data()[ord];
     if (dst.empty()) continue;
-    fill_op(dst.data(), tot_outer.data()[ord].data(), scalar_outer.data()[ord],
-            dst.size());
+    fill_op(
+        dst.data(), tot_outer.data()[arena_cell_offset(tot_outer, ord)].data(),
+        scalar_outer.data()[arena_cell_offset(scalar_outer, ord)], dst.size());
   }
   return result;
 }

src/TiledArray/tensor/tensor.h

@@ -2167,6 +2167,39 @@ class Tensor {
     }
   }
 
+  /// An owning deep copy of \p right, safe to mutate without touching it.
+
+  /// `detail::clone_or_cast` deep-copies only when `Right` *is* `Tensor`, in
+  /// which case it calls `clone()`. For any other nested-tensor `Right` it
+  /// falls into a generic branch that `std::copy`s the *cells*, and a cell is
+  /// a shallow handle -- an `ArenaTensor` is a non-owning view, and
+  /// `TA::Tensor` is itself reference-counted. The copy would then share
+  /// storage with \p right, so negating or scaling it afterwards writes
+  /// through to the source, corrupting a const operand.
+  ///
+  /// Route nested tensors through the arena unary kernel instead: it allocates
+  /// a fresh slab and copies *elements*, and being range-addressed it handles
+  /// a strided (block-view) source correctly.
+  /// \tparam Right The source tensor type
+  /// \param right The tensor to copy
+  /// \return An owning `Tensor` sharing no storage with \p right
+  template <typename Right>
+  static Tensor owning_copy_of(const Right& right) {
+    using RightT = std::remove_cv_t<std::remove_reference_t<Right>>;
+    if constexpr (detail::is_tensor_of_tensor_v<Tensor> &&
+                  !std::is_same_v<RightT, Tensor> &&
+                  std::is_same_v<value_type, typename RightT::value_type>) {
+      auto fill = [](typename value_type::value_type* dst,
+                     const typename value_type::value_type* src,
+                     std::size_t n) {
+        for (std::size_t i = 0; i < n; ++i) dst[i] = src[i];
+      };
+      return detail::arena_trivial_unary<Tensor>(right, fill);
+    } else {
+      return detail::clone_or_cast<Tensor>(right);
+    }
+  }
+
   // Addition operations
 
   /// Element-wise add for `Tensor<ArenaTensor>` ToT operands. Routes through
@@ -2175,7 +2208,7 @@ class Tensor {
     requires(is_arena_tensor_v<value_type> &&
              is_arena_tensor_v<typename Right::value_type>)
   Tensor add(const Right& right) const {
-    if (empty()) return detail::clone_or_cast<Tensor>(right);
+    if (empty()) return owning_copy_of(right);
     if (right.empty()) return this->clone();
     auto fill = [](typename value_type::value_type* dst,
                    const typename value_type::value_type* l,
@@ -2290,7 +2323,7 @@ class Tensor {
     if (right.empty()) return this->clone();
 
     // early exit for empty this
-    if (empty()) return detail::clone_or_cast<Tensor>(right);
+    if (empty()) return owning_copy_of(right);
 
     if constexpr (detail::is_tensor_of_tensor_v<Tensor> &&
                   detail::is_ta_tensor_v<value_type> &&
@@ -2461,7 +2494,7 @@ class Tensor {
 
     // early exit for empty this
     if (empty()) {
-      *this = detail::clone_or_cast<Tensor>(right);
+      *this = owning_copy_of(right);
       return *this;
     }
 
@@ -2497,7 +2530,7 @@ class Tensor {
 
     // early exit for empty this: (0 + right) * factor
     if (empty()) {
-      *this = detail::clone_or_cast<Tensor>(right);
+      *this = owning_copy_of(right);
       this->scale_to(factor);
       return *this;
     }
@@ -2529,7 +2562,7 @@ class Tensor {
   Tensor& axpy_to(const Right& right, const Scalar factor) {
     if (right.empty()) return *this;
     if (empty()) {
-      *this = detail::clone_or_cast<Tensor>(right);
+      *this = owning_copy_of(right);
       this->scale_to(factor);
       return *this;
     }
@@ -2585,7 +2618,7 @@ class Tensor {
         // result inner cell): initialize to factor * (perm ^ arg) rather
         // than asserting non-empty in inplace_binary -- mirrors the
         // non-permuting axpy_to overload above.
-        *this = detail::clone_or_cast<Tensor>(permuted);
+        *this = owning_copy_of(permuted);
         this->scale_to(factor);
         return *this;
       }
@@ -2837,7 +2870,7 @@ class Tensor {
 
     // early exit for empty this: 0 - right == -right
     if (empty()) {
-      *this = detail::clone_or_cast<Tensor>(right);
+      *this = owning_copy_of(right);
       this->neg_to();
       return *this;
     }
@@ -2874,7 +2907,7 @@ class Tensor {
 
     // early exit for empty this: (0 - right) * factor
     if (empty()) {
-      *this = detail::clone_or_cast<Tensor>(right);
+      *this = owning_copy_of(right);
       this->neg_to();
       this->scale_to(factor);
       return *this;

tests/arena_tot_trivial.cpp

@@ -694,4 +694,85 @@ BOOST_AUTO_TEST_CASE(arena_mult_to_does_not_densify_null_cells) {
   }
 }
 
+// --- an empty-`this` fast path must not alias, let alone mutate, the source
+// `detail::clone_or_cast` deep-copies only when `Right` *is* `Tensor`; for any
+// other nested-tensor `Right` (e.g. the `TensorInterface<..., BlockRange>` that
+// `block()` yields) it copies the cell *handles*, which for a view cell aliases
+// the source's storage. Scaling or negating the copy then writes through to a
+// const operand.
+
+BOOST_AUTO_TEST_CASE(subt_to_empty_left_must_not_mutate_source) {
+  arena_outer_t P =
+      make_arena_tot_2d(2, 3, 4, 1.0, {true, true, true, true, true, true});
+  const double before = P.data()[0].data()[0];
+  arena_outer_t t;
+  t.subt_to(P.block({0l, 0l}, {2l, 2l}));
+  BOOST_CHECK_EQUAL(P.data()[0].data()[0], before);
+  BOOST_CHECK_EQUAL(t.data()[0].data()[0], -before);
+}
+
+BOOST_AUTO_TEST_CASE(subt_to_scaled_empty_left_must_not_mutate_source) {
+  arena_outer_t P =
+      make_arena_tot_2d(2, 3, 4, 1.0, {true, true, true, true, true, true});
+  const double before = P.data()[0].data()[0];
+  arena_outer_t t;
+  t.subt_to(P.block({0l, 0l}, {2l, 2l}), 2.0);
+  BOOST_CHECK_EQUAL(P.data()[0].data()[0], before);
+  BOOST_CHECK_EQUAL(t.data()[0].data()[0], -before * 2.0);
+}
+
+BOOST_AUTO_TEST_CASE(add_to_scaled_empty_left_must_not_mutate_source) {
+  arena_outer_t P =
+      make_arena_tot_2d(2, 3, 4, 1.0, {true, true, true, true, true, true});
+  const double before = P.data()[0].data()[0];
+  arena_outer_t t;
+  t.add_to(P.block({0l, 0l}, {2l, 2l}), 3.0);
+  BOOST_CHECK_EQUAL(P.data()[0].data()[0], before);
+  BOOST_CHECK_EQUAL(t.data()[0].data()[0], before * 3.0);
+}
+
+BOOST_AUTO_TEST_CASE(axpy_to_empty_left_must_not_mutate_source) {
+  arena_outer_t P =
+      make_arena_tot_2d(2, 3, 4, 1.0, {true, true, true, true, true, true});
+  const double before = P.data()[0].data()[0];
+  arena_outer_t t;
+  t.axpy_to(P.block({0l, 0l}, {2l, 2l}), 2.0);
+  BOOST_CHECK_EQUAL(P.data()[0].data()[0], before);
+}
+
+// the value-returning empty-left exits have the same hazard: the result must
+// own its cells, or the caller's later mutation corrupts the source
+BOOST_AUTO_TEST_CASE(add_empty_left_result_must_not_alias_source) {
+  arena_outer_t P =
+      make_arena_tot_2d(2, 3, 4, 1.0, {true, true, true, true, true, true});
+  const double before = P.data()[0].data()[0];
+  const arena_outer_t empty_left;
+  arena_outer_t r = empty_left.add(P.block({0l, 0l}, {2l, 2l}));
+  r.scale_to(-1.0);
+  BOOST_CHECK_EQUAL(P.data()[0].data()[0], before);
+}
+
+// --- intersection sparsity also governs the plain owning ToT --------------
+// `Tensor::mult`'s intersection rule is applied in the `is_ta_tensor_v` branch
+// too, so it is a behavior change for `TA::Tensor<TA::Tensor<T>>`, not only for
+// arena tiles. Pin it: the pre-existing mult_mismatched_null_inners is written
+// permissively and passes under either rule.
+BOOST_AUTO_TEST_CASE(plain_tot_mult_uses_intersection_sparsity) {
+  outer_t L = make_tot_sparse(5, 4, 2.0, nz_L);
+  outer_t R = make_tot_sparse(5, 4, 0.5, nz_R);
+  outer_t prod = L.mult(R);
+  for (std::size_t ord = 0; ord < 5; ++ord) {
+    const inner_t& l = *(L.data() + ord);
+    const inner_t& r = *(R.data() + ord);
+    const inner_t& d = *(prod.data() + ord);
+    if (!l.empty() && !r.empty()) {
+      BOOST_REQUIRE(!d.empty());
+      for (std::size_t i = 0; i < d.range().volume(); ++i)
+        BOOST_CHECK_EQUAL(d.at_ordinal(i), l.at_ordinal(i) * r.at_ordinal(i));
+    } else {
+      BOOST_CHECK(d.empty());
+    }
+  }
+}
+
 BOOST_AUTO_TEST_SUITE_END()