tensor: fix dropped rows/cols in the ToT x scalar strided scale paths

The GEMM-based ToT x scalar scale paths of Tensor::gemm (and the
T x ToT mirror) probe each row (column) for cleanliness; the presence
probe stops at the first ABSENT cell, leaving the probed inner size
A == -1 when the leading cell is absent. The subsequent 'A <= 0 =>
empty row, nothing to do' shortcut then dropped the ENTIRE row's
contributions even when later cells were present. Rows whose leading
contracted cell is absent (common for screened tensor-of-tensors, e.g.
PNO-CC CSV intermediates) silently lost their contraction.

Fix: when the probe ends with A <= 0, scan the full row (column) for
any present cell; only a fully absent row is skipped, anything else
takes the per-cell AXPY fallback. Also guard the engine's scale
fallback element op against absent cells.

This bug predates the general-product work but was masked on the
legacy einsum route, whose canonical operand layout fails the
NoTranspose gate of the strided path; the expression route's
GEMM-canonical layout exposed it. Found with TA_EINSUM_DIFFERENTIAL on
c6h14 PNO-CCSD: the opt-in expression-route energy error drops from
1.7 Eh to 6.9e-7 (the small residual is a screening-semantics
difference -- the legacy route hard-zeroes sub-threshold result tiles
that the engine route keeps -- plus a small systematic difference in
phantom-unit denest products, both under review).

Adds the CSV-like reproduction test (arena view cells, SparsePolicy,
variable inner extents, screened cells, non-leading fused index,
interleaved target): expression route and einsum routes agree to
1e-10, deterministically.

Author: evaleev

src/TiledArray/einsum/tiledarray.h

@@ -1096,6 +1096,51 @@ auto einsum(expressions::TsrExpr<ArrayA_> A, expressions::TsrExpr<ArrayB_> B,
                     << (std::string)b + inner.b << " -> " << c_annot
                     << " : diff2 = " << d2 << ", legacy2 = " << ref2
                     << std::endl;
+        // per-tile forensics: compare tile norms of the two routes
+        auto tile_norm2 = [](auto const &tile) -> double {
+          using TileT =
+              std::remove_cv_t<std::remove_reference_t<decltype(tile)>>;
+          double n2 = 0;
+          if constexpr (TiledArray::detail::is_tensor_of_tensor_v<TileT>) {
+            for (std::size_t o = 0; o < tile.range().volume(); ++o) {
+              auto const &cell = tile.data()[o];
+              if (cell.empty()) continue;
+              for (std::size_t e = 0; e < cell.range().volume(); ++e)
+                n2 += std::abs(cell.data()[e]) * std::abs(cell.data()[e]);
+            }
+          } else {
+            for (std::size_t e = 0; e < tile.range().volume() * tile.nbatch();
+                 ++e)
+              n2 += std::abs(tile.data()[e]) * std::abs(tile.data()[e]);
+          }
+          return n2;
+        };
+        const auto ntiles = C.array.trange().tiles_range().volume();
+        std::size_t n_diff = 0, n_expr_zero = 0, n_legacy_zero = 0,
+                    n_printed = 0;
+        for (std::size_t ord = 0; ord < ntiles; ++ord) {
+          const bool ez = expr_route_result->is_zero(ord);
+          const bool lz = C.array.is_zero(ord);
+          double en2 =
+              ez ? 0.0 : tile_norm2(expr_route_result->find(ord).get());
+          double ln2 = lz ? 0.0 : tile_norm2(C.array.find(ord).get());
+          const double dd = std::abs(en2 - ln2);
+          if (dd <= 1e-14 * std::max(std::max(en2, ln2), 1.0)) continue;
+          ++n_diff;
+          if (en2 == 0.0) ++n_expr_zero;
+          if (ln2 == 0.0) ++n_legacy_zero;
+          if (world.rank() == 0 && n_printed < 8) {
+            ++n_printed;
+            std::cerr << "    tile " << ord << " ("
+                      << C.array.trange().tiles_range().idx(ord)
+                      << "): expr_n2 = " << en2 << ", legacy_n2 = " << ln2
+                      << std::endl;
+          }
+        }
+        if (world.rank() == 0)
+          std::cerr << "    summary: " << n_diff << "/" << ntiles
+                    << " tiles differ by norm; expr-zero " << n_expr_zero
+                    << ", legacy-zero " << n_legacy_zero << std::endl;
       }
     }
 

