From ac300ebeaab2a6f3307aae4e893df01486c3aa2c Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 25 May 2026 20:20:44 -0700 Subject: [PATCH 1/3] update CLAUDE.md --- CLAUDE.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CLAUDE.md b/CLAUDE.md index bd549470..17a47240 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -104,7 +104,7 @@ The sparse triangular solve function (`trsm`) handles transposition similarly: - `lskgeX`: **Left** sketching, **ge**neral (dense) data, variant X - `rskgeX`: **Right** sketching, **ge**neral (dense) data, variant X - `lskges`: Left sketching with **s**parse operator -- `lskge3`: Left sketching with dense operator (calls GEMM, "3" for 3-argument GEMM-like) +- `lskge3`: Left sketching with dense operator (calls GEMM, a "BLAS 3" operation) - `lsksp3`: Left sketching **sp**arse data (where "left" refers to operator position) **Counterintuitive detail**: In `lsksp3` and `rsksp3`, the "left/right" refers to the operator's position. But these functions call `right_spmm`/`left_spmm` respectively, where "left/right" refers to the sparse data matrix position. See `sparse_data/DevNotes.md` lines 59-74. From 894712e2fe936b2f06e27d0cdfc1935affabe8ea Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 25 May 2026 20:21:06 -0700 Subject: [PATCH 2/3] initial resolution of GitHub issue #149 --- RandBLAS/sparse_skops.hh | 29 ++++++++++++++++++++++++++--- 1 file changed, 26 insertions(+), 3 deletions(-) diff --git a/RandBLAS/sparse_skops.hh b/RandBLAS/sparse_skops.hh index 8d5c5b61..1a8df3c1 100644 --- a/RandBLAS/sparse_skops.hh +++ b/RandBLAS/sparse_skops.hh @@ -110,6 +110,23 @@ static state_t repeated_fisher_yates( T *vals ) { randblas_error_if(vec_nnz > dim_major); + if (vec_nnz == 1) { + // Sampling 1 element without replacement is the same as with replacement, + // so delegate to the cheaper i.i.d. sampler. Per-matvec call (k=1) matches + // the per-matvec increment formula in compute_next_state. + state_t cur_state = state; + for (int64_t i = 0; i < dim_minor; ++i) { + if (vals != nullptr) { + cur_state = sample_indices_iid_uniform( + dim_major, 1, idxs_major + i, vals + i, cur_state); + } else { + cur_state = sample_indices_iid_uniform( + dim_major, 1, idxs_major + i, cur_state); + } + if (idxs_minor != nullptr) idxs_minor[i] = i; + } + return cur_state; + } std::vector vec_work(dim_major); std::iota(vec_work.begin(), vec_work.end(), 0); std::vector pivots(vec_nnz); @@ -296,9 +313,15 @@ RNGState compute_next_state(SparseDist dist, RNGState state) { int64_t num_mavec, incrs_per_mavec; if (dist.major_axis == Axis::Short) { num_mavec = std::max(dist.n_rows, dist.n_cols); - incrs_per_mavec = dist.vec_nnz; - // ^ SASOs don't try to be frugal with CBRNG increments. - // See repeated_fisher_yates. + if (dist.vec_nnz == 1) { + // SASO with vec_nnz=1 delegates to sample_indices_iid_uniform + // (per-matvec, with rademachers), which is frugal with CBRNG increments. + incrs_per_mavec = (int64_t) std::ceil((double) dist.vec_nnz / ((double) state.len_c/2)); + } else { + incrs_per_mavec = dist.vec_nnz; + // ^ SASOs don't try to be frugal with CBRNG increments. + // See repeated_fisher_yates. + } } else { num_mavec = std::min(dist.n_rows, dist.n_cols); incrs_per_mavec = (int64_t) std::ceil((double) dist.vec_nnz / ((double) state.len_c/2)); From 5be803520addd4ff7e5292ab64c95ee4a1dc59e0 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 25 May 2026 21:12:40 -0700 Subject: [PATCH 3/3] revise sample_indices_iid_uniform and use it as a short-circuit for repeated_fisher_yates with vec_nnz==1 --- RandBLAS/sparse_skops.hh | 53 ++++++++++++---------------------------- RandBLAS/util.hh | 38 +++++++++++++--------------- 2 files changed, 33 insertions(+), 58 deletions(-) diff --git a/RandBLAS/sparse_skops.hh b/RandBLAS/sparse_skops.hh index 1a8df3c1..7a1afdc8 100644 --- a/RandBLAS/sparse_skops.hh +++ b/RandBLAS/sparse_skops.hh @@ -48,10 +48,6 @@ namespace RandBLAS::sparse { -inline std::uint64_t promote_uint_pair(std::uint32_t &a, std::uint32_t &b) { - return static_cast(a) + (static_cast(b) << 32); -} - template > void _considerate_fisher_yates( const state_t &state, @@ -112,20 +108,16 @@ static state_t repeated_fisher_yates( randblas_error_if(vec_nnz > dim_major); if (vec_nnz == 1) { // Sampling 1 element without replacement is the same as with replacement, - // so delegate to the cheaper i.i.d. sampler. Per-matvec call (k=1) matches - // the per-matvec increment formula in compute_next_state. - state_t cur_state = state; - for (int64_t i = 0; i < dim_minor; ++i) { - if (vals != nullptr) { - cur_state = sample_indices_iid_uniform( - dim_major, 1, idxs_major + i, vals + i, cur_state); - } else { - cur_state = sample_indices_iid_uniform( - dim_major, 1, idxs_major + i, cur_state); - } - if (idxs_minor != nullptr) idxs_minor[i] = i; + // so delegate to the cheaper i.i.d. sampler in a single batched call. + if (idxs_minor != nullptr) + std::iota(idxs_minor, idxs_minor + dim_minor, sint_t{0}); + if (vals != nullptr) { + return sample_indices_iid_uniform( + dim_major, dim_minor, idxs_major, vals, state); + } else { + return sample_indices_iid_uniform( + dim_major, dim_minor, idxs_major, state); } - return cur_state; } std::vector vec_work(dim_major); std::iota(vec_work.begin(), vec_work.end(), 0); @@ -310,26 +302,13 @@ inline state_t repeated_fisher_yates( template RNGState compute_next_state(SparseDist dist, RNGState state) { - int64_t num_mavec, incrs_per_mavec; - if (dist.major_axis == Axis::Short) { - num_mavec = std::max(dist.n_rows, dist.n_cols); - if (dist.vec_nnz == 1) { - // SASO with vec_nnz=1 delegates to sample_indices_iid_uniform - // (per-matvec, with rademachers), which is frugal with CBRNG increments. - incrs_per_mavec = (int64_t) std::ceil((double) dist.vec_nnz / ((double) state.len_c/2)); - } else { - incrs_per_mavec = dist.vec_nnz; - // ^ SASOs don't try to be frugal with CBRNG increments. - // See repeated_fisher_yates. - } - } else { - num_mavec = std::min(dist.n_rows, dist.n_cols); - incrs_per_mavec = (int64_t) std::ceil((double) dist.vec_nnz / ((double) state.len_c/2)); - // ^ LASOs do try to be frugal with CBRNG increments. - // See sample_indices_iid_uniform. - } - int64_t full_incr = num_mavec * incrs_per_mavec; - state.counter.incr(full_incr); + // Both _considerate_fisher_yates (SASO with vec_nnz > 1) and + // sample_indices_iid_uniform (SASO with vec_nnz == 1, and LASO) consume + // exactly one CBRNG counter increment per nonzero. + int64_t num_major_axis_vec = (dist.major_axis == Axis::Short) + ? std::max(dist.n_rows, dist.n_cols) + : std::min(dist.n_rows, dist.n_cols); + state.counter.incr(num_major_axis_vec * dist.vec_nnz); return state; } diff --git a/RandBLAS/util.hh b/RandBLAS/util.hh index 597b1e55..dcff6b5a 100644 --- a/RandBLAS/util.hh +++ b/RandBLAS/util.hh @@ -513,35 +513,31 @@ state_t sample_indices_iid(int64_t n, const T* cdf, int64_t k, sint_t* samples, return state_t(ctr, key); } +inline std::uint64_t promote_uint_pair(std::uint32_t a, std::uint32_t b) { + return static_cast(a) + (static_cast(b) << 32); +} + template > state_t sample_indices_iid_uniform(int64_t n, int64_t k, sint_t* samples, T* rademachers, const state_t &state) { - using RNG = typename state_t::generator; - auto [ctr, key] = state; - RNG gen; - auto rv_array = r123ext::uneg11::generate(gen, ctr, key); - int64_t len_c = static_cast(state.len_c); if constexpr (WriteRademachers) { - len_c = 2*(len_c/2); - // ^ round down to the nearest multiple of two. + randblas_require(state.len_c >= 4); + } else { + randblas_require(state.len_c >= 2); } - int64_t rv_index = 0; - double dN = (double) n; + using RNG = typename state_t::generator; + RNG gen; + auto ctr = state.counter; + auto key = state.key; + std::uint64_t n_64 = static_cast(n); for (int64_t i = 0; i < k; ++i) { - auto random_unif01 = uneg11_to_u01(rv_array[rv_index]); - sint_t sample_index = (sint_t) dN * random_unif01; - samples[i] = sample_index; - rv_index += 1; + auto rv = gen(ctr, key); + ctr.incr(); + std::uint64_t s = promote_uint_pair(rv[0], rv[1]); + samples[i] = static_cast(s % n_64); if constexpr (WriteRademachers) { - rademachers[i] = (rv_array[rv_index] >= 0) ? (T) 1 : (T) -1; - rv_index += 1; - } - if (rv_index == len_c) { - ctr.incr(1); - rv_array = r123ext::uneg11::generate(gen, ctr, key); - rv_index = 0; + rademachers[i] = (rv[2] % 2 == 0) ? (T) 1 : (T) -1; } } - if (0 < rv_index) ctr.incr(1); return state_t(ctr, key); }