review: hoist the ld bound out of the segment loops; build the blas unit test

tests/math_blas.cpp is tracked but was never listed in ta_test_src_files, so
nothing in it has ever been compiled -- including the measured_ld_fits case
added for this guard. Registering it next to linalg.cpp brings the whole
blas_suite (17 cases) into the run; it passes as it stands.

ld_fits() divides, and the two ce+ce walkers called it twice per admitted
cell in the innermost segment-growth loop, whose other work is a handful of
loads and compares. The strides and extents are fixed for a segment once
off == 1, so take the bound once there instead: max_ld_offset(ld, extent) is
the same inequality solved for the slab offset, and admission then costs a
compare. The long arena-contiguous segments this kernel exists to hit --
g_seg_len9p_ce_ce -- were paying two divisions per cell for a check that can
never fire on them. ld_fits() is now expressed through max_ld_offset(), so
the two cannot drift apart; measure_segments takes the same hoist, since it
exists to simulate the walker's segmentation exactly.

Also, per review:

- spell math::blas::ld_fits at the call sites rather than pulling it into
  TiledArray::detail with a using-declaration: arena_einsum.h is widely
  included, and the name is of no interest to anything else in detail. The
  rationale it carried becomes a section comment.
- ce_e_stride_past_blas_int_falls_back returned before any assertion on both
  its skip paths (ILP64, refused reservation), so on those configurations it
  passed while checking nothing and only Boost's easily-missed "did not check
  any assertions" note said so. It now asserts that an ordinary contiguous
  2-cell run is still accepted before either branch, and the ILP64 skip
  explains itself like the reservation one already did.
- the reservation and the two placement-new'd cells get one RAII owner: a
  failing BOOST_REQUIRE skipped ~Cell(), leaking each cell's TA::Range
  buffer, which under this repo's ASan debug build would bury the real
  failure in leak reports.
- drop two lambda captures of constants (`len`, `K`) that -Wall reports as
  unnecessary. CI runs TA_WERROR=ON without -Wall so they did not gate it,
  but cmake/toolchains/travis.cmake does add -Wall, and that combination
  makes them errors.

Author: evaleev

src/TiledArray/math/blas.h

@@ -61,19 +61,30 @@ inline constexpr integer max_ld() {
   return static_cast<integer>(std::numeric_limits<blas_int>::max());
 }
 
+/// The largest row offset `o` for which element o*ld + extent-1 of a
+/// row-major matrix at leading dimension `ld` is still indexable by a
+/// `blas_int`, i.e. the largest `o` with ld_fits(ld, o+1, extent). Phrased as
+/// a division so the bound itself cannot overflow when `blas_int` is 64-bit
+/// (ILP64), where max_ld() is INT64_MAX. Callers that grow a run one row at a
+/// time should take this bound ONCE, while `ld` and `extent` are fixed,
+/// rather than calling ld_fits() per row -- that divides every time round the
+/// loop.
+inline constexpr integer max_ld_offset(integer ld, integer extent) {
+  if (ld <= 0 || extent < 1 || extent - 1 > max_ld()) return 0;
+  return (max_ld() - (extent - 1)) / ld;
+}
+
 /// True if a row-major matrix of `nslab` rows and `extent` columns at leading
 /// dimension `ld` addresses only elements a `blas_int` can index. max_ld()
 /// alone bounds ONE row step; the GEMM reaches element
 /// (nslab-1)*ld + extent-1, and a 32-bit-indexed BLAS must be able to form
 /// THAT index too, so a stride comfortably under the cap still overflows it
-/// once there are enough rows. Phrased as a division so the bound itself
-/// cannot overflow when `blas_int` is 64-bit (ILP64), where max_ld() is
-/// INT64_MAX.
+/// once there are enough rows.
 inline constexpr bool ld_fits(integer ld, integer nslab, integer extent) {
   if (ld < 0 || extent < 1 || extent - 1 > max_ld() || ld > max_ld())
     return false;
   if (nslab <= 1) return true;
-  return ld <= (max_ld() - (extent - 1)) / (nslab - 1);
+  return nslab - 1 <= max_ld_offset(ld, extent);
 }
 
 /// converts Op to ints in manner useful for bit manipulations

src/TiledArray/tensor/arena_einsum.h

@@ -294,20 +294,27 @@ inline std::atomic<std::uint64_t> g_fall_op_acrossk_ce_ce{0};      // 16
 inline std::atomic<std::uint64_t> g_fall_both_clean_ce_ce{0};      // 17
 inline std::atomic<std::uint64_t> g_fall_op_ldrange_ce_ce{0};      // 18
 
-/// Does a measured strided run address only elements the BLAS can index?
-/// The strided kernels below MEASURE their leading dimensions as pointer
-/// differences between neighbouring cells, so -- unlike a kernel that derives
-/// them from extents -- they can hand blaspp an ld it cannot convert to
-/// `blas_int` (`blas::Error: ldb, in function to_blas_int_`). Two present
-/// cells that are not arena neighbours -- e.g. individually allocated inner
-/// tensors the allocator placed gigabytes apart -- satisfy every other
-/// precondition of a 2-cell run (uniform size, stride >= cell size, trivially
-/// constant) yet yield exactly such an address delta; and a stride well inside
-/// `blas_int` still overflows the BLAS index once the run is long enough,
-/// because the GEMM reaches element (nslab-1)*ld + extent-1. Such a run is not
-/// strided-GEMM material: the walkers break it up and the per-cell path
-/// handles it.
-using math::blas::ld_fits;
+// ---------------------------------------------------------------------------
+// Measured leading dimensions vs. the BLAS index range.
+//
+// The strided kernels below MEASURE their leading dimensions as pointer
+// differences between neighbouring cells, so -- unlike a kernel that derives
+// them from extents -- they can hand blaspp an ld it cannot convert to
+// `blas_int` (`blas::Error: ldb, in function to_blas_int_`). Two present
+// cells that are not arena neighbours -- e.g. individually allocated inner
+// tensors the allocator placed gigabytes apart -- satisfy every other
+// precondition of a 2-cell run (uniform size, stride >= cell size, trivially
+// constant) yet yield exactly such an address delta; and a stride well inside
+// `blas_int` still overflows the BLAS index once the run is long enough,
+// because the GEMM reaches element (nslab-1)*ld + extent-1. Such a run is not
+// strided-GEMM material: the walkers break it up and the per-cell path
+// handles it.
+//
+// math::blas::ld_fits() is the whole-run test. The segment walkers instead
+// take math::blas::max_ld_offset() once per segment -- the same bound solved
+// for the slab offset -- because ld_fits() divides, and the admission loop
+// runs once per cell.
+// ---------------------------------------------------------------------------
 
 /// Classify a strided run: 1=absent, 2=nonuniform size, 3=bad stride,
 /// 4=stride past the BLAS index range (see ld_fits), 0=clean.
@@ -334,7 +341,7 @@ inline int classify_run(GetCell getcell, std::size_t n) {
     if (getcell(i).data() != base + static_cast<std::ptrdiff_t>(i) * st)
       return 3;  // non-constant stride
   // ...and the whole run must stay inside the BLAS index range
-  if (!ld_fits(st, static_cast<long>(n), s0)) return 4;
+  if (!math::blas::ld_fits(st, static_cast<long>(n), s0)) return 4;
   return 0;
 }
 
@@ -375,7 +382,7 @@ inline int classify_operand(GetR getR, GetL getL, std::size_t nrun,
       for (std::size_t i = 0; i < nrun; ++i)
         if (getR(k, i).data() != base + static_cast<std::ptrdiff_t>(i) * sk)
           return 15;
-      if (!ld_fits(sk, static_cast<long>(nrun), Q)) return 18;
+      if (!math::blas::ld_fits(sk, static_cast<long>(nrun), Q)) return 18;
       if (k == 0)
         sR = sk;
       else if (sk != sR)
@@ -513,6 +520,7 @@ inline void measure_segments(GetC getC, GetR getR, GetL getL, std::size_t nrun,
       const auto* rb = r0.data();
       std::size_t end = mu + 1;
       long sC = -1, sR = -1;
+      long max_off = 0;  // set with the strides, below
       while (end < nrun) {
         const auto& ce = getC(end);
         const auto& re = getR(k, end);
@@ -527,12 +535,15 @@ inline void measure_segments(GetC getC, GetR getR, GetL getL, std::size_t nrun,
           sC = dc;
           sR = dr;
           if (sC < P || sR < Q) break;
+          // A segment of off+1 slabs reaches element off*sC+P-1 of C and
+          // off*sR+Q-1 of R, both of which the BLAS has to be able to index.
+          // The strides and extents are fixed from here on, so bound off once.
+          max_off = std::min(math::blas::max_ld_offset(sC, P),
+                             math::blas::max_ld_offset(sR, Q));
         } else if (dc != off * sC || dr != off * sR) {
           break;
         }
-        // admitting this cell makes the segment off+1 slabs long; the GEMM
-        // would then reach element off*sC+P-1 of C and off*sR+Q-1 of R.
-        if (!ld_fits(sC, off + 1, P) || !ld_fits(sR, off + 1, Q)) break;
+        if (off > max_off) break;
         ++end;
       }
       const std::size_t len = end - mu;
@@ -1293,8 +1304,8 @@ void arena_strided_gemm_ce_e(ResultOuter& C, const LeftOuter& L,
           ldA = static_cast<long>(lc[lbase + a_off(m, 1)].data() - l0.data());
           ldB = static_cast<long>(rc[rbase + b_off(1, n)].data() - r0.data());
           if (ldA < P || ldB < Q) clean = false;
-          if (!ld_fits(ldA, static_cast<long>(K), P) ||
-              !ld_fits(ldB, static_cast<long>(K), Q))
+          if (!math::blas::ld_fits(ldA, static_cast<long>(K), P) ||
+              !math::blas::ld_fits(ldB, static_cast<long>(K), Q))
             clean = false;
           for (std::size_t k = 0; clean && k < K; ++k) {
             if (lc[lbase + a_off(m, k)].data() !=
@@ -1583,6 +1594,7 @@ void arena_strided_gemm_ce_ce_right(ResultOuter& C, const LeftOuter& L,
           // reuse a run-wide stale stride).
           std::size_t end = mu + 1;
           long sR = -1, sC = -1;
+          long max_off = 0;  // set with the strides, below
           while (end < Mmu) {
             const auto& rce = rc[rbase + r_off(k, end)];
             const auto& cce = cc[cbase + c_off(m, end)];
@@ -1597,13 +1609,16 @@ void arena_strided_gemm_ce_ce_right(ResultOuter& C, const LeftOuter& L,
               sR = dR;
               sC = dC;
               if (sR < Q || sC < P) break;  // page-jump / overlap
+              // A segment of off+1 slabs reaches element off*sR+Q-1 of A and
+              // off*sC+P-1 of C, both of which the BLAS has to be able to
+              // index. The strides and extents are fixed from here on, so
+              // bound off once rather than dividing per admitted cell.
+              max_off = std::min(blas::max_ld_offset(sR, Q),
+                                 blas::max_ld_offset(sC, P));
             } else if (dR != off * sR || dC != off * sC) {
               break;
             }
-            // admitting this cell makes the segment off+1 slabs long; the
-            // GEMM would then reach element off*sR+Q-1 of A and off*sC+P-1
-            // of C.
-            if (!ld_fits(sR, off + 1, Q) || !ld_fits(sC, off + 1, P)) break;
+            if (off > max_off) break;
             ++end;
           }
           const std::size_t Mseg = end - mu;
@@ -1830,6 +1845,7 @@ void arena_strided_gemm_ce_ce_left(ResultOuter& C, const LeftOuter& L,
           // reuse a run-wide stale stride).
           std::size_t end = m + 1;
           long sA = -1, sC = -1;
+          long max_off = 0;  // set with the strides, below
           while (end < Mo) {
             const auto& lce = lc[lbase + l_off(end, k)];
             const auto& cce = cc[cbase + c_off(end, n)];
@@ -1844,13 +1860,16 @@ void arena_strided_gemm_ce_ce_left(ResultOuter& C, const LeftOuter& L,
               sA = dA;
               sC = dC;
               if (sA < Q || sC < P) break;  // page-jump / overlap
+              // A segment of off+1 slabs reaches element off*sA+Q-1 of A and
+              // off*sC+P-1 of C, both of which the BLAS has to be able to
+              // index. The strides and extents are fixed from here on, so
+              // bound off once rather than dividing per admitted cell.
+              max_off = std::min(blas::max_ld_offset(sA, Q),
+                                 blas::max_ld_offset(sC, P));
             } else if (dA != off * sA || dC != off * sC) {
               break;
             }
-            // admitting this cell makes the segment off+1 slabs long; the
-            // GEMM would then reach element off*sA+Q-1 of A and off*sC+P-1
-            // of C.
-            if (!ld_fits(sA, off + 1, Q) || !ld_fits(sC, off + 1, P)) break;
+            if (off > max_off) break;
             ++end;
           }
           const std::size_t Mseg = end - m;

tests/CMakeLists.txt

@@ -97,6 +97,7 @@ set(ta_test_src_files  ta_test.cpp
     dot_inner.cpp
     general_product.cpp
     linalg.cpp
+    math_blas.cpp
     cp.cpp
     btas.cpp
     arena.cpp

tests/arena_strided_gemm.cpp

@@ -391,8 +391,17 @@ BOOST_AUTO_TEST_CASE(ce_e_stride_past_blas_int_falls_back) {
   // second cell sits in lazily reserved address space, so only the two cells'
   // pages are ever touched; if the reservation is refused the case is skipped.
   namespace blas = TiledArray::math::blas;
-  if (blas::max_ld() >= std::numeric_limits<blas::integer>::max()) return;
   const std::size_t M = 1, N = 1, K = 2, P = 3, Q = 4;
+  // Always run, so neither skip below leaves the case without a single
+  // assertion: an ordinary contiguous 2-cell run is still accepted.
+  BOOST_CHECK(blas::ld_fits(static_cast<blas::integer>(Q), 2,
+                            static_cast<blas::integer>(Q)));
+  if (blas::max_ld() >= std::numeric_limits<blas::integer>::max()) {
+    BOOST_TEST_MESSAGE(
+        "ILP64: blas_int is as wide as integer, so no stride can outrun it; "
+        "skipped");
+    return;
+  }
   const std::size_t stride = static_cast<std::size_t>(blas::max_ld()) + 1;
   const std::size_t align = Inner::cell_alignment();
   const std::size_t len = stride * sizeof(double) + Inner::cell_size(P) + align;
@@ -403,14 +412,28 @@ BOOST_AUTO_TEST_CASE(ce_e_stride_past_blas_int_falls_back) {
                                          << " bytes of address space; skipped");
     return;
   }
-  const auto unmap = [len](void* p) { ::munmap(p, len); };
-  std::unique_ptr<void, decltype(unmap)> mapping(mem, unmap);
+  using Cell = Inner::Cell;
+  // Owns the reservation and the placement-new'd cells together: a failing
+  // BOOST_REQUIRE below must not skip ~Cell(), whose TA::Range owns a heap
+  // buffer -- under this repo's ASan debug build that would bury the real
+  // failure in leak reports.
+  struct Reservation {
+    void* mem;
+    std::size_t len;
+    Cell* c0 = nullptr;
+    Cell* c1 = nullptr;
+    ~Reservation() {
+      if (c1) c1->~Cell();
+      if (c0) c0->~Cell();
+      ::munmap(mem, len);
+    }
+  } reservation{mem, len};
   auto* base = reinterpret_cast<std::byte*>(
       (reinterpret_cast<std::uintptr_t>(mem) + align - 1) & ~(align - 1));
   // stride * sizeof(double) is a power of two >= align, so cell 1 is aligned
-  using Cell = Inner::Cell;
-  Cell* c0 = ::new (base) Cell{TA::Range{P}};
-  Cell* c1 = ::new (base + stride * sizeof(double)) Cell{TA::Range{P}};
+  Cell* c0 = reservation.c0 = ::new (base) Cell{TA::Range{P}};
+  Cell* c1 = reservation.c1 =
+      ::new (base + stride * sizeof(double)) Cell{TA::Range{P}};
   {
     Inner l0(c0), l1(c1);
     BOOST_REQUIRE_EQUAL(static_cast<std::size_t>(l1.data() - l0.data()),
@@ -423,9 +446,9 @@ BOOST_AUTO_TEST_CASE(ce_e_stride_past_blas_int_falls_back) {
     Inner* raw = alloc.allocate(K);
     ::new (raw) Inner(l0);
     ::new (raw + 1) Inner(l1);
-    std::shared_ptr<Inner[]> ldata(raw, [alloc, K](Inner* p) mutable {
-      for (std::size_t i = 0; i < K; ++i) (p + i)->~Inner();
-      alloc.deallocate(p, K);
+    std::shared_ptr<Inner[]> ldata(raw, [alloc, n = K](Inner* p) mutable {
+      for (std::size_t i = 0; i < n; ++i) (p + i)->~Inner();
+      alloc.deallocate(p, n);
     });
     Outer L(TA::Range{M, K}, /*nbatch=*/1, std::move(ldata));
     Outer R = make_filled(
@@ -438,9 +461,7 @@ BOOST_AUTO_TEST_CASE(ce_e_stride_past_blas_int_falls_back) {
     const double* got = C.data()[0].data();
     for (std::size_t e = 0; e < P * Q; ++e)
       BOOST_CHECK_CLOSE(got[e], ref[e], 1e-12);
-  }  // the views die before their cells
-  c1->~Cell();
-  c0->~Cell();
+  }  // the views die before their cells, which ~Reservation then destroys
 }
 #ifdef TA_STRIDED_GEMM_COUNT
 BOOST_AUTO_TEST_CASE(ce_e_fires_clean_path) {