src/TiledArray/expressions/cont_engine.h

@@ -1679,11 +1679,13 @@ class ContEngine : public BinaryEngine<Derived> {
             if (outer_uses_summa) {
               using TiledArray::axpy_to;
               if constexpr (tot_x_t) {
+                if (left.empty()) return;  // absent cell: no contribution
                 if (perm)
                   axpy_to(result, left, right, perm);
                 else
                   axpy_to(result, left, right);
               } else {
+                if (right.empty()) return;  // absent cell: no contribution
                 if (perm)
                   axpy_to(result, right, left, perm);
                 else

src/TiledArray/tensor/tensor.h

@@ -3374,7 +3374,19 @@ class Tensor {
               }
               detail::scale_phase_stop(detail::g_scale[0].check_str_ns,
                                        _scale_tcs);
-              if (A <= 0) continue;  // empty row -> nothing to do
+              if (A <= 0) {
+                // The presence probe stops at the first empty cell, so
+                // A <= 0 only means no non-empty cell was seen BEFORE the
+                // probe stopped -- NOT that the row is empty. The row
+                // contributes nothing only if EVERY left cell is absent;
+                // else it must take the per-cell AXPY fallback (whose
+                // element op skips absent cells).
+                bool row_empty = true;
+                for (integer k = 0; row_empty && k != K; ++k)
+                  if (!lc0[k].empty()) row_empty = false;
+                if (row_empty) continue;
+                clean = false;
+              }
               if (clean) {
                 // result[m,n][a] += sum_k left[m,k][a] * right[k,n].
                 // Row-major gemm: C2(N x A) += right^T(N x K) * L2(K x A),
@@ -3533,7 +3545,16 @@ class Tensor {
               }
               detail::scale_phase_stop(detail::g_scale[1].check_str_ns,
                                        _scale_tcs);
-              if (A <= 0) continue;
+              if (A <= 0) {
+                // see the ToT x scalar mirror above: A <= 0 does not imply
+                // an empty column when the probe stopped at the first
+                // absent cell
+                bool col_empty = true;
+                for (integer k = 0; col_empty && k != K; ++k)
+                  if (!right_data[k * N + n].empty()) col_empty = false;
+                if (col_empty) continue;
+                clean = false;
+              }
               if (clean) {
                 // C_n(M x A) += left(M x K) * B_n(K x A). Row-major gemm.
                 if (detail::scale_gemm_timing_enabled()) {

src/TiledArray/tile_op/batched_contract_reduce.h

@@ -193,8 +193,18 @@ class BatchedContractReduce {
           typename result_type::range_type(lobounds, upbounds));
     } else {
       // accumulate through a folded, zero-copy view of the result
-      auto result_folded = result.reshape(fold_range(result.range()), batch);
+      const auto full_range = result.range();
+      auto result_folded = result.reshape(fold_range(full_range), batch);
       op_(result_folded, left_folded, right_folded);
+      // the wrapped op may REBIND the result instead of writing in place:
+      // the arena grow-to-cover path (a later K-panel touching inner cells
+      // that earlier panels left null, for contracted-dimension-sparse
+      // tensor-of-tensor operands) replaces the tile wholesale
+      // (arena_tot_grow_inplace ends with `result = std::move(grown)`).
+      // That rebinding lands on the folded view; propagate it to the
+      // full-range result, else every post-growth contribution is lost.
+      if (result_folded.data() != result.data())
+        result = result_folded.reshape(full_range);
     }
   }