From 4f1336eb0427aafb6dc7d9a37e4f6fadf9f53767 Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Wed, 9 Sep 2026 13:32:19 -0400 Subject: [PATCH 1/3] Tensor::gemm: arena scale fast paths serve both operand orientations The two strided-GEMM fast paths for the "scale" contraction of a nested (tensor-of-tensors) tile with a plain tile -- nested-left, one GEMM per result row, and its nested-right mirror, one GEMM per result column -- were gated on NoTranspose/NoTranspose. The GEMM permutation optimizer hands a nested-right operand whose contracted mode trails to the GEMM as an implicit transpose (right_is_trans), so that operand order dropped to the per-cell AXPY loop: the same product measured ~2.5x slower flat-first than nested-first on two tiles per mode and ~4.5x on one tile, for any result order and with or without an imposed result shape. Both paths now take either orientation: the plain matrix's transpose is passed to BLAS as its op (with the matching leading dimension), and a transposed nested tile only changes which cells form a row (cells (k,m) stored k-major, M cells apart) or a column (cells (n,k) contiguous); the stride is measured from the cell addresses as before. The per-cell fallbacks index both operands by the same rule as the generic loop. tests/arena_mixed_product_order.cpp: a DistArray-level regression that asserts all spellings of the mixed product agree (both operand orders, plain and shaped, both result orders) and reports their wall times; its second case documents a separate, still-open defect (a permuted read of an arena nested array inside a binary expression corrupts the source). --- src/TiledArray/tensor/tensor.h | 94 +++++++++++------ tests/CMakeLists.txt | 1 + tests/arena_mixed_product_order.cpp | 155 ++++++++++++++++++++++++++++ 3 files changed, 219 insertions(+), 31 deletions(-) create mode 100644 tests/arena_mixed_product_order.cpp diff --git a/src/TiledArray/tensor/tensor.h b/src/TiledArray/tensor/tensor.h index bc7a4de9d3..f8abcf4856 100644 --- a/src/TiledArray/tensor/tensor.h +++ b/src/TiledArray/tensor/tensor.h @@ -3568,8 +3568,19 @@ class Tensor { is_tensor_view_v) { using Real = std::remove_cv_t; if constexpr (std::is_same_v, Real>) { - if (gemm_helper.left_op() == TiledArray::math::blas::NoTranspose && - gemm_helper.right_op() == TiledArray::math::blas::NoTranspose) { + { + // 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 + // stored k-major, so the row's cells are M cells apart; the stride + // is measured from the cell addresses either way). + const bool left_no_trans = + gemm_helper.left_op() == TiledArray::math::blas::NoTranspose; + const bool right_no_trans = + gemm_helper.right_op() == TiledArray::math::blas::NoTranspose; + auto lcell = [&](integer m, integer k) -> integer { + return left_no_trans ? m * lda + k : k * lda + m; + }; // kernel-total timer: destroyed at `return *this;` below, so it // captures the whole for-b/for-m loop. loop-residual is derived from // it minus the sub-phases. @@ -3583,7 +3594,9 @@ class Tensor { auto left_data = left.batch_data(b); auto right_data = right.batch_data(b); // K x N row-major scalars for (integer m = 0; m != M; ++m) { - auto* lc0 = left_data + (m * K); // left cells (m,0..K-1) + auto lc = [&](integer k) -> decltype(auto) { + return left_data[lcell(m, k)]; // left cell (m,k) + }; auto* rc0 = this_data + (m * N); // result cells (m,0..N-1) // A "clean" row has all cells present, uniform inner size A, and // laid out as one contiguous stride-A block (so the GEMM can run @@ -3593,7 +3606,7 @@ class Tensor { bool clean = true; const auto _scale_tcp = detail::scale_phase_start(); for (integer k = 0; k != K && clean; ++k) { - const auto& c = lc0[k]; + const auto& c = lc(k); if (c.empty()) { clean = false; break; @@ -3630,13 +3643,13 @@ class Tensor { integer ldc = static_cast(A); if (clean && A > 0) { if (K > 1) - ldb = static_cast(lc0[1].data() - lc0[0].data()); + ldb = static_cast(lc(1).data() - lc(0).data()); if (N > 1) ldc = static_cast(rc0[1].data() - rc0[0].data()); if (ldb < A || ldc < A) clean = false; // sanity const std::ptrdiff_t sb = ldb, sc = ldc; for (integer k = 0; clean && k != K; ++k) - if (lc0[k].data() != lc0[0].data() + k * sb) clean = false; + if (lc(k).data() != lc(0).data() + k * sb) clean = false; for (integer n = 0; clean && n != N; ++n) if (rc0[n].data() != rc0[0].data() + n * sc) clean = false; } @@ -3651,7 +3664,7 @@ class Tensor { // 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 (!lc(k).empty()) row_empty = false; if (row_empty) continue; clean = false; } @@ -3672,12 +3685,15 @@ class Tensor { } const integer Ai = static_cast(A); detail::ScopedScaleTimer _scale_gt(detail::g_scale[0].gemm_ns); + // right is K x N (ld N) when NoTranspose -> pass it + // transposed; stored N x K (ld K) when Transpose -> as is. TiledArray::math::blas::gemm( - TiledArray::math::blas::Transpose, + right_no_trans ? TiledArray::math::blas::Transpose + : TiledArray::math::blas::NoTranspose, TiledArray::math::blas::NoTranspose, /*M=*/N, /*N=*/Ai, /*K=*/K, Real(1), - /*A=*/right_data, /*lda=*/N, - /*B=*/lc0[0].data(), /*ldb=*/ldb, Real(1), + /*A=*/right_data, /*lda=*/right_no_trans ? N : K, + /*B=*/lc(0).data(), /*ldb=*/ldb, Real(1), /*C=*/rc0[0].data(), /*ldc=*/ldc); } else { // per-cell AXPY fallback for this row if (detail::scale_gemm_timing_enabled()) { @@ -3686,7 +3702,7 @@ class Tensor { bool absent = false, ragged = false; long a0 = -1; for (integer k = 0; k != K; ++k) { - const auto& c = lc0[k]; + const auto& c = lc(k); if (c.empty()) { absent = true; break; @@ -3727,9 +3743,11 @@ class Tensor { for (integer n = 0; n != N; ++n) { auto c_offset = m * N + n; for (integer k = 0; k != K; ++k) - elem_muladd_op(*(this_data + c_offset), - *(left_data + (m * K + k)), - *(right_data + (k * N + n))); + // N.B. `ldb` is shadowed here by the cell stride; index + // the plain right matrix by its own leading dimension. + elem_muladd_op(*(this_data + c_offset), lc(k), + *(right_data + + (right_no_trans ? k * N + n : n * K + k))); } } } @@ -3748,8 +3766,18 @@ class Tensor { is_tensor_view_v) { using Real = std::remove_cv_t; if constexpr (std::is_same_v, Real>) { - if (gemm_helper.left_op() == TiledArray::math::blas::NoTranspose && - gemm_helper.right_op() == TiledArray::math::blas::NoTranspose) { + { + // 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 + // then one contiguous run (the friendlier layout). + const bool left_no_trans = + gemm_helper.left_op() == TiledArray::math::blas::NoTranspose; + const bool right_no_trans = + gemm_helper.right_op() == TiledArray::math::blas::NoTranspose; + auto rcell = [&](integer k, integer n) -> integer { + return right_no_trans ? k * ldb + n : n * ldb + k; + }; // kernel-total timer (see tot_x_t block); destroyed at `return`. detail::ScopedScaleTimer _scale_kt(detail::g_scale[1].kernel_ns); if (detail::scale_gemm_timing_enabled()) @@ -3765,7 +3793,7 @@ class Tensor { bool clean = true; const auto _scale_tcp = detail::scale_phase_start(); for (integer k = 0; k != K && clean; ++k) { - const auto& c = right_data[k * N + n]; + const auto& c = right_data[rcell(k, n)]; if (c.empty()) { clean = false; break; @@ -3791,20 +3819,20 @@ class Tensor { detail::scale_phase_stop(detail::g_scale[1].check_pres_ns, _scale_tcp); const auto _scale_tcs = detail::scale_phase_start(); - integer ldb = static_cast(A); // k-stride, right col n + integer sbc = static_cast(A); // k-stride, right col n integer ldc = static_cast(A); // m-stride, result col n if (clean && A > 0) { if (K > 1) - ldb = static_cast(right_data[N + n].data() - - right_data[n].data()); + sbc = static_cast(right_data[rcell(1, n)].data() - + right_data[rcell(0, n)].data()); if (M > 1) ldc = static_cast(this_data[N + n].data() - this_data[n].data()); - if (ldb < A || ldc < A) clean = false; - const std::ptrdiff_t sb = ldb, sc = ldc; + if (sbc < A || ldc < A) clean = false; + const std::ptrdiff_t sb = sbc, sc = ldc; for (integer k = 0; clean && k != K; ++k) - if (right_data[k * N + n].data() != - right_data[n].data() + k * sb) + if (right_data[rcell(k, n)].data() != + right_data[rcell(0, n)].data() + k * sb) clean = false; for (integer m = 0; clean && m != M; ++m) if (this_data[m * N + n].data() != @@ -3819,7 +3847,7 @@ class Tensor { // 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 (!right_data[rcell(k, n)].empty()) col_empty = false; if (col_empty) continue; clean = false; } @@ -3836,12 +3864,15 @@ class Tensor { } const integer Ai = static_cast(A); detail::ScopedScaleTimer _scale_gt(detail::g_scale[1].gemm_ns); + // left is M x K (ld K) when NoTranspose; stored K x M (ld M) + // when Transpose -> pass BLAS the op. TiledArray::math::blas::gemm( - TiledArray::math::blas::NoTranspose, + left_no_trans ? TiledArray::math::blas::NoTranspose + : TiledArray::math::blas::Transpose, TiledArray::math::blas::NoTranspose, /*M=*/M, /*N=*/Ai, /*K=*/K, Real(1), - /*A=*/left_data, /*lda=*/K, - /*B=*/right_data[n].data(), /*ldb=*/ldb, Real(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); } else { // per-cell AXPY fallback for this column if (detail::scale_gemm_timing_enabled()) { @@ -3850,7 +3881,7 @@ class Tensor { bool absent = false, ragged = false; long a0 = -1; for (integer k = 0; k != K; ++k) { - const auto& c = right_data[k * N + n]; + const auto& c = right_data[rcell(k, n)]; if (c.empty()) { absent = true; break; @@ -3893,8 +3924,9 @@ class Tensor { auto c_offset = m * N + n; for (integer k = 0; k != K; ++k) elem_muladd_op(*(this_data + c_offset), - *(left_data + (m * K + k)), - *(right_data + (k * N + n))); + *(left_data + (left_no_trans ? m * lda + k + : k * lda + m)), + *(right_data + rcell(k, n))); } } } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5b71ee6c53..e0da3296b7 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -107,6 +107,7 @@ set(ta_test_src_files ta_test.cpp arena_tensor.cpp arena_tensor_kernels.cpp arena_strided_gemm.cpp + arena_mixed_product_order.cpp tot_construction.cpp ) diff --git a/tests/arena_mixed_product_order.cpp b/tests/arena_mixed_product_order.cpp new file mode 100644 index 0000000000..36c187ea29 --- /dev/null +++ b/tests/arena_mixed_product_order.cpp @@ -0,0 +1,155 @@ +// tests/arena_mixed_product_order.cpp +// +// Two DistArray-level regressions for MIXED products (an arena tensor-of- +// tensors operand times a plain tensor operand), found on the CSV-CC +// per-batch intermediate (occ,occ,PAO;PNO) * (PAO,PAO,DF) -> +// (occ,occ,PAO,DF;PNO): +// +// 1. operand order: the expression layer's cost of a mixed product depended +// on which operand came first (nested-first ~2.5x faster on two tiles +// per mode, ~4.5x on one tile, regardless of the result order or an +// imposed shape). Cause: Tensor::gemm's two arena "scale" strided-GEMM +// fast paths (nested-left rows, nested-right columns) accepted only the +// NoTranspose/NoTranspose orientation, and the permutation optimizer +// hands a nested-right operand whose contracted mode trails to the GEMM +// as an implicit transpose -- so that order fell to the per-cell AXPY +// loop. Both paths now serve either orientation (the plain matrix's +// transpose goes to BLAS as its op; a transposed nested tile only changes +// which cells form a row/column). The test asserts that all spellings +// (both orders, plain and shaped, both result orders) agree and reports +// the wall time of each (informational; no timing assertion). +// +// 2. permuted operand in a binary op: reading an arena ToT through a +// PERMUTED annotation inside a binary expression (e.g. `x(perm) - y`) +// destroyed the source: the permuted operand tile is a SHALLOW permute +// (cells alias the source slab) and the in-place binary op wrote through +// it. A permute-COPY (`y(perm) = x(...)`) was unaffected, and owning +// (Tensor) inners were unaffected. The test asserts the source's +// norm is unchanged after such a read. + +#include "TiledArray/tensor/arena_einsum.h" +#include "TiledArray/tensor/arena_tensor.h" +#include "tiledarray.h" +#include "unit_test_config.h" + +#include +#include +#include + +namespace TA = TiledArray; + +namespace { + +using ArenaInner = TA::ArenaTensor; +using ArenaOuter = TA::Tensor; +using ToTArray = TA::DistArray; +using FlatArray = TA::DistArray, TA::SparsePolicy>; + +struct MixedFixture { + long const I = 8, P = 24, K = 10, A = 6; + TA::TiledRange1 const tr_i{0l, I / 2, I}, tr_p{0l, P / 2, P}, tr_k{0l, K}; + TA::TiledRange const ij_p{tr_i, tr_i, tr_p}; + TA::TiledRange const q_p_k{tr_p, tr_p, tr_k}; + TA::TiledRange const ijqk{tr_i, tr_i, tr_p, tr_k}; + std::string const ta = "i,j,p;a", tb = "q,p,k", tc = "i,j,q,k;a", + tc_alt = "q,k,i,j;a"; + ToTArray tot; + FlatArray flat; + + MixedFixture() { + auto& world = TA::get_default_world(); + tot = ToTArray(world, ij_p); + tot.init_tiles([&](TA::Range const& tr) { + ArenaOuter t = TA::detail::arena_outer_init( + tr, 1, [=](std::size_t) { return TA::Range{A}; }); + for (std::size_t o = 0; o < t.range().volume(); ++o) { + ArenaInner& c = t.data()[o]; + if (!c) continue; + for (long a = 0; a < A; ++a) + c.data()[a] = 1.0 + 0.001 * static_cast((o * 7 + a) % 13); + } + return t; + }); + flat = FlatArray(world, q_p_k); + flat.init_tiles([&](TA::Range const& tr) { + TA::Tensor t(tr); + for (std::size_t o = 0; o < t.range().volume(); ++o) + t.data()[o] = 0.5 + 0.001 * static_cast(o % 17); + return t; + }); + world.gop.fence(); + } + + template + static double timed(F&& f) { + auto& world = TA::get_default_world(); + world.gop.fence(); + auto const t0 = std::chrono::steady_clock::now(); + f(); + world.gop.fence(); + return std::chrono::duration(std::chrono::steady_clock::now() - t0) + .count(); + } + + double diff(ToTArray const& x, ToTArray const& ref, std::string const& xa) { + ToTArray d; + d(tc) = x(xa) - ref(tc); + return TA::norm2(d); + } +}; + +} // namespace + +BOOST_FIXTURE_TEST_SUITE(arena_mixed_product_order_suite, MixedFixture, + TA_UT_LABEL_SERIAL) + +BOOST_AUTO_TEST_CASE(operand_order_is_canonicalized) { + ToTArray ref, ff, nf_sh, ff_sh, nf_alt, ff_alt; + double const t_nf = timed([&] { ref(tc) = tot(ta) * flat(tb); }); + double const t_ff = timed([&] { ff(tc) = flat(tb) * tot(ta); }); + TA::Tensor norms(ijqk.tiles_range(), 1.0f); + TA::SparseShape const shape(norms, ijqk, /*do_not_scale=*/true); + double const t_nf_sh = + timed([&] { nf_sh(tc) = (tot(ta) * flat(tb)).set_shape(shape); }); + double const t_ff_sh = + timed([&] { ff_sh(tc) = (flat(tb) * tot(ta)).set_shape(shape); }); + double const t_nf_alt = timed([&] { nf_alt(tc_alt) = tot(ta) * flat(tb); }); + double const t_ff_alt = timed([&] { ff_alt(tc_alt) = flat(tb) * tot(ta); }); + double const n = TA::norm2(ref); + BOOST_REQUIRE(n > 0.0); + double const tol = 1e-12 * n; + BOOST_CHECK_SMALL(diff(ff, ref, tc), tol); + BOOST_CHECK_SMALL(diff(nf_sh, ref, tc), tol); + BOOST_CHECK_SMALL(diff(ff_sh, ref, tc), tol); + BOOST_CHECK_SMALL(diff(nf_alt, ref, tc_alt), tol); + BOOST_CHECK_SMALL(diff(ff_alt, ref, tc_alt), tol); + BOOST_TEST_MESSAGE("mixed product wall time (s): nested*flat=" + << t_nf << " flat*nested=" << t_ff + << " +shape: " << t_nf_sh << " / " << t_ff_sh + << " result q,k,i,j: " << t_nf_alt << " / " << t_ff_alt); +} + +BOOST_AUTO_TEST_CASE(permuted_arena_operand_in_binary_op_keeps_source) { + ToTArray ref, src_copy, src_binary; + ref(tc) = tot(ta) * flat(tb); + src_copy(tc) = tot(ta) * flat(tb); + src_binary(tc) = tot(ta) * flat(tb); + TA::get_default_world().gop.fence(); + double const n0 = TA::norm2(src_copy); + BOOST_REQUIRE(n0 > 0.0); + + ToTArray y; + y(tc_alt) = src_copy(tc); // permute-copy: must not touch the source + TA::get_default_world().gop.fence(); + BOOST_CHECK_CLOSE(TA::norm2(src_copy), n0, 1e-10); + + ToTArray z; + z(tc) = src_binary(tc_alt) - ref(tc); // permuted operand of a binary op + TA::get_default_world().gop.fence(); + BOOST_CHECK_CLOSE(TA::norm2(src_binary), n0, 1e-10); + // and the binary op itself was right: |z| is |ref_perm - ref| = 0 only if + // the permutation was applied; here they differ, so just require finite. + BOOST_CHECK(std::isfinite(TA::norm2(z))); +} + +BOOST_AUTO_TEST_SUITE_END() From 282a8765ad8b53280ed7cead46c2edd78022bb21 Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Wed, 9 Sep 2026 14:16:20 -0400 Subject: [PATCH 2/3] tests/arena_mixed_product_order: capture the inner extent explicitly (implicit this capture is an error under the CI warning policy) --- tests/arena_mixed_product_order.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/arena_mixed_product_order.cpp b/tests/arena_mixed_product_order.cpp index 36c187ea29..df374aecab 100644 --- a/tests/arena_mixed_product_order.cpp +++ b/tests/arena_mixed_product_order.cpp @@ -61,7 +61,7 @@ struct MixedFixture { tot = ToTArray(world, ij_p); tot.init_tiles([&](TA::Range const& tr) { ArenaOuter t = TA::detail::arena_outer_init( - tr, 1, [=](std::size_t) { return TA::Range{A}; }); + tr, 1, [A = this->A](std::size_t) { return TA::Range{A}; }); for (std::size_t o = 0; o < t.range().volume(); ++o) { ArenaInner& c = t.data()[o]; if (!c) continue; From cafebd6cd7ad0329b0a18689e6799264e68d46c8 Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Wed, 9 Sep 2026 16:27:05 -0400 Subject: [PATCH 3/3] review: ConjTranspose takes the generic loop; test includes ; permuted-operand case moves to its fix Neither scale fast path conjugates (nor did the per-cell op they replace), so a ConjTranspose operand op now bails out to the generic loop instead of being treated as a plain transpose. tests/arena_mixed_product_order.cpp: explicit for std::isfinite; the second case (a permuted arena operand consumed inside a binary expression) is dropped here -- as written it also relabeled modes of different extents, which a Debug build rejects at the binary engine's TiledRange check -- and lands, corrected, with the fix for that defect. --- src/TiledArray/tensor/tensor.h | 27 +++++++++++++--------- tests/arena_mixed_product_order.cpp | 35 ++++------------------------- 2 files changed, 20 insertions(+), 42 deletions(-) diff --git a/src/TiledArray/tensor/tensor.h b/src/TiledArray/tensor/tensor.h index f8abcf4856..895e97864c 100644 --- a/src/TiledArray/tensor/tensor.h +++ b/src/TiledArray/tensor/tensor.h @@ -3568,12 +3568,15 @@ class Tensor { is_tensor_view_v) { using Real = std::remove_cv_t; if constexpr (std::is_same_v, Real>) { - { - // 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 - // stored k-major, so the row's cells are M cells apart; the stride - // is measured from the cell addresses either way). + // 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 + // stored k-major, so the row's cells are M cells apart; the stride + // is measured from the cell addresses either way). A conjugating op + // is not: neither this path nor the per-cell op below conjugates, + // so ConjTranspose takes the generic loop. + if (gemm_helper.left_op() != TiledArray::math::blas::ConjTranspose && + gemm_helper.right_op() != TiledArray::math::blas::ConjTranspose) { const bool left_no_trans = gemm_helper.left_op() == TiledArray::math::blas::NoTranspose; const bool right_no_trans = @@ -3766,11 +3769,13 @@ class Tensor { is_tensor_view_v) { using Real = std::remove_cv_t; if constexpr (std::is_same_v, Real>) { - { - // 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 - // then one contiguous run (the friendlier layout). + // 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 + // then one contiguous run (the friendlier layout). ConjTranspose + // takes the generic loop, as above. + if (gemm_helper.left_op() != TiledArray::math::blas::ConjTranspose && + gemm_helper.right_op() != TiledArray::math::blas::ConjTranspose) { const bool left_no_trans = gemm_helper.left_op() == TiledArray::math::blas::NoTranspose; const bool right_no_trans = diff --git a/tests/arena_mixed_product_order.cpp b/tests/arena_mixed_product_order.cpp index df374aecab..39663a7e85 100644 --- a/tests/arena_mixed_product_order.cpp +++ b/tests/arena_mixed_product_order.cpp @@ -1,6 +1,6 @@ // tests/arena_mixed_product_order.cpp // -// Two DistArray-level regressions for MIXED products (an arena tensor-of- +// DistArray-level regression for MIXED products (an arena tensor-of- // tensors operand times a plain tensor operand), found on the CSV-CC // per-batch intermediate (occ,occ,PAO;PNO) * (PAO,PAO,DF) -> // (occ,occ,PAO,DF;PNO): @@ -19,13 +19,8 @@ // (both orders, plain and shaped, both result orders) agree and reports // the wall time of each (informational; no timing assertion). // -// 2. permuted operand in a binary op: reading an arena ToT through a -// PERMUTED annotation inside a binary expression (e.g. `x(perm) - y`) -// destroyed the source: the permuted operand tile is a SHALLOW permute -// (cells alias the source slab) and the in-place binary op wrote through -// it. A permute-COPY (`y(perm) = x(...)`) was unaffected, and owning -// (Tensor) inners were unaffected. The test asserts the source's -// norm is unchanged after such a read. +// (A second case, a permuted arena operand consumed inside a binary +// expression, lands with the fix for that defect.) #include "TiledArray/tensor/arena_einsum.h" #include "TiledArray/tensor/arena_tensor.h" @@ -33,6 +28,7 @@ #include "unit_test_config.h" #include +#include #include #include @@ -129,27 +125,4 @@ BOOST_AUTO_TEST_CASE(operand_order_is_canonicalized) { << " result q,k,i,j: " << t_nf_alt << " / " << t_ff_alt); } -BOOST_AUTO_TEST_CASE(permuted_arena_operand_in_binary_op_keeps_source) { - ToTArray ref, src_copy, src_binary; - ref(tc) = tot(ta) * flat(tb); - src_copy(tc) = tot(ta) * flat(tb); - src_binary(tc) = tot(ta) * flat(tb); - TA::get_default_world().gop.fence(); - double const n0 = TA::norm2(src_copy); - BOOST_REQUIRE(n0 > 0.0); - - ToTArray y; - y(tc_alt) = src_copy(tc); // permute-copy: must not touch the source - TA::get_default_world().gop.fence(); - BOOST_CHECK_CLOSE(TA::norm2(src_copy), n0, 1e-10); - - ToTArray z; - z(tc) = src_binary(tc_alt) - ref(tc); // permuted operand of a binary op - TA::get_default_world().gop.fence(); - BOOST_CHECK_CLOSE(TA::norm2(src_binary), n0, 1e-10); - // and the binary op itself was right: |z| is |ref_perm - ref| = 0 only if - // the permutation was applied; here they differ, so just require finite. - BOOST_CHECK(std::isfinite(TA::norm2(z))); -} - BOOST_AUTO_TEST_SUITE_END()