Skip to content

perf(builder): replace hash-map dedup with sort-based CSR assembly #8

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
deltaeecs merged 1 commit into from
Sep 7, 2026
Merged

perf(builder): replace hash-map dedup with sort-based CSR assembly #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions bench/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
85 changes: 85 additions & 0 deletions bench/bench_builder.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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 <algorithm>
#include <chrono>
#include <cstdio>
#include <random>
#include <vector>

using namespace csr4mpi;
using Scalar = double;

static double NowSec()
{
using clk = std::chrono::steady_clock;
return std::chrono::duration<double>(clk::now().time_since_epoch()).count();
}

static std::vector<cTriplet<Scalar>> MakeNearFieldTriplets(int nBlocks, int blockSize, int nbr,
double duplicateFraction, unsigned seed)
{
std::mt19937_64 rng(seed);
std::uniform_real_distribution<double> val(-1.0, 1.0);
auto idx = [blockSize](int blk, int i) { return static_cast<iIndex>(blk) * blockSize + i; };

std::vector<cTriplet<Scalar>> 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<Scalar>(val(rng)));
}
}
const std::size_t nDup = static_cast<std::size_t>(static_cast<double>(v.size()) * duplicateFraction);
std::uniform_int_distribution<std::size_t> 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<Scalar>(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<iSize>(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<Scalar> 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<double>(triplets.size()) / best / 1e6);
return 0;
}
142 changes: 77 additions & 65 deletions include/CSRMatrixBuilder.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@
#include "CSRMatrix.h"
#include "Global.h"
#include <algorithm>
#include <unordered_map>
#include <vector>

namespace csr4mpi {
Expand Down Expand Up @@ -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<Scalar> 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<iIndex>()(k.r);
seed ^= std::hash<iIndex>()(k.c) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
return seed;
}
};

iSize iLocalRows = iGlobalRowEnd - iGlobalRowBegin;

// Deduplicate using hash map
std::unordered_map<cKey, Scalar, cKeyHash> mEntries;
mEntries.reserve(m_vTriplets.size());
const iSize iLocalRows = iGlobalRowEnd - iGlobalRowBegin;

// Pass 1: count in-range entries per row (duplicates included)
std::vector<iIndex> vRowPtr(static_cast<std::size_t>(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<std::size_t>(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<std::size_t>(r + 1)] += vRowPtr[static_cast<std::size_t>(r)];
}

const std::size_t nTotal = static_cast<std::size_t>(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<cEntry> vSortedEntries;
vSortedEntries.reserve(mEntries.size());

for (const auto& p : mEntries) {
vSortedEntries.push_back({ p.first.r, p.first.c, p.second });
std::vector<cKV> vTmp(nTotal);
std::vector<iIndex> 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<std::size_t>(triplet.m_iRow - iGlobalRowBegin);
vTmp[static_cast<std::size_t>(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<iIndex> vRowPtr(static_cast<std::size_t>(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<iIndex> vColInd;
std::vector<Scalar> vValues;

vColInd.reserve(vSortedEntries.size());
vValues.reserve(vSortedEntries.size());

for (const auto& e : vSortedEntries) {
iIndex iLocalRow = e.r - iGlobalRowBegin;
vRowPtr[static_cast<std::size_t>(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<std::size_t>(r + 1)] += vRowPtr[static_cast<std::size_t>(r)];
const std::size_t s = static_cast<std::size_t>(vRowPtr[static_cast<std::size_t>(r)]);
const std::size_t e = static_cast<std::size_t>(vRowPtr[static_cast<std::size_t>(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<std::size_t>(r)] = static_cast<iIndex>(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<std::size_t>(iLocalRows)] = static_cast<iIndex>(vColInd.size());

return cCSRMatrix<Scalar>(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount,
return cCSRMatrix<Scalar>(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount,
std::move(vRowPtr), std::move(vColInd), std::move(vValues));
}

Expand Down
Loading
Loading