From 801dd2655f583724f59b3481f87c0f1b26864234 Mon Sep 17 00:00:00 2001 From: xyhe <1225385871@qq.com> Date: Mon, 7 Sep 2026 20:39:26 +0800 Subject: [PATCH] perf(builder): replace hash-map dedup with sort-based CSR assembly cCSRMatrixBuilder::Build() previously deduplicated triplets through a global std::unordered_map: one heap node per unique entry (~48-64B), a hashed lookup per triplet, a reserve() sized on the full triplet count, and a final global O(n log n) sort. For near-field assemblies with tens of millions of triplets this dominated ILU-preconditioner build time (issue #7: 210s at the 3-lambda case) and peak memory (~4x the final CSR size). Replace it with a streaming four-pass pipeline: 1. count in-range entries per row -> rowPtr prefix sum O(N) 2. scatter (col, value) into row-segmented slots O(N) 3. per-row std::stable_sort by column + merge of equal- column runs (accumulate or keep-last) 4. final row offsets while emitting the CSR arrays std::stable_sort preserves insertion order within equal columns, so duplicates accumulate in the same order - and to the same bits - as the hash-map implementation. The public API is unchanged. Measured (MinGW g++ -O2, near-field pattern, 10% duplicates): 7.7M triplets: 1.06s -> 0.106s (~10x) 31M triplets: 0.439s (~71M triplets/s) Tests: all 17 existing CSRBuilderTest cases pass unchanged; add five bit-exact equivalence tests against a verbatim copy of the old hash-map pipeline (accumulate/keep-last, near-field pattern with out-of-range rows, edge cases, complex scalars). Add bench/bench_builder.cpp reproducing the near-field assembly. Fixes #7 --- bench/CMakeLists.txt | 6 + bench/bench_builder.cpp | 85 ++++++++++++++ include/CSRMatrixBuilder.h | 142 ++++++++++++---------- tests/test_csr_builder.cpp | 235 ++++++++++++++++++++++++++++++++++++- 4 files changed, 400 insertions(+), 68 deletions(-) create mode 100644 bench/bench_builder.cpp diff --git a/bench/CMakeLists.txt b/bench/CMakeLists.txt index 343587e..7253b1d 100644 --- a/bench/CMakeLists.txt +++ b/bench/CMakeLists.txt @@ -9,4 +9,10 @@ add_executable(csr4mpi_bench_spmv ) target_link_libraries(csr4mpi_bench_spmv PRIVATE csr4mpi) target_compile_definitions(csr4mpi_bench_spmv PRIVATE CSR4MPI_SOURCE_DIR="${CMAKE_SOURCE_DIR}" CSR4MPI_VALUE_TYPE=${CSR4MPI_VALUE_TYPE}) + +add_executable(csr4mpi_bench_builder + bench_builder.cpp +) +target_link_libraries(csr4mpi_bench_builder PRIVATE csr4mpi) +target_compile_definitions(csr4mpi_bench_builder PRIVATE CSR4MPI_SOURCE_DIR="${CMAKE_SOURCE_DIR}" CSR4MPI_VALUE_TYPE=${CSR4MPI_VALUE_TYPE}) # MPI is required by the library; no extra conditionals needed here \ No newline at end of file diff --git a/bench/bench_builder.cpp b/bench/bench_builder.cpp new file mode 100644 index 0000000..6af1540 --- /dev/null +++ b/bench/bench_builder.cpp @@ -0,0 +1,85 @@ +// Benchmark for cCSRMatrixBuilder::Build() on an MLFMA-like near-field +// assembly pattern (issue #7): leaf boxes of blockSize interacting with a +// (2*nbr+1)^3 stencil neighborhood as dense blocks, plus a configurable +// fraction of re-emitted duplicate triplets. +#include "CSRMatrixBuilder.h" +#include +#include +#include +#include +#include + +using namespace csr4mpi; +using Scalar = double; + +static double NowSec() +{ + using clk = std::chrono::steady_clock; + return std::chrono::duration(clk::now().time_since_epoch()).count(); +} + +static std::vector> MakeNearFieldTriplets(int nBlocks, int blockSize, int nbr, + double duplicateFraction, unsigned seed) +{ + std::mt19937_64 rng(seed); + std::uniform_real_distribution val(-1.0, 1.0); + auto idx = [blockSize](int blk, int i) { return static_cast(blk) * blockSize + i; }; + + std::vector> v; + for (int bi = 0; bi < nBlocks; ++bi) { + for (int bj = std::max(0, bi - nbr); bj <= std::min(nBlocks - 1, bi + nbr); ++bj) { + for (int i = 0; i < blockSize; ++i) + for (int j = 0; j < blockSize; ++j) + v.emplace_back(idx(bi, i), idx(bj, j), static_cast(val(rng))); + } + } + const std::size_t nDup = static_cast(static_cast(v.size()) * duplicateFraction); + std::uniform_int_distribution pick(0, v.size() - 1); + v.reserve(v.size() + nDup); + for (std::size_t k = 0; k < nDup; ++k) { + const auto& t = v[pick(rng)]; + v.emplace_back(t.m_iRow, t.m_iCol, static_cast(val(rng))); + } + return v; +} + +int main(int argc, char** argv) +{ + int nBlocks = 1024; + int blockSize = 16; + int nbr = 13; + double dupFrac = 0.10; + int repeats = 3; + if (argc > 1) nBlocks = std::atoi(argv[1]); + if (argc > 2) blockSize = std::atoi(argv[2]); + if (argc > 3) nbr = std::atoi(argv[3]); + if (argc > 4) dupFrac = std::atof(argv[4]); + if (argc > 5) repeats = std::atoi(argv[5]); + + const iSize n = static_cast(nBlocks) * blockSize; + std::printf("generating near-field triplets: %d blocks x %d, stencil +/-%d, %.0f%% dups ...\n", + nBlocks, blockSize, nbr, dupFrac * 100.0); + const auto triplets = MakeNearFieldTriplets(nBlocks, blockSize, nbr, dupFrac, 2024); + std::printf("triplets: %zu\n", triplets.size()); + + cCSRMatrixBuilder builder; + builder.Reserve(triplets.size()); + builder.AddEntries(triplets); + + double best = 1e300; + std::size_t nnz = 0; + for (int r = 0; r < repeats; ++r) { + double t0 = NowSec(); + auto m = builder.Build(0, n, n, true); + double t1 = NowSec(); + nnz = m.vColInd().size(); + best = std::min(best, t1 - t0); + if (!ValidateCSRMatrix(m)) { + std::printf("ERROR: built matrix failed validation\n"); + return 1; + } + } + std::printf("Build(): nnz=%zu best=%.3fs (%.1fM triplets/s)\n", + nnz, best, static_cast(triplets.size()) / best / 1e6); + return 0; +} diff --git a/include/CSRMatrixBuilder.h b/include/CSRMatrixBuilder.h index 2772697..98ad643 100644 --- a/include/CSRMatrixBuilder.h +++ b/include/CSRMatrixBuilder.h @@ -3,7 +3,6 @@ #include "CSRMatrix.h" #include "Global.h" #include -#include #include namespace csr4mpi { @@ -59,95 +58,108 @@ class cCSRMatrixBuilder { m_vTriplets.clear(); } - // Build the CSR matrix with automatic deduplication + // Build the CSR matrix with automatic deduplication. + // + // Sort-based pipeline (issue #7): the previous implementation deduplicated + // through a global std::unordered_map, which allocated one heap node per + // unique entry (~48-64B), hashed every triplet, and ended with a global + // O(n log n) sort — for near-field assemblies with tens of millions of + // triplets this dominated both build time and peak memory. The pipeline + // below replaces it with: + // + // Pass 1 count in-range entries per row -> rowPtr prefix sum O(N) + // Pass 2 scatter (col, value) into row-segmented slots O(N) + // Pass 3 per-row std::stable_sort by column + merge of equal- + // column runs (accumulate or keep-last) O(sum(nnz_r log nnz_r)) + // Pass 4 final row offsets while emitting the CSR arrays + // + // All passes stream sequentially over memory; no hash nodes, no rehash, + // no global sort. std::stable_sort preserves insertion order within equal + // columns, so duplicates accumulate in the same order (and to the same + // bits) as the previous hash-map implementation. + // // iGlobalRowBegin: starting global row index (typically 0 for a standalone local matrix) // iGlobalRowEnd: ending global row index (exclusive) // iGlobalColCount: total number of columns // bAccumulateDuplicates: if true, accumulate duplicate entries; if false, keep last value cCSRMatrix Build(iSize iGlobalRowBegin, iSize iGlobalRowEnd, iSize iGlobalColCount, bool bAccumulateDuplicates = true) const { - struct cKey { - iIndex r; - iIndex c; - bool operator==(const cKey& other) const { return r == other.r && c == other.c; } - }; - struct cKeyHash { - std::size_t operator()(const cKey& k) const noexcept - { - // Use a robust hash combining function to reduce collisions - // Based on boost::hash_combine algorithm - std::size_t seed = std::hash()(k.r); - seed ^= std::hash()(k.c) + 0x9e3779b9 + (seed << 6) + (seed >> 2); - return seed; - } - }; - - iSize iLocalRows = iGlobalRowEnd - iGlobalRowBegin; - - // Deduplicate using hash map - std::unordered_map mEntries; - mEntries.reserve(m_vTriplets.size()); + const iSize iLocalRows = iGlobalRowEnd - iGlobalRowBegin; + // Pass 1: count in-range entries per row (duplicates included) + std::vector vRowPtr(static_cast(iLocalRows + 1), 0); for (const auto& triplet : m_vTriplets) { // Validate row is within range if (triplet.m_iRow < iGlobalRowBegin || triplet.m_iRow >= iGlobalRowEnd) { continue; // Skip out-of-range entries } - - cKey key { triplet.m_iRow, triplet.m_iCol }; - auto it = mEntries.find(key); - if (it == mEntries.end()) { - mEntries.emplace(key, triplet.m_vValue); - } else { - if (bAccumulateDuplicates) { - it->second += triplet.m_vValue; - } else { - it->second = triplet.m_vValue; // Keep last value - } - } + ++vRowPtr[static_cast(triplet.m_iRow - iGlobalRowBegin + 1)]; } - // Convert to sorted triplets - struct cEntry { - iIndex r; + // Prefix sum -> segment bounds per row + for (iSize r = 0; r < iLocalRows; ++r) { + vRowPtr[static_cast(r + 1)] += vRowPtr[static_cast(r)]; + } + + const std::size_t nTotal = static_cast(vRowPtr.back()); + + // Pass 2: scatter (col, value) pairs into their row segment, + // preserving insertion order within each row + struct cKV { iIndex c; Scalar v; }; - std::vector vSortedEntries; - vSortedEntries.reserve(mEntries.size()); - - for (const auto& p : mEntries) { - vSortedEntries.push_back({ p.first.r, p.first.c, p.second }); + std::vector vTmp(nTotal); + std::vector vCursor(vRowPtr.begin(), vRowPtr.end() - 1); + for (const auto& triplet : m_vTriplets) { + if (triplet.m_iRow < iGlobalRowBegin || triplet.m_iRow >= iGlobalRowEnd) { + continue; + } + const std::size_t r = static_cast(triplet.m_iRow - iGlobalRowBegin); + vTmp[static_cast(vCursor[r]++)] = cKV { triplet.m_iCol, triplet.m_vValue }; } - // Sort by row then column - std::sort(vSortedEntries.begin(), vSortedEntries.end(), [](const cEntry& a, const cEntry& b) { - if (a.r != b.r) - return a.r < b.r; - return a.c < b.c; - }); - - // Build CSR arrays - std::vector vRowPtr(static_cast(iLocalRows + 1), 0); + // Pass 3 + 4: per-row stable sort by column, merge equal-column runs, + // and emit the compacted CSR arrays. Slot r of vRowPtr is overwritten + // with the final offset after its segment [vRowPtr[r], vRowPtr[r+1]) + // has been consumed, so the old prefix sums stay readable while the + // final ones are written in place. std::vector vColInd; std::vector vValues; - - vColInd.reserve(vSortedEntries.size()); - vValues.reserve(vSortedEntries.size()); - - for (const auto& e : vSortedEntries) { - iIndex iLocalRow = e.r - iGlobalRowBegin; - vRowPtr[static_cast(iLocalRow + 1)]++; - vColInd.push_back(e.c); - vValues.push_back(e.v); - } + vColInd.reserve(nTotal); + vValues.reserve(nTotal); - // Prefix sum to get row pointers for (iSize r = 0; r < iLocalRows; ++r) { - vRowPtr[static_cast(r + 1)] += vRowPtr[static_cast(r)]; + const std::size_t s = static_cast(vRowPtr[static_cast(r)]); + const std::size_t e = static_cast(vRowPtr[static_cast(r + 1)]); + + std::stable_sort(vTmp.begin() + s, vTmp.begin() + e, + [](const cKV& a, const cKV& b) { return a.c < b.c; }); + + vRowPtr[static_cast(r)] = static_cast(vColInd.size()); + + std::size_t k = s; + while (k < e) { + std::size_t m = k + 1; + while (m < e && vTmp[m].c == vTmp[k].c) { + ++m; + } + Scalar vAgg = vTmp[k].v; + for (std::size_t q = k + 1; q < m; ++q) { + if (bAccumulateDuplicates) { + vAgg = vAgg + vTmp[q].v; + } else { + vAgg = vTmp[q].v; // Keep last value (stable order -> last inserted) + } + } + vColInd.push_back(vTmp[k].c); + vValues.push_back(vAgg); + k = m; + } } + vRowPtr[static_cast(iLocalRows)] = static_cast(vColInd.size()); - return cCSRMatrix(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount, + return cCSRMatrix(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount, std::move(vRowPtr), std::move(vColInd), std::move(vValues)); } diff --git a/tests/test_csr_builder.cpp b/tests/test_csr_builder.cpp index 37f9b0e..69cc353 100644 --- a/tests/test_csr_builder.cpp +++ b/tests/test_csr_builder.cpp @@ -2,6 +2,10 @@ #include "CSRMatrix.h" #include "Operations.h" #include +#include +#include +#include +#include #include using namespace csr4mpi; @@ -437,15 +441,240 @@ TEST(CSRBuilderTest, BuiltMatrixAlwaysValid) { // Any matrix built with the builder should be valid cCSRMatrixBuilder builder; - + // Add random entries builder.AddEntry(5, 3, 1.0); builder.AddEntry(2, 7, 2.0); builder.AddEntry(8, 1, 3.0); builder.AddEntry(0, 9, 4.0); builder.AddEntry(2, 7, 5.0); // duplicate - + auto matrix = builder.Build(0, 10, 10); - + EXPECT_TRUE(ValidateCSRMatrix(matrix)); } + +// --------------------------------------------------------------------------- +// Bit-exact equivalence against the previous (hash-map) implementation. +// The reference below is a verbatim copy of the pre-issue-7 Build() pipeline; +// it exists only inside the tests to guarantee the sort-based rewrite produces +// identical structure AND identical floating-point bit patterns. +// --------------------------------------------------------------------------- +namespace refimpl { + +template +cCSRMatrix BuildReference(const std::vector>& vTriplets, + iSize iGlobalRowBegin, iSize iGlobalRowEnd, iSize iGlobalColCount, bool bAccumulateDuplicates) +{ + struct cKey { + iIndex r; + iIndex c; + bool operator==(const cKey& other) const { return r == other.r && c == other.c; } + }; + struct cKeyHash { + std::size_t operator()(const cKey& k) const noexcept + { + std::size_t seed = std::hash()(k.r); + seed ^= std::hash()(k.c) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + return seed; + } + }; + + iSize iLocalRows = iGlobalRowEnd - iGlobalRowBegin; + + std::unordered_map mEntries; + mEntries.reserve(vTriplets.size()); + + for (const auto& triplet : vTriplets) { + if (triplet.m_iRow < iGlobalRowBegin || triplet.m_iRow >= iGlobalRowEnd) + continue; + cKey key { triplet.m_iRow, triplet.m_iCol }; + auto it = mEntries.find(key); + if (it == mEntries.end()) { + mEntries.emplace(key, triplet.m_vValue); + } else { + if (bAccumulateDuplicates) + it->second += triplet.m_vValue; + else + it->second = triplet.m_vValue; + } + } + + struct cEntry { + iIndex r; + iIndex c; + Scalar v; + }; + std::vector vSortedEntries; + vSortedEntries.reserve(mEntries.size()); + for (const auto& p : mEntries) + vSortedEntries.push_back({ p.first.r, p.first.c, p.second }); + + std::sort(vSortedEntries.begin(), vSortedEntries.end(), [](const cEntry& a, const cEntry& b) { + if (a.r != b.r) + return a.r < b.r; + return a.c < b.c; + }); + + std::vector vRowPtr(static_cast(iLocalRows + 1), 0); + std::vector vColInd; + std::vector vValues; + vColInd.reserve(vSortedEntries.size()); + vValues.reserve(vSortedEntries.size()); + + for (const auto& e : vSortedEntries) { + iIndex iLocalRow = e.r - iGlobalRowBegin; + vRowPtr[static_cast(iLocalRow + 1)]++; + vColInd.push_back(e.c); + vValues.push_back(e.v); + } + for (iSize r = 0; r < iLocalRows; ++r) { + vRowPtr[static_cast(r + 1)] += vRowPtr[static_cast(r)]; + } + + return cCSRMatrix(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount, + std::move(vRowPtr), std::move(vColInd), std::move(vValues)); +} + +} // namespace refimpl + +template +static void ExpectBitwiseEqual(const cCSRMatrix& a, const cCSRMatrix& b) +{ + ASSERT_EQ(a.iGlobalRowBegin(), b.iGlobalRowBegin()); + ASSERT_EQ(a.iGlobalRowEnd(), b.iGlobalRowEnd()); + ASSERT_EQ(a.iGlobalColCount(), b.iGlobalColCount()); + ASSERT_EQ(a.vRowPtr(), b.vRowPtr()); + ASSERT_EQ(a.vColInd(), b.vColInd()); + ASSERT_EQ(a.vValues().size(), b.vValues().size()); + // Bitwise comparison: identical accumulation order must give identical bits + for (std::size_t k = 0; k < a.vValues().size(); ++k) { + if constexpr (csr4mpi::is_complex_v) { + EXPECT_EQ(std::memcmp(&a.vValues()[k], &b.vValues()[k], sizeof(Scalar)), 0) + << "value bits differ at " << k; + } else { + EXPECT_EQ(a.vValues()[k], b.vValues()[k]) << "value differs at " << k; + } + } +} + +// Randomized triplets with duplicates: accumulate mode +TEST(CSRBuilderTest, BitwiseEqualsReferenceAccumulate) +{ + std::mt19937_64 rng(12345); + std::uniform_int_distribution rowDist(0, 199); + std::uniform_int_distribution colDist(0, 299); + std::uniform_real_distribution valDist(-10.0, 10.0); + + std::vector> triplets; + triplets.reserve(20000); + for (int k = 0; k < 20000; ++k) + triplets.emplace_back(rowDist(rng), colDist(rng), static_cast(valDist(rng))); + + cCSRMatrixBuilder builder; + builder.AddEntries(triplets); + auto got = builder.Build(0, 200, 300, true); + auto want = refimpl::BuildReference(triplets, 0, 200, 300, true); + ExpectBitwiseEqual(got, want); + EXPECT_TRUE(ValidateCSRMatrix(got)); +} + +// Randomized triplets with duplicates: keep-last mode +TEST(CSRBuilderTest, BitwiseEqualsReferenceKeepLast) +{ + std::mt19937_64 rng(54321); + std::uniform_int_distribution rowDist(0, 99); + std::uniform_int_distribution colDist(0, 149); + std::uniform_real_distribution valDist(-10.0, 10.0); + + std::vector> triplets; + triplets.reserve(15000); + for (int k = 0; k < 15000; ++k) + triplets.emplace_back(rowDist(rng), colDist(rng), static_cast(valDist(rng))); + + cCSRMatrixBuilder builder; + builder.AddEntries(triplets); + auto got = builder.Build(0, 100, 150, false); + auto want = refimpl::BuildReference(triplets, 0, 100, 150, false); + ExpectBitwiseEqual(got, want); +} + +// MLFMA-like near-field pattern: block-dense bands + injected duplicates, +// including out-of-range rows that must be filtered identically +TEST(CSRBuilderTest, BitwiseEqualsReferenceNearFieldPattern) +{ + std::mt19937_64 rng(777); + std::uniform_real_distribution valDist(-1.0, 1.0); + const int nBlocks = 64; + const int blockSize = 8; + const int nRows = nBlocks * blockSize; + + auto idx = [](int blk, int i) { return static_cast(blk) * blockSize + i; }; + std::vector> triplets; + for (int bi = 0; bi < nBlocks; ++bi) { + for (int bj = std::max(0, bi - 7); bj <= std::min(nBlocks - 1, bi + 7); ++bj) { + for (int i = 0; i < blockSize; ++i) + for (int j = 0; j < blockSize; ++j) + triplets.emplace_back(idx(bi, i), idx(bj, j), static_cast(valDist(rng))); + } + } + // Re-emit 10% duplicates (exercises accumulation) and add out-of-range rows + std::uniform_int_distribution pick(0, triplets.size() - 1); + const std::size_t nDup = triplets.size() / 10; + triplets.reserve(triplets.size() + nDup + 100); + for (std::size_t k = 0; k < nDup; ++k) { + const auto& t = triplets[pick(rng)]; + triplets.emplace_back(t.m_iRow, t.m_iCol, static_cast(valDist(rng))); + } + for (int k = 0; k < 100; ++k) + triplets.emplace_back(nRows + k, 0, static_cast(valDist(rng))); // out of range + + for (bool accumulate : { true, false }) { + cCSRMatrixBuilder builder; + builder.AddEntries(triplets); + auto got = builder.Build(0, nRows, nRows, accumulate); + auto want = refimpl::BuildReference(triplets, 0, nRows, nRows, accumulate); + ExpectBitwiseEqual(got, want); + EXPECT_TRUE(ValidateCSRMatrix(got)); + } +} + +// Empty input and empty row range remain trivially equal +TEST(CSRBuilderTest, BitwiseEqualsReferenceEdgeCases) +{ + // No entries at all + std::vector> empty; + auto gotA = cCSRMatrixBuilder().Build(0, 5, 5); + auto wantA = refimpl::BuildReference(empty, 0, 5, 5, true); + ExpectBitwiseEqual(gotA, wantA); + + // All entries out of range + std::vector> outOfRange = { + cTriplet(100, 0, static_cast(1.0)), + cTriplet(200, 1, static_cast(2.0)), + }; + auto gotB = cCSRMatrixBuilder().Build(0, 10, 10); + auto wantB = refimpl::BuildReference(outOfRange, 0, 10, 10, true); + ExpectBitwiseEqual(gotB, wantB); + EXPECT_EQ(gotB.vColInd().size(), 0u); +} + +// Complex scalar accumulation order must also match bit-for-bit +TEST(CSRBuilderTest, BitwiseEqualsReferenceComplex) +{ + using CScalar = std::complex; + std::mt19937_64 rng(2468); + std::uniform_int_distribution rowDist(0, 49); + std::uniform_int_distribution colDist(0, 79); + std::uniform_real_distribution valDist(-1.0, 1.0); + + std::vector> triplets; + for (int k = 0; k < 8000; ++k) + triplets.emplace_back(rowDist(rng), colDist(rng), CScalar(valDist(rng), valDist(rng))); + + cCSRMatrixBuilder builder; + builder.AddEntries(triplets); + auto got = builder.Build(0, 50, 80, true); + auto want = refimpl::BuildReference(triplets, 0, 50, 80, true); + ExpectBitwiseEqual(got, want); +}