Sparse sampling performance benchmarks and parallelism

[rileyjmurray]

The following is written by a robot. PR changes were reviewed by me (Riley) extensively over the course of several rounds. The agent-driven work was by GPT 5.6 Sol. I ran adversarial review through Fable 5 at high effort.

---

This PR began as an effort to bring sparse sketching-operator generation in line with dense generation: SASOs and LASOs should use the ambient OpenMP configuration without making their samples depend on the number of threads. It now resolves #139 and includes a broader cleanup of sampling benchmarks, dimension validation, empty-submatrix behavior, and sparse matrix multiplication edge cases.

Sampling semantics

Sparse sampling assigns randomness to logical major-axis vectors rather than to physical threads. If the first requested vector starts at counter c, then vector i starts at c + i * vec_nnz and writes to a deterministic lane. The returned state advances by num_major_axis_vectors * vec_nnz, independent of the number of active threads or the OpenMP schedule.

Changing the OpenMP thread count does not change vals, rows, cols, nnz, or the returned RNGState. This applies to full and submatrix requests, tall and wide operators, and both choices of major axis. A build without OpenMP retains the serial behavior.

No execution-policy or thread-count argument was added. RandBLAS consults the ambient OpenMP configuration and applies an internal policy that may select fewer threads when there is too little work to amortize parallel startup or per-thread workspace.

Main changes

Parallel SASO sampling

repeated_fisher_yates now processes major-axis vectors in parallel. The vec_nnz == 1 path reduces to batched uniform sampling and does not allocate permutation workspace. For vec_nnz > 1, every active thread owns a restored permutation of length dim_major and a pivot array of length vec_nnz. This keeps the sampling work at O(num_major_axis_vectors * vec_nnz), but increases permutation workspace from O(dim_major) to O(T * dim_major) for T active threads.

The same OpenMP implementation handles both the single-threaded and multithreaded cases. The sparse-operator path asks repeated_fisher_yates to sort each completed major-axis vector before its sampling lane ends, while the public support-only function preserves the original unsorted order. This removes a separate sorting pass without changing its output.

Parallel LASO sampling

LASO vectors are sampled with replacement into fixed-width lanes. Each thread uses private hash maps to merge duplicate locations, sorts the surviving entries within its vector, and records the survivor count. A serial pass then packs those lanes in increasing major-axis-vector order. Every destination is at or before its source, so the pack is safe in place and preserves canonical COO order without allocating another O(nnz) buffer.

Allocation failures during parallel merging are captured inside the OpenMP region and rethrown after the team exits. I deliberately left a parallel prefix-sum and scatter out of this PR. The serial pack is simpler, preserves the existing representation, and should only be replaced if end-to-end measurements show that it is material.

Benchmark and testing infrastructure

This PR adds saso_sampling_performance, following the conventions in sketch_general_performance.cc. It compares repeated_fisher_yates with five plausible C++ implementations:

• std::sample over an iota-filled vector;
• partial Fisher-Yates with a restored length-n population;
• std::shuffle followed by the first k entries;
• uniform draws with duplicate rejection; and
• Floyd's algorithm with std::unordered_set.

The natural comparison uses std::mt19937_64 for those implementations and Philox for RandBLAS. A second comparison adapts Philox to the standard URBG interface. That adapter uses one 64-bit result per counter by default, matching RandBLAS' sampling stream, and can optionally consume both 64-bit results from each Philox counter.

The baseline samplers now live in RandBLAS/testing/samplers.hh; future discrete-uniform baseline samplers have a common home there. Shared benchmark parsing and OpenMP-setting helpers live in RandBLAS/testing/benchmarking.hh. The std::sample baseline uses an ordinary iota-filled vector rather than a ranges sentinel pair, which keeps the code portable to MSVC.

Comparison mode forces RandBLAS to one OpenMP thread because each competing implementation owns one serial engine. Scaling mode times RandBLAS alone and reports both requested and effective thread counts. sketch_general_performance also has a separate sparse-sampling scaling section alongside its existing sparse-operator application measurements.

Correctness and portability hardening

The implementation now checks products that determine matrix sizes, output capacities, workspace lengths, and pointer offsets before using them. Those checks cover dense and sparse operator generation. They are intentionally about addressable matrix data; wrap-around of the extended-width unsigned CBRNG counter remains accepted behavior.

Submatrix bounds checking is centralized in validate_submat_dims and reused by dense sketching, sparse sketching, and sparse matrix operations. Boundary windows with a zero dimension are valid. offset_and_ldim still reports a leading dimension of at least one for such matrices, as required by BLAS conventions, and sparse submatrix materialization can now return an empty COO matrix.

randblas_get_thread_num and randblas_get_num_threads provide the OpenMP queries RandBLAS needs while returning serial answers when OpenMP is absent. This lets the sampling and sparse-matrix kernels use one code path in both builds.

A small layout-aware lascl utility now centralizes in-place dense matrix scaling. This helper was adapted from Max Melnichenko's work on the add-symm-kernels branch. The MKL sparse wrapper uses it when an empty contraction dimension or an all-zero sparse operand leaves beta * C; the RandBLAS sparse dispatcher uses the same operation before its own accumulation kernels. This avoids depending on MKL's treatment of sparse empty-matrix edge cases.

Performance

I ran the required SASO sweep on an Apple M3 MacBook Air with Clang 19.1.3 and OpenMP 5.1. The problem size was dim_major = 2000, num_major_axis_vectors = 100000, with 10 trials for each vec_nnz.

vec_nnz	1-thread median	8-thread minimum	minimum-time speedup
1	610 us	182 us	3.31x
2	1,716 us	471 us	3.64x
4	3,226 us	906 us	3.56x
8	6,070 us	1,386 us	4.37x
16	11,975 us	2,662 us	4.49x

Every row passed the exact output/state check. I did not observe a one-thread regression relative to the pre-implementation baseline; the final medians were lower for all five values of vec_nnz.

I also ran the sampling section of sketch_general_performance for both major-axis orientations at vec_nnz in {1, 4, 16}. The 8-thread minimum-time speedups ranged from 2.12x to 4.19x, and every configuration passed its exact output/state check.

Tests and documentation

The tests compare exact sparse buffers and returned states across thread counts. They cover the specialized and general Fisher-Yates paths; full, strict-submatrix, and empty-boundary requests; tall and wide SASOs and LASOs; nondefault value and index types; thread-policy boundaries; overflow checks; and validation errors that must be reported before parallel execution.

Separate tests exercise the baseline samplers, their Philox adapter, the benchmark parser, dense and sparse submatrix validation, BLAS-compatible leading dimensions, and the sparse multiplication cases that leave beta * C. The more specialized tests now include comments explaining the behavior under test and how the chosen setup exposes it.

The documentation describes logical-vector counter partitioning, SASO workspace costs, LASO packing, accepted CBRNG wrap-around, and the distinction between random-stream keys and physical thread assignment. The testing developer notes document the shared sampler and benchmark helpers. The style guide was also clarified while bringing the changed files into compliance.

Size of the change

The table below groups the origin/main...HEAD diff into nonoverlapping directories. The full PR changes 34 files, with +3,057 / -297 lines.

Directory	Insertions (+)	Deletions (-)
RandBLAS/ (top level)	+486	-168
RandBLAS/sparse_data/	+43	-61
RandBLAS/testing/	+368	0
examples/ (top level)	+11	0
examples/simple-kernel-benchmarks/	+928	-49
examples/sparse-low-rank-approx/	+7	-2
test/ (top level)	+184	-2
test/basic_rng/	+387	0
test/datastructures/	+304	-2
test/linops/	+202	-2
test/meta/	+84	0
rtd/	+6	-4
Repository root	+47	-7

[rileyjmurray]

Reviewing this PR has me thinking about what the safe_int_product helper function buys us.

[rileyjmurray]

Leaving a comment to confirm that I've reviewed all tests.

[rileyjmurray]

Leaving tons of comments to prove that I actually reviewed each file. Some comments request changes.