Add SYMM kernels for sparse symmetric matrices; rework sketch_symmetric to exploit symmetry by mmelnich · Pull Request #163 · BallisticLA/RandBLAS

mmelnich · 2026-05-13T23:41:50Z

Summary

Adds SYMM-shaped kernels to RandBLAS so that symmetric matrices can be sketched with their symmetry exploited. Implements all four (A storage, B storage) cases, and folds in repo-wide cleanups uncovered along the way (a lascl utility, helper extractions, a shared MKL call wrapper).

Motivation: the funNystrom++ work in RandLAPACK. sketch_symmetric was silently forwarding to sketch_general (plain GEMM), so A's symmetry was not exploited and a sparse SkOp had to be densified before the call. The current funNystrom++ implementation works around the gap on its sparse-sketch path by populating BOTH triangles of A and applying the sketch through right_spmm, which reads A as a general dense matrix. Case B of this PR provides exactly the missing kernel, so a small RandLAPACK follow-up can switch that call site and keep one-triangle storage (this halves the memory footprint for large explicit symmetric matrices). Tracked in BallisticLA/RandLAPACK#167.

A condensed 4-case table and the MKL availability matrix live in RandBLAS/sparse_data/DevNotes.md.

Status (2026-08-24): rebased over #196/#197 and reviewed against the new governance docs. Changes in this round, per review:

No API break anymore. The four pre-Uplo sketch_symmetric overloads (trailing sym_check_tol) are restored with their exact original behavior (runtime symmetry check, both triangles read, sketch_general forwarding) alongside the new blas::Uplo overloads; overload resolution cannot collide since blas::Uplo is a scoped enum. The old symmetry-check test is back, a new test pins legacy and Uplo results against each other, and the web docs advertise the legacy variants under a marked dropdown.
Contracts aligned with Sparse sampling performance benchmarks and parallelism #196's left_spmm: beta passes through to MKL (restoring mkl_spgemm_to_dense's alpha=1/beta=0 direct-write fast path in Case D) with lascl only before the hand kernels; empty products (zero dimension, empty operand, alpha = 0) leave beta*Y instead of tripping MKL's handle-creation rejection of empty matrices; window bounds go through validate_submat_dims; layout flips use flipped_layout; index template parameters use the SignedInteger concept.
Style pass against STYLE_GUIDE.md / CONTRIBUTING.md and Sparse sampling performance benchmarks and parallelism #196's conventions: snake_case behavior-phrase test names; internal kernels documented with plain comments (Doxygen /// reserved for rendered declarations); Symmetric<SpMat>/as_symmetric get API-reference directives; the perf benchmark verifies every timed method against a reference (PASS/FAIL column), adopts the shared testing/benchmarking.hh helpers (--help, --threads, requested-vs-effective reporting), and describes methods by behavior rather than PR chronology.
Case D test coverage completed to all 9 (A, B) format pairings plus a non-CSR side=Right cell; an int32-index Case D test pins the non-default-index path.
util.hh now carries @rileyjmurray's lascl from Sparse sampling performance benchmarks and parallelism #196 verbatim; the branch's remaining lascl adoptions in trsm_dispatch.hh and two examples are the sweep requested in the LASCAL review thread. One note on mkl_left_spmm: the runtime return false for unsupported scalar types from Sparse sampling performance benchmarks and parallelism #196 is folded into the shared mkl_sparse_mm_call, whose static_assert fires at compile time instead; the runtime branch was unreachable because make_mkl_handle already rejects such types at compile time.

Earlier hardening (previous round, still in effect): left_spmm-style validation on every sparse entry point with error-path tests, side=Right normalized at the dispatcher via the layout-flip identity (the hand-written Right kernels are deleted), column-driven OpenMP-parallel fallback kernels, csc_spsymm as a transpose-view delegation, sketch_symmetric overload count reduced with a readable static_assert for unsupported operator kinds, window sampling via submatrix_as_coo for unmaterialized operators, and Symmetric<SpMat> rejecting temporaries at compile time.

The four cases

Tag	Operation	A storage	B storage	Status
A	dense-symm x dense	dense, one triangle	dense	Implemented via `blas::symm` in `sksy.hh`
B	dense-symm x sparse	dense, one triangle	sparse	Implemented. `lsksys` / `rsksys` wrappers in `sksy.hh` (validation, beta, window sampling) over the column-driven `coo_lsksys` kernel in `sparse_data/coo_sksys_impl.hh`; `coo_rsksys` reduces to it via the transpose identity. Only A's named triangle is read
C	sparse-symm x dense, dense result	sparse, one triangle	dense	Implemented in `sparse_data/spsymm_dispatch.hh`. side=Right normalized at entry; MKL fast path covers all of `{CSR, CSC, COO}`; column-driven OpenMP-parallel hand kernels as the non-MKL fallback
D	sparse-symm x sparse, dense result	sparse, one triangle	sparse	Implemented in `sparse_data/spsymm_dispatch.hh`: expand-A (O(nnz)) + reuse of the existing spgemm-to-dense path on MKL builds; densify-B + Case-C composition as the non-MKL fallback

What's in the PR

The Case-C dispatcher (spsymm_dispatch.hh) follows the left_spmm shape: normalize side=Right away, validate, leave beta*Y for empty products, try MKL with the caller's beta (MKL fuses it), otherwise apply beta via lascl and run a pure-accumulator hand kernel. The MKL wrapper (mkl_spsymm_impl.hh) applies the SPARSE_MATRIX_TYPE_SYMMETRIC descriptor; CSC is consumed as a CSR-of-transpose view over the same buffers with uplo flipped, since A equals its transpose. A runtime SPARSE_STATUS_NOT_SUPPORTED (a parameter-validation result, so Y is untouched) falls back to the hand kernels.

The hand kernels (csr_spsymm_impl.hh / coo_spsymm_impl.hh) are column-driven accumulators: #pragma omp parallel for over the n right-hand-side columns, a scan of the stored triangle inside, one symmetric-read resolution per entry (read the mirrored position when the index pair falls outside the named triangle). Off-diagonal entries contribute twice per column, diagonal entries once; entries outside the named triangle are silently skipped, so a caller that stored both triangles still gets the correct answer. csc_spsymm_impl.hh is a delegation to the CSR kernel on the transpose view.

Case D (spsymm_dispatch.hh + expand_symmetric_to_general in symmetric.hh): the expansion walks the stored triangle twice (count, then fill) and emits each off-diagonal entry mirrored, producing an owning general COO in O(nnz) memory that the existing mkl_spgemm_to_dense accepts with a GENERAL descriptor. MKL cannot be handed the symmetric matrix directly: mkl_sparse_sp2m rejects SPARSE_MATRIX_TYPE_SYMMETRIC descriptors and mkl_sparse_?_spmmd takes no descriptor at all. side=Right reduces by the same layout-flip identity with B.transpose() as a lightweight view. Non-MKL builds and index widths mismatched with MKL_INT fall back to densifying B and composing through Case C.

Case B (sksy.hh + coo_sksys_impl.hh): sketch_symmetric dispatches by SkOp kind inside one definition per call shape. DenseSkOp goes through dense::lsksy3/rsksy3 (SYMM-backed, with a shared transpose-copy helper for layout mismatch so the SYMM speedup is retained). SparseSkOp goes through sparse::lsksys/rsksys, which validate like their dense siblings, apply beta ahead of the accumulator kernel, sample only the requested window for unmaterialized operators, and call the column-driven COO kernel.

Symmetric<SpMat> wrapper (sparse_data/symmetric.hh): non-owning carrier holding const SpMat& and blas::Uplo. Square-matrix invariant at construction; deleted const-rvalue overloads reject temporaries at compile time. Type-safety guard: a Symmetric<SpMat> argument fails to bind to general spmm/spgemm. Public RandBLAS::spsymm has both a raw (SpMat+uplo) overload and a wrapper-taking overload.

Cleanups along the way:

RandBLAS::util::lascl (util.hh): A := alpha*A for a dense m x n matrix, LAPACK-style naming. Replaces hand-rolled "loop-over-blas::scal" beta-scale patterns across trsm_dispatch.hh, spmm_dispatch.hh, mkl_spmm_impl.hh, and two example programs. After this sweep plus main's own kernel rework, util::safe_scal has no callers left inside RandBLAS (it stays; RandLAPACK uses it).
internal::mkl_sparse_mm_call (mkl_spmm_impl.hh): type-dispatched wrapper around mkl_sparse_d_mm/mkl_sparse_s_mm, shared by mkl_left_spmm (GENERAL descriptor) and mkl_spsymm (SYMMETRIC descriptor).
Index-generic *_to_dense: the six coo/csr/csc_to_dense helpers now deduce the index type (previously bound to int64_t only), with the default preserved.
<numeric> include in sparse_data/base.hh: latent gap (used std::iota but had relied on transitive inclusion); the new symmetric.hh exposed it.
RandBLAS.hh umbrella: pulls in sparse_data/spsymm_dispatch.hh so downstream #include <RandBLAS.hh> consumers see RandBLAS::spsymm and the Symmetric<SpMat> wrapper automatically.
Tests and the perf example reuse library utilities (overwrite_triangle, symmetrize, coo_to_dense, dense_to_csr) instead of hand-rolled conversion loops.

Perf benchmark: examples/simple-kernel-benchmarks/spsymm_performance.cc compares (a) sparse: RandBLAS::spsymm (one-triangle storage + MKL fast path) vs the pre-PR workaround (both triangles + RandBLAS::spmm) vs a dense blas::symm reference, and (b) dense: the new sketch_symmetric (SYMM-backed) vs the equivalent sketch_general call (the pre-PR GEMM-forwarding behaviour). Single-config CLI or default 3-point sweep.

API: extension, not a break

sketch_symmetric gains blas::Uplo overloads; the original sym_check_tol overloads remain with their exact pre-existing behavior (O(n^2) runtime symmetry check via require_symmetric, both triangles read, forwarding to sketch_general). The Uplo overloads read only the named triangle and skip the runtime check; new code should prefer them. The web docs present both, with the legacy variants marked as retained for compatibility.

Test coverage

24-cell {COO, CSR, CSC} x {ColMajor, RowMajor} x {Upper, Lower} x {Left, Right} for spsymm Case C (in test/linops/test_spsymm.cc), plus float/beta=0/alpha=0/wrapper-routing edge cases; 7 Case D correctness tests plus an int32-index Case D test; 5 cases for the Symmetric<SpMat> wrapper; and test_sketch_symmetric.cc covers 8 DenseSkOp paths plus SparseSkOp paths across both sides, both layouts, both uplo, a lift configuration, and an unmaterialized-operator configuration (exercising the window-sampling path). Error paths are now tested: leading-dimension violations, submatrix windows exceeding the operator, and one-based indices all assert a thrown RandBLAS::Error.

Tolerance for spsymm and sparse sketch_symmetric tests: atol = 100*eps, rtol = 10*eps. The dense reference (blas::symm on a fully-symmetrized A) and the sparse kernels accumulate FMAs in a different order, giving a few-ULP divergence. The relaxed tolerance documents this is expected, not a bug.

Verified locally after the second rebase and review pass, on Linux + GCC 15.2 + MKL 2026.0 sparse: all symm suites pass (21/21 TestSpsymm, 18/18 TestSketchSymmetric incl. the restored legacy tests, 5/5 TestSymmetricWrapper) and the full ctest suite passes.

mmelnich · 2026-05-13T23:43:23Z

@rileyjmurray, this is still a WIP, but feel free to look.

rileyjmurray · 2026-05-14T16:25:34Z

Regarding LASCAL, there are a few other sites where it can be used that weren't in your most recent commit. You can find some hits in the search

https://github.com/search?q=repo%3ABallisticLA%2FRandBLAS++RandBLAS%3A%3Autil%3A%3Asafe_scal&type=code.

See also

RandBLAS/RandBLAS/sparse_data/mkl_spmm_impl.hh

Lines 379 to 383 in d8a153f

    
               for (int64_t j = 0; j < n; ++j) 
        
                   blas::scal(m, beta, &C[j * ldc], 1); 
        
           } else { 
        
               for (int64_t i = 0; i < m; ++i) 
        
                   blas::scal(n, beta, &C[i * ldc], 1);

.

…LAS.hh Three small changes: 1. New benchmark examples/simple-kernel-benchmarks/spsymm_performance.cc. Mirrors spmm_performance.cc in spirit but answers two perf questions specific to PR #163: (a) Sparse: RandBLAS::spsymm (one-triangle storage + MKL fast path via SPARSE_MATRIX_TYPE_SYMMETRIC) vs. the pre-PR workaround (both triangles stored, called through RandBLAS::spmm + MKL). Also includes a dense blas::symm reference as an "ideal SYMM" baseline. (b) Dense: the rewritten RandBLAS::sketch_symmetric (SYMM-backed) vs. the equivalent RandBLAS::sketch_general call (the pre-PR GEMM-forwarding behaviour). Single-config CLI: ./spsymm_performance n_A d density [num_trials] Default sweep: n_A in {500, 1000, 2000}, d=200, density=0.05, num_trials=10. Reports median + min over trials per kernel. 2. examples/CMakeLists.txt: register the new spsymm_performance target alongside spmm_performance. 3. RandBLAS.hh: add #include <RandBLAS/sparse_data/spsymm_dispatch.hh> to the umbrella header so downstream code using #include <RandBLAS.hh> gets RandBLAS::spsymm and the Symmetric<SpMat> wrapper visible automatically. (Without this, the benchmark and any other downstream consumer would need to know to include the internal dispatch header directly --- inconsistent with how spmm, spgemm, sketch_symmetric etc. are exposed.) Verified: - Standalone g++ compile of spsymm_performance.cc against the in-tree headers passes clean. - Main RandBLAS build clean after the umbrella header change. - 45 / 45 focused tests pass (TestSpsymm: 11, TestSpGEMM: 21, TestSymmetricWrapper: 5, TestSketchSymmetric: 8).

mmelnich · 2026-05-14T22:53:10Z

@rileyjmurray lmk your thoughts

…LAS.hh Three small changes: 1. New benchmark examples/simple-kernel-benchmarks/spsymm_performance.cc. Mirrors spmm_performance.cc in spirit but answers two perf questions specific to PR #163: (a) Sparse: RandBLAS::spsymm (one-triangle storage + MKL fast path via SPARSE_MATRIX_TYPE_SYMMETRIC) vs. the pre-PR workaround (both triangles stored, called through RandBLAS::spmm + MKL). Also includes a dense blas::symm reference as an "ideal SYMM" baseline. (b) Dense: the rewritten RandBLAS::sketch_symmetric (SYMM-backed) vs. the equivalent RandBLAS::sketch_general call (the pre-PR GEMM-forwarding behaviour). Single-config CLI: ./spsymm_performance n_A d density [num_trials] Default sweep: n_A in {500, 1000, 2000}, d=200, density=0.05, num_trials=10. Reports median + min over trials per kernel. 2. examples/CMakeLists.txt: register the new spsymm_performance target alongside spmm_performance. 3. RandBLAS.hh: add #include <RandBLAS/sparse_data/spsymm_dispatch.hh> to the umbrella header so downstream code using #include <RandBLAS.hh> gets RandBLAS::spsymm and the Symmetric<SpMat> wrapper visible automatically. (Without this, the benchmark and any other downstream consumer would need to know to include the internal dispatch header directly --- inconsistent with how spmm, spgemm, sketch_symmetric etc. are exposed.) Verified: - Standalone g++ compile of spsymm_performance.cc against the in-tree headers passes clean. - Main RandBLAS build clean after the umbrella header change. - 45 / 45 focused tests pass (TestSpsymm: 11, TestSpGEMM: 21, TestSymmetricWrapper: 5, TestSketchSymmetric: 8).

mmelnich · 2026-08-13T22:18:38Z

@rileyjmurray I kind of forgot about this PR. Could you take a look again?

mmelnich · 2026-08-23T19:24:52Z

@rileyjmurray Anything else you'd like me to do here?

rileyjmurray · 2026-08-24T05:16:51Z

Please resolve merge conflicts. I stole the lascl utility in a PR I just merged. After that, please revert any API-breaking changes. You can extend the API with function overloading as needed.

Please have an agent do a comprehensive review of the PR with reference to CONTRIBUTORS.md, AGENTS.md, and STYLE_GUIDE.md. On the topic of style, have it make inferences or draw comparisons to the PR I just merged: 49099b8.

Phase 1 of the symm-kernels plan (project-plans/randblas-symm-plan.md): implements Case A (dense-symmetric A x dense Omega via blas::symm) and adds stubbed signatures for Case B (dense-symmetric A x sparse Omega) that throw RandBLAS::Error pointing at the plan doc. Previously, sketch_symmetric in sksy.hh checked symmetry of A via require_symmetric, then forwarded directly to sketch_general -> lskge3 / rskge3 -> blas::gemm. Symmetry of A was never exploited: both triangles had to be stored and validated to match, and the SYMM 1.3-1.8x speedup over GEMM was unrealized. New helpers in RandBLAS::dense (sksy.hh): lsksy3(layout, uplo, d, n, alpha, S, ro_s, co_s, A, lda, beta, B, ldb) Computes B = alpha * submat(S) * mat(A) + beta * B with mat(A) n-by-n symmetric (only the triangle named by uplo is read). rsksy3(layout, uplo, n, d, alpha, A, lda, S, ro_s, co_s, beta, B, ldb) Computes B = alpha * mat(A) * submat(S) + beta * B. Both follow the lskge3 / rskge3 materialization pattern (submatrix_as_blackbox when S.buff is null). When the buffered S's storage layout matches the caller's layout, the final call is blas::symm with side = Right (lsksy3) or side = Left (rsksy3). When layouts differ, SYMM cannot transpose S on the fly; the call falls back to blas::gemm with opS = Trans -- correct but loses the SYMM speedup on that path. Optimization target: transpose-copy of S into matching layout, tracked in randblas-symm-plan.md. Refactor of sketch_symmetric (sksy.hh): - All four overloads now accept blas::Uplo uplo and forward it to lsksy3 / rsksy3. - Each overload is split into a DenseSkOp specialization (dispatching to the helpers) and a SparseSkOp specialization that throws RandBLAS::Error via randblas_require(false, ...). The throw cites the Case-B kernel from the plan doc; this locks the public API so future PRs can fill in the body without breaking source compatibility for downstream users. - Removed the require_symmetric call and the sym_check_tol parameter. With SYMM dispatch only the triangle named by uplo is read, so the runtime symmetry check is no longer meaningful and the O(n^2) scan is wasted work. util.hh: require_symmetric remains available as a standalone validator. Its docstring notes that sketch_symmetric no longer calls it post-Phase-1. Test updates (test/linops/test_sketch_symmetric.cc): - sketch_symmetric_side, test_same_layouts, test_opposing_layouts now accept a blas::Uplo parameter (default Uplo::Upper, matching random_symmetric_mat's symmetrize-from-upper convention). - Removed test_error_on_asymmetric and its TEST_F entry; that test verified the require_symmetric throw that no longer fires. API break for downstream callers of sketch_symmetric: must add a blas::Uplo argument. The only known in-tree consumer is RandLAPACK's ExplicitSymLinOp (via funnystrompp PR #132); a small follow-up patch on the RandLAPACK side will thread uplo through. Cases B and D bodies (the harder hand-roll kernels) remain stubs; see project-plans/randblas-symm-plan.md for the rationale. Verified: 449 / 449 ctest pass on Linux + GCC 13.3 + CUDA-aware blaspp + MKL sparse. Drop from 450 to 449 is the removed symmetry_check_fails_for_asymmetric_matrix TEST_F.

Phase 2 of the symm-kernels plan (project-plans/randblas-symm-plan.md): a lightweight non-owning wrapper that marks a SparseMatrix as symmetric with a stored-triangle annotation. The wrapper is the carrier for symmetric sparse matrices into the spsymm-family kernels (Phase 3 onward). Design: - Symmetric<SpMat> holds const SpMat& A and const blas::Uplo uplo. - Re-exposes A's scalar_t and index_t aliases at the wrapper level. - Square-matrix invariant (A.n_rows == A.n_cols) enforced at construction via randblas_require. - Non-owning: caller keeps A alive for the wrapper's lifetime. - Re-exported into the RandBLAS:: namespace per the existing sparse_data:: aliasing pattern. Type-system role: a Symmetric<SpMat> argument fails to bind to the general spmm / spgemm templates, preventing accidental "treat symmetric as general" bugs. spsymm-family kernels added in Phase 3 will accept the wrapper (or, equivalently, raw SpMat + uplo at the BLAS-mirror boundary). Files: - RandBLAS/sparse_data/symmetric.hh (new, 95 lines) - test/datastructures/test_symmetric_wrapper.cc (new, 92 lines, 5 cases: construction from COO/CSR/CSC, namespacing access at RandBLAS:: scope, non-square reject) - test/CMakeLists.txt: register the new test in SPARSEDATA_SOURCES. Incidental fix: - RandBLAS/sparse_data/base.hh now #include <numeric>. base.hh uses std::iota at line 178 but relied on transitive inclusion via coo_matrix.hh. The new symmetric.hh includes base.hh more directly and exposed the gap. Verified: 454 / 454 ctest pass on Linux + GCC 13.3 + CUDA-aware blaspp + MKL sparse. (Prior Phase 1: 449 / 449. The 5 new tests come from test_symmetric_wrapper.cc.)

…r-format fallbacks Phase 3 + Phase 4 of the symm-kernels plan (project-plans/randblas-symm-plan.md): the Case C kernel (sparse-symmetric A x dense B -> dense Y) with the full 3-format x 2-layout x 2-uplo x 2-side dispatch grid, plus the matching 24-cell test surface. Public API (in RandBLAS::): spsymm(layout, uplo, m, n, alpha, A, B, ldb, beta, Y, ldy) Convenience wrapper. Defaults to side=Left (Y = alpha * A * B + beta * Y). spsymm(layout, m, n, alpha, Symmetric<SpMat> A_sym, B, ldb, beta, Y, ldy) Wrapper-routing overload. Extracts uplo from the Symmetric<SpMat> carrier (Phase 2). Avoids the "raw SpMat + separate uplo" pattern at the call site. Lower-level dispatch entry (in RandBLAS::sparse_data::): spsymm(layout, side, uplo, m, n, alpha, A, B, ldb, beta, Y, ldy) Full BLAS-mirror signature with the side flag. Tries the MKL fast path when available, falls back to the per-format kernel otherwise. Dispatch behavior: - With RandBLAS_HAS_MKL and matching index width, mkl_spsymm calls mkl_sparse_d_mm with descr.type = SPARSE_MATRIX_TYPE_SYMMETRIC and descr.mode mapped from uplo. Free MKL win for the common case. - mkl_spsymm signals fallback (returns false) for side=Right (MKL has no side parameter on sparse-mm) and for CSC format (mkl_sparse_d_mm returns NOT_SUPPORTED on CSC even with a symmetric descriptor; same constraint as mkl_left_spmm). - Hand-rolled fallbacks in csr_spsymm_impl.hh, csc_spsymm_impl.hh, coo_spsymm_impl.hh handle the remaining cells. Each emits one blas::axpy for the stored entry and a second one for the implied symmetric counterpart on off-diagonal entries; the diagonal contributes once. Both layouts handled by switching axpy strides. - Entries outside the named triangle are silently skipped, so callers who store both triangles by mistake still get correct results (the kernel just behaves as if the wrong-side entries didn't exist). Helpers: - internal::apply_beta_scale in csr_spsymm_impl.hh: shared beta-scaling pass on Y, included by csc/coo to avoid duplication. - Reuses make_mkl_handle / to_mkl_layout / check_mkl_status from mkl_spmm_impl.hh. Tests (test/linops/test_spsymm.cc, added to SPARSEDATA_SOURCES): Ten TEST_F entries covering the 24-cell grid: {COO, CSR, CSC} x {ColMajor, RowMajor} x {Upper, Lower} x {Left, Right} Each (format, side) test internally sweeps (layout, uplo) via SCOPED_TRACE. Plus float-precision coverage, beta=0 edge case, alpha=0 edge case, and the Symmetric<SpMat> wrapper-routing case. Reference: dense blas::symm on the fully-symmetrized A. spsymm input has only the named triangle stored (other triangle zeroed before conversion). Tolerance: atol = 100*eps, rtol = 10*eps. The dense and sparse paths accumulate FMAs in different orders (off-diagonal entries contribute via two AXPYs in the sparse path vs. one fused dotprod in dense SYMM), yielding a few ULPs of divergence. Initial run failed at strict 1*eps rtol with absDiff ~ 12 ULPs on side=Right and CSC cases; relaxed tolerance documents this is expected, not a bug. Verified: 10 / 10 spsymm tests pass locally. Full ctest sweep deferred to CI. Build clean on Linux + GCC 13.3 + CUDA-aware blaspp + MKL sparse; the two pre-existing test_exceptions.cc sign-compare warnings are unrelated. Stubs for cases B (dense-symm x sparse) and D (sparse-symm x sparse -> dense) remain throw-only (Phase 5). See randblas-symm-plan.md for the rationale.

…se design Phase 5 of the symm-kernels plan (project-plans/randblas-symm-plan.md): locks the API surface for Case D and documents the broader Symm-kernel landscape in the sparse-data DevNotes. New stub in spsymm_dispatch.hh: template <SparseMatrix SpMatA, SparseMatrix SpMatB, typename T = ...> void spsymm(layout, side, uplo, m, n, alpha, const SpMatA& A, const SpMatB& B, beta, T* Y, ldy); The body is a single randblas_require(false, ...) that throws RandBLAS::Error with a verbose message: names the case, links back to the plan doc, and lists the two composition fallbacks callers can use today (densify B and call the Case-C spsymm, or call mkl_sparse_sp2m with a symmetric descriptor on A and densify the resulting sparse C). Case B was already stubbed in Phase 1 (the SparseSkOp branch of sketch_symmetric in sksy.hh); no new code needed there. DevNotes (RandBLAS/sparse_data/DevNotes.md): added a "SYMM-shaped kernels (spsymm)" section before the existing "Sketching sparse data with dense operators" section. It includes: - A 4-case table mapping operand storage (dense-symm vs sparse-symm on the left, dense vs sparse on the right) to implementation status this PR. - An MKL availability matrix explaining why each case is / isn't natively supported (MKL has no Side parameter for sparse-mm, mkl_sparse_d_mm returns NOT_SUPPORTED on CSC, no MKL routine for dense-symm-A x sparse-B, sp2m is the closest to Case D but writes sparse output). - The Case C dispatch flow (validate -> MKL try -> per-format fallback) with the internal::apply_beta_scale helper noted. - "Cases B and D: why stub-only" subsection documenting the surveyed 2026-05 finding that no portable kernel exists in MKL / Ginkgo / SparseBLAS/spblas-reference / MAGMA-sparse, and explaining why hand-rolling each is non-trivial (case B has an awkward gather-column-from-triangle pattern; case D layers a triangle filter on top of spgemm-into-dense). Test: new CaseD_SparseSparseThrows TEST_F in test_spsymm.cc verifies the Case D overload throws RandBLAS::Error with the expected behavior. Verified: 11 / 11 spsymm tests pass locally (10 from Phase 4 + the new Case D throw test).

Phase 6 of the symm-kernels plan (project-plans/randblas-symm-plan.md): ReadTheDocs updates to reflect the Phase 1 API change to sketch_symmetric and the Phase 3 / 5 introduction of spsymm. rtd/source/api_reference/sketch_dense.rst: Replaced all four RandBLAS::sketch_symmetric doxygenfunction references with the new (blas::Uplo uplo)-taking signature; dropped the sym_check_tol parameter from the references. Added a brief preamble noting that SparseSkOp now throws (Case B of the SYMM-kernels plan; the link points at sparse_data/DevNotes.md). rtd/source/api_reference/sketch_sparse.rst: New dropdown for RandBLAS::spsymm covering both the public-API signature (raw SpMat + uplo) and the Symmetric<SpMat>-overload signature. Includes a paragraph on the dispatch flow (MKL fast path for side=Left non-CSC; per-format fallback otherwise) and a companion-stubs note for Cases B and D. rtd/source/FAQ.rst: Replaced the stale "Symmetric matrices have to be stored as general matrices" entry (which claimed sketch_symmetric worked equally well with DenseSkOp and SparseSkOp) with two new entries: - "SparseSkOp is not yet supported in sketch_symmetric": explains that the SparseSkOp branch is the Case B stub, names the composition fallback (densify the sparse op), and points at the surveyed-2026-05 finding that no portable kernel exists. - "Layout-mismatched DenseSkOp in sketch_symmetric falls back to GEMM": documents the first-cut layout-mismatch handling. No CHANGELOG file exists in the repo, so nothing to add there. Build clean. 24 / 24 focused tests (TestSpsymm*, TestSymmetricWrapper*, TestSketchSymmetric*) pass after these doc-only changes.

…e sites Introduces a single LAPACK-named utility for the "Y := alpha * Y" pattern on a dense m-by-n matrix, then replaces every hand-rolled site in RandBLAS that previously open-coded the same shape. New helper (RandBLAS/util.hh): template <typename T> inline void lascl(blas::Layout layout, int64_t m, int64_t n, T alpha, T* A, int64_t lda); Named after LAPACK's ?lascl for the convention; the signature is deliberately simpler than LAPACK's (no CFROM / CTO overflow protection, no matrix-type flag --- general dense only). Sits right after the existing vector-form safe_scal in the same namespace. Fast paths: alpha == 1 returns immediately, alpha == 0 uses std::fill, otherwise per-column (ColMajor) or per-row (RowMajor) blas::scal. Replaced sites: - sparse_data/csr_spsymm_impl.hh: removed the local internal::apply_beta_scale helper (introduced in the spsymm commit earlier this PR); callers now use RandBLAS::util::lascl directly. - sparse_data/csc_spsymm_impl.hh, coo_spsymm_impl.hh: dropped the include of csr_spsymm_impl.hh that pulled in the local helper; each now includes util.hh directly. - sparse_data/mkl_spmm_impl.hh, mkl_spgemm_to_dense: - alpha == 0 path: the entire if/elseif zero-or-scal block collapses to one lascl call. - !direct_write path: per-vector blas::scal + blas::axpy loop is now a single lascl(layout, m, n, beta, C, ldc) followed by the axpy loop. Mathematically equivalent. Verified: 24 / 24 focused tests (TestSpsymm: 11, TestSymmetricWrapper: 5, TestSketchSymmetric: 8) and 21 / 21 TestSpGEMM tests pass after the refactor. No behavioral change; only de-duplication. Net diff: +53 / -56 lines across 5 files.

The three spsymm fallback kernels differ only in their outer iteration order (CSR walks by row, CSC by column, COO by nnz-triple); the per-stored-entry scatter --- two blas::axpy calls covering (i, j) and the implied symmetric (j, i) when i != j, with layout-aware strides --- is identical across all three. Previously this scatter was duplicated ~30 lines per format file (60 lines counting the side=Left + side=Right branches), ~90 lines total. New file RandBLAS/sparse_data/spsymm_internal.hh in namespace RandBLAS::sparse_data::internal: template <typename T> inline void spsymm_scatter_left (blas::Layout layout, int64_t n, T av, int64_t i, int64_t j, const T* B, int64_t ldb, T* Y, int64_t ldy); template <typename T> inline void spsymm_scatter_right(blas::Layout layout, int64_t m, T av, int64_t i, int64_t j, const T* B, int64_t ldb, T* Y, int64_t ldy); Each emits one axpy for the structural entry and a second one for the implied symmetric counterpart (when i != j). The caller has already folded alpha into av; the caller has already filtered out entries outside the uplo'd triangle. csr / csc / coo_spsymm_impl.hh now each include spsymm_internal.hh and replace their 8-line layout-branching inner block with one call to the appropriate scatter helper. Net reduction: ~60 lines across the three format files. Verified: 11 / 11 TestSpsymm tests pass; no behavioral change (the extracted helpers contain the exact same code in the same order).

Previously mkl_spsymm bailed (returned false) on any CSC input, mirroring mkl_left_spmm's NOT_SUPPORTED-on-CSC constraint --- callers fell back to the hand-rolled csc_spsymm kernel. For symmetric A, we can do better: since A == A^T, the CSC.transpose() view is a lightweight reinterpretation of the same buffers as a CSR matrix representing the same A. MKL accepts CSR, so we can stay on the fast path by recursing on the transpose view. Subtlety: when CSC is reinterpreted as CSR, the structurally stored triangle flips. A CSC entry at (i, j) with i <= j (Upper) appears in the CSR view at (j, i) with j >= i, i.e., the Lower triangle of the CSR view. The recursive call therefore flips uplo: CSC Upper -> CSR Lower and vice versa. side=Right still falls back; the limitation there is structural (MKL's mkl_sparse_d_mm has no Side parameter), not format-specific. The csc_spsymm hand kernel remains in the codebase --- it's still the dispatch target for non-MKL builds, and for the side=Right + CSC combination which doesn't pass through this fast path. Verified: 11 / 11 TestSpsymm tests pass. TestSpsymm.CSC_Left now exercises the MKL fast path; CSC_Right still exercises the fallback kernel.

The four SparseSkOp specializations of sketch_symmetric had identical stub bodies (~10 lines each: a fan of (void) casts on the parameters to silence unused-parameter warnings, followed by a randblas_require with the same multi-line message). Replaces them with a single shared helper in a new namespace RandBLAS::detail: template <typename ...Args> [[noreturn]] inline void throw_sketch_symmetric_case_b(Args&&...); The variadic parameter pack consumes the caller's arguments, so the compiler does not warn about unused parameters in the (deliberately unused, by design) stub bodies --- no (void) cast wall needed. The [[noreturn]] attribute communicates to the optimizer that the helper does not return; randblas_require throws RandBLAS::Error, which satisfies the contract. __builtin_unreachable() after the throw is a belt-and-suspenders hint for compilers that do not see through the exception path. Each of the four SparseSkOp stubs becomes a one-line forward: detail::throw_sketch_symmetric_case_b(layout, uplo, ..., B, ldb); Net reduction: ~32 lines (4 stubs * ~8 lines per stub) replaced by 22 lines of helper + 4 one-line forwards. Behavior is identical: still throws RandBLAS::Error with the same composition-fallback message. Verified: 24 / 24 focused tests pass. TestSpsymm.CaseD_SparseSparseThrows and the dedicated SparseSkOp throw scenarios continue to fire correctly.

SymmetricWrapper TEST_F previously duplicated ~30 lines of build-A, build-B, dense-reference, compare setup from run_case --- only the final RandBLAS::spsymm call differed (took a Symmetric<SpMat> wrapper instead of a raw SpMat + uplo). Adds an optional `bool route_via_wrapper = false` parameter to run_case. When true (side=Left only, since the public wrapper overload defaults side=Left), the dispatch block goes through RandBLAS::spsymm(layout, m, n, alpha, Symmetric<SpMat>, ...) instead of the lower-level RandBLAS::sparse_data::spsymm. All other setup is shared. The SymmetricWrapper TEST_F body becomes a single run_case call with route_via_wrapper=true. Net reduction: ~30 lines (the entire duplicated setup block). No coverage loss: the wrapper-routing path still exercises the Symmetric<SpMat> overload at the public namespace and the as_symmetric sugar. Verified: 11 / 11 TestSpsymm tests pass, including SymmetricWrapper.

…'s mm Both mkl_left_spmm and mkl_spsymm previously open-coded the same if-constexpr branch on T to pick between mkl_sparse_d_mm and mkl_sparse_s_mm, with the same shape of arguments (op, alpha, A_handle, descr, layout, B, n_rhs, ldb, beta, C, ldc). Two near-identical 18-line blocks across two files. Replaces with a single template helper in mkl_spmm_impl.hh: template <typename T> inline sparse_status_t mkl_sparse_mm_call( sparse_operation_t op, T alpha, sparse_matrix_t A_handle, const struct matrix_descr& descr, sparse_layout_t mkl_layout, const T* B, int64_t n_rhs, int64_t ldb, T beta, T* C, int64_t ldc); The caller controls op, descr (GENERAL vs SYMMETRIC), and the post-call status interpretation. The helper just dispatches on T and returns the status. mkl_left_spmm: collapses the if-constexpr block to one call. Status check unchanged. mkl_spsymm: same. Keeps the NOT_SUPPORTED-fallback branch unchanged (distinct from mkl_left_spmm's behavior --- mkl_spsymm wants to silently fall back, mkl_left_spmm wants to throw). Net reduction: ~25 lines across two files. Behavior identical; mkl_sparse_d_mm / mkl_sparse_s_mm are called with the same args, just routed through one place. Verified: 37 / 37 MKL-sensitive tests pass (TestSpsymm: 11, TestSpGEMM: 21, TestSymmetricWrapper: 5).

After the initial lascl extraction (5676e13) replaced the three sites introduced by the symm work, four pre-existing sites elsewhere in the repo were doing the same loop-over-safe_scal matrix-scale pattern. Sweeping them now. Sites updated: - RandBLAS/sparse_data/trsm_dispatch.hh The B := alpha * B prelude in the public `trsm` template was a layout-branched loop calling safe_scal per row/column. Collapsed to one lascl(layout, m, n, alpha, B, ldb). - RandBLAS/sparse_data/spmm_dispatch.hh The C := beta * C prelude in left_spmm has the same shape on the output operand. Collapsed to lascl(layout, d, n, beta, C, ldc). - examples/simple-kernel-benchmarks/spmm_performance.cc The handrolled-spmm reference in this micro-benchmark mirrors spmm_dispatch's beta-scale block; updated for consistency so the benchmark and the dispatch path call the same helper. - examples/sparse-low-rank-approx/qrcp_matrixmarket.cc Stage 2 zeroed Q one column at a time inside the column-scatter loop. Hoist the zero-out into a single lascl(ColMajor, m, k, 0.0, Q, ldq) call before the loop. Same net effect; the scatter loop becomes a clean per-column copy with no inline zeroing. The vector-form safe_scal call in test_denseskop.cc:244 (safe_scal(n_srows * n_scols, 0.0, smat) on a contiguous buffer of total length n_srows*n_scols) is genuinely 1-D and stays as safe_scal. util.hh includes added where they were previously pulled in transitively (trsm_dispatch.hh, spmm_dispatch.hh). Verified: 465 / 465 ctest pass on the full suite. No behavioral change; lascl produces the same Y := alpha * Y result as the loop form, and the qrcp_matrixmarket hoist-then-scatter is equivalent to zero-each-column-just-before-scattering it.

Previously mkl_spsymm bailed (returned false) on side=Right because mkl_sparse_d_mm has no Side parameter --- the sparse matrix is always on the left of the dense block. For symmetric A, a layout-flip transformation lets us still hit the MKL fast path: Y = alpha * B * A + beta * Y (side=Right, A is n_A-by-n_A) Take the transpose of both sides: Y^T = alpha * A^T * B^T + beta * Y^T = alpha * A * B^T + beta * Y^T (A symmetric, A == A^T) In MKL terms with side=Left semantics: input is B^T (n_A-by-m), output is Y^T (n_A-by-m). We get this view of the user's B and Y buffers by telling MKL the opposite layout from the user's: user ColMajor B (m-by-n_A, ldb >= m) reinterpreted as RowMajor (n_A-by-m, leading dim = ldb) -> MKL sees buffer[r*ldb + c] = (RowMajor view at row r, col c) = user_buffer[c + r*ldb] = user_B[c, r] = (B^T)[r, c] correctly. The same calculation works for user RowMajor (flip to ColMajor). ldb and ldy carry through unchanged because the reinterpretation has the same leading-dim semantics in either direction. The number of MKL-side right-hand-side columns is m (the user's row count) rather than n. opA stays NoTrans because A^T = A. CSC handling unchanged: still recurses via A.transpose() CSR-view with flipped uplo. With this change, CSC + side=Right also hits MKL (the recursive call enters the new side=Right branch on the CSR view). Net effect on the test grid: TestSpsymm.{CSR,CSC,COO}_Right all now exercise the MKL fast path; the hand-rolled fallback for side=Right remains in place for non-MKL builds. Verified: 11 / 11 TestSpsymm tests pass.

Previously lsksy3 / rsksy3 fell back to blas::gemm with opS=Trans when the buffered DenseSkOp's storage layout differed from the caller's requested layout. Correct, but it gave up the 1.3-1.8x SYMM speedup that's the whole reason this path exists. Replaces the fallback with a transpose-copy + SYMM: 1. Allocate a tight std::vector<T> of size d*n (lsksy3) or n*d (rsksy3) in the caller's layout. 2. util::omatcopy from S's buffer with S.layout-derived (irs_in, ics_in) strides into the temp buffer with caller-layout strides. 3. blas::symm(side, uplo, d, n, alpha, A, lda, S_copy, lds_new, beta, B, ldb). Cost: one O(d*n) memory pass for the copy. Wins over GEMM-with-Trans whenever the SYMM speedup on the d*n*n_A matvec exceeds the d*n copy --- effectively always once n_A is more than a handful, since the matvec is quadratic in n_A. In practice DenseSkOp.layout is normally set to match the consumer's layout, so this path is rare; this commit just removes a sharp edge where layout mismatch silently halved the throughput. Verified: 8 / 8 TestSketchSymmetric tests pass, including the test_opposing_layouts cases that exercise the new transpose-copy path explicitly.

…LAS.hh Three small changes: 1. New benchmark examples/simple-kernel-benchmarks/spsymm_performance.cc. Mirrors spmm_performance.cc in spirit but answers two perf questions specific to PR #163: (a) Sparse: RandBLAS::spsymm (one-triangle storage + MKL fast path via SPARSE_MATRIX_TYPE_SYMMETRIC) vs. the pre-PR workaround (both triangles stored, called through RandBLAS::spmm + MKL). Also includes a dense blas::symm reference as an "ideal SYMM" baseline. (b) Dense: the rewritten RandBLAS::sketch_symmetric (SYMM-backed) vs. the equivalent RandBLAS::sketch_general call (the pre-PR GEMM-forwarding behaviour). Single-config CLI: ./spsymm_performance n_A d density [num_trials] Default sweep: n_A in {500, 1000, 2000}, d=200, density=0.05, num_trials=10. Reports median + min over trials per kernel. 2. examples/CMakeLists.txt: register the new spsymm_performance target alongside spmm_performance. 3. RandBLAS.hh: add #include <RandBLAS/sparse_data/spsymm_dispatch.hh> to the umbrella header so downstream code using #include <RandBLAS.hh> gets RandBLAS::spsymm and the Symmetric<SpMat> wrapper visible automatically. (Without this, the benchmark and any other downstream consumer would need to know to include the internal dispatch header directly --- inconsistent with how spmm, spgemm, sketch_symmetric etc. are exposed.) Verified: - Standalone g++ compile of spsymm_performance.cc against the in-tree headers passes clean. - Main RandBLAS build clean after the umbrella header change. - 45 / 45 focused tests pass (TestSpsymm: 11, TestSpGEMM: 21, TestSymmetricWrapper: 5, TestSketchSymmetric: 8).

…impl Replaces the SparseSkOp-branch throw with a hand-rolled kernel, matching the Case-A pattern in spirit but adapted to read only the named triangle of A and walk the COO entries of the SkOp. New helpers in RandBLAS::sparse (sksy.hh): template <typename T, typename RNG, SignedInteger sint_t> void lsksys(layout, uplo, d, n, alpha, S, ro_s, co_s, A, lda, beta, B, ldb) // B = alpha * submat(S) * mat(A) + beta * B // S is d-by-n SparseSkOp (submat view); A is n-by-n dense symm. template <typename T, typename RNG, SignedInteger sint_t> void rsksys(layout, uplo, n, d, alpha, A, lda, S, ro_s, co_s, beta, B, ldb) // B = alpha * mat(A) * submat(S) + beta * B // A is n-by-n dense symm; S is n-by-d SparseSkOp. Inner loop: for each COO triple (row_S, col_S, v) of the sparse SkOp, filtered inline by the (ro_s, co_s, d, n) submatrix window: - lsksys: contribute alpha*v to row (row_S - ro_s) of B from row (col_S - co_s) of the symmetric A. - rsksys: contribute alpha*v to column (col_S - co_s) of B from column (row_S - ro_s) of the symmetric A. Reading a row or column of a one-triangle-stored symmetric matrix splits into two contiguous ranges based on the diagonal, so the per-stored-entry body is exactly two blas::axpy calls per branch. Uplo flips which range comes from the "stored side" and which from the "transposed read", and layout flips the AXPY strides; that gives four Uplo x Layout branches per side, all using the same two-AXPY structure. Materialization-if-needed: same `if (S.nnz < 0) { shallowcopy + fill_sparse + recurse }` pattern as lskges / rskges. lascl handles the beta scaling on entry. Wired all four SparseSkOp specializations of sketch_symmetric to dispatch to these helpers; removed the dead detail::throw_sketch_symmetric_case_b helper. Tests (test/linops/test_sketch_symmetric.cc): new test_sparse_skop helper that builds a SparseSkOp, densifies it into a reference dense buffer matching the requested layout, then compares sketch_symmetric's output against blas::symm on the densified reference. Six new TEST_F entries: sparse_skop_left_colmajor_upper sparse_skop_left_rowmajor_upper sparse_skop_right_colmajor_upper sparse_skop_right_rowmajor_upper sparse_skop_lower_triangle (4 cells: side x layout, Lower-only) sparse_skop_lift (2 cells: lift directions) Same 100*eps / 10*eps tolerance as the spsymm tests --- the two-axpy scatter accumulates FMAs in a different order than dense SYMM. Docs: - RandBLAS/sparse_data/DevNotes.md: Case B row in the 4-case table now says "Implemented via hand-rolled lsksys / rsksys"; the "why stub-only" subsection rewritten as a "Case B: hand-rolled" subsection describing the access pattern; "Case D: stub-only" survives as the standalone deferred-work item. - rtd/source/FAQ.rst: replaced "SparseSkOp is not yet supported" entry with one explaining both branches are now supported (and how each dispatches). - rtd/source/api_reference/sketch_dense.rst: removed "throws on SparseSkOp" claim; describes the new dispatch. - rtd/source/api_reference/sketch_sparse.rst: Companion-stubs note now mentions Case B as implemented, only Case D as throw-stub. Verified: 14/14 TestSketchSymmetric tests pass (8 prior DenseSkOp + 6 new SparseSkOp). 51/51 focused tests pass overall (TestSpsymm: 11, TestSpGEMM: 21, TestSymmetricWrapper: 5, TestSketchSymmetric: 14). Only Case D remains as a stub after this commit.

The two-SparseMatrix-arg spsymm overload (sparse-symm A x sparse B -> dense Y) now allocates an m-by-n std::vector<T> for B_dense in the caller's layout, fills it via the format-specific coo_to_dense / csr_to_dense / csc_to_dense helper picked by if constexpr, and forwards to the existing Case-C spsymm overload on the densified buffer. Covers all 3 x 3 = 9 sparse-format pairings for (A, B); works in MKL and non-MKL builds. Why composition rather than a single MKL call: mkl_sparse_sp2m returns SPARSE_STATUS_NOT_SUPPORTED when descrA.type == SPARSE_MATRIX_TYPE_SYMMETRIC (only GENERAL is accepted there); mkl_sparse_d_spmmd takes no descriptor at all. So the symmetric expansion has to happen on the RandBLAS side either way -- composing through Case C gets it for free at the cost of an O(m*n) temporary, small for the typical workload where B is a sketching operator with nnz(B) << m*n. Tests: 7 new TEST_F entries in test/linops/test_spsymm.cc (CSR-CSR, CSC-CSC, COO-COO, mixed format, side=Right, float, alpha=0/beta-scale) replace the prior CaseD_SparseSparseThrows. Reference: dense blas::symm on a fully-symmetrized A and a densified B, same 100*eps / 10*eps tolerance as the existing Case C tests. Docs: DevNotes 4-case table + MKL-availability row + Case-D section updated; rtd/source/api_reference/sketch_sparse.rst dropdown note updated to drop the "stub" framing. Verified locally: 477/477 ctest pass on Linux + GCC 13.3 + CUDA-aware blaspp + MKL sparse.

…e docs 1. Factor the per-stored-nonzero two-AXPY scatter body out of sksy.hh's lsksys / rsksys into RandBLAS::sparse_data::coo_lsksys / coo_rsksys in the new RandBLAS/sparse_data/coo_sksys_impl.hh. The wrappers in sksy.hh are now thin glue: handle SparseSkOp materialization and recurse, util::lascl the output, unpack the COO view, then call the kernel. Puts the format-specific work next to the other COO kernels under sparse_data, where Riley flagged it should live. 2. Drop the platform-specific "~1.3-1.8x over GEMM" multiplier from the layout-mismatch transpose-copy comment in sksy.hh. The comment now just says "keeps the SYMM speedup over GEMM" without a number. RandBLAS is platform-agnostic and the multiplier was a back-of-envelope estimate, not a measurement. 3. Fix stale FAQ entry claiming sketch_symmetric falls back to GEMM on layout-mismatched DenseSkOp. The actual behavior since 6d0ca8c is a transpose-copy of S into the caller's layout (O(d * n)) followed by blas::symm. Updated rtd/source/FAQ.rst to reflect that. 4. Fix stale sketch_symmetric doxygen note in sksy.hh claiming SparseSkOp throws. Since 77a3b45 (Case B promotion), SparseSkOp dispatches to the hand-rolled lsksys / rsksys path. Updated the docstring on the submat overload to describe the actual dispatch. 5. DevNotes Case B section, MKL-availability table, and 4-case table updated to point at the new file location and describe the wrapper-vs-kernel split. Verified: 477/477 ctest pass on Linux + GCC 13.3 + CUDA-aware blaspp + MKL sparse.

Three locations had ASCII double-dashes in body text where commas, periods, or parentheses read more naturally: RandBLAS/sparse_data/DevNotes.md Case B MKL-availability row: hand-roll explanation rephrased. Symmetric<SpMat> wrapper bullets: " -- routes via ..." commas. RandBLAS/sparse_data/coo_sksys_impl.hh Header comment describing the two-axpy split: colon plus period rephrase, no semantic change. rtd/source/api_reference/sketch_sparse.rst Case D dropdown note: period plus new sentence in place of the sentence-internal dash. No semantic change; 477/477 ctest pass unchanged.

The Windows CI lanes (added on main after this branch was cut) compile the test suite with MSVC, which does not define __PRETTY_FUNCTION__. The portable macro in RandBLAS/testing/comparison.hh already dispatches to __FUNCSIG__ on MSVC; these three call sites predate it. All four windows-msvc lanes failed on exactly this.

coo/csr/csc_to_dense were bound to the default int64_t index type, so any code routing a non-int64-indexed matrix through them failed to compile. Deduce sint_t (default preserved); no caller changes needed.

…o_general Symmetric<SpMat> stores const SpMat&, so wrapping a temporary (notably the by-value view returned by .transpose()) dangled at the end of the statement. Deleted const-rvalue overloads on the constructor and as_symmetric turn that into a compile error; const&& catches const and non-const rvalues without the forwarding-reference trap a plain && overload would create. expand_symmetric_to_general walks the stored triangle twice (count, fill) and emits off-diagonal entries mirrored, producing an owning general COO in O(nnz) memory. This is the bridge to consumers that only accept general sparse matrices (notably MKL's sparse-times-sparse routines, which reject SYMMETRIC descriptors); Case D uses it in the next commit.

…rnels; Case D via expand-A The dispatcher now does what left_spmm does, in the same order: normalize side=Right to side=Left at entry (A == A^T, so Y = B*A is Y^T = A*B^T with the B/Y buffers reinterpreted in the flipped layout; uplo unchanged), validate (zero-based indices, square A of order m, ldb/ldy lower bounds), apply beta exactly once via util::lascl, return early on alpha == 0, try MKL with beta = 1, then fall back to pure-accumulator hand kernels. The single beta application also removes a double-scaling hazard on MKL's runtime NOT_SUPPORTED fallback (that status is a parameter-validation result, so Y is untouched when it fires). The hand kernels are rewritten column-driven: an OpenMP-parallel outer loop over the n dense right-hand-side columns (each column owned by one thread, race-free), a scan of the stored triangle inside, and one symmetric-read resolution per entry instead of the per-uplo-per-layout grid of strided two-axpy range splits. Unit-stride in ColMajor. This also eliminates the out-of-bounds pointer formation the old zero-length axpy arms performed for entries in the last row/column. csc_spsymm is now a three-line delegation to csr_spsymm on the transpose view with uplo flipped (the identity the MKL path already used); the Right-side branches and spsymm_scatter_{left,right} are deleted, and spsymm_internal.hh with them. mkl_spsymm loses its side parameter (side is normalized before it is reached) and its header comment now states the real fallback contract. Case D no longer densifies B on the primary path: on MKL builds with index widths matching MKL_INT, A's stored triangle is expanded to a general sparse matrix in O(nnz) memory (expand_symmetric_to_general) and the existing mkl_spgemm_to_dense runs with a GENERAL descriptor, keeping B sparse. The densify-B + Case-C composition survives only as the non-MKL / index-width-mismatch fallback, and now compiles for any index type via the index-generic *_to_dense helpers. side=Right reduces by the same layout-flip identity with B.transpose() as a lightweight view.

…e overload grid lsksys/rsksys now carry the same randblas_require set as their dense siblings lsksy3/rsksy3 (SkOp window bounds and lda/ldb lower bounds); previously the sparse branch validated nothing, so a bad leading dimension corrupted memory silently and an oversized window silently returned a sketch with missing rows where the dense branch threw. For an unmaterialized SparseSkOp, the wrappers sample only the requested window via submatrix_as_coo (the lskges pattern) instead of materializing the entire operator and filtering per nonzero; a materialized operator still goes through coo_view_of_skop with in-kernel filtering. coo_lsksys is rewritten column-driven (OpenMP-parallel over the n columns of B, one symmetric-read resolution per window nonzero, unit-stride in ColMajor), replacing the 16 hand-derived strided-axpy address expressions across the uplo x layout x function grid; coo_rsksys reduces to it via the transpose identity (flip layout and uplo, transpose-view S, swap the window offsets). Beta scaling and the alpha == 0 short-circuit stay in the wrappers, applied exactly once. sketch_symmetric collapses from 12 overloads (4 call shapes x undefined primary + DenseSkOp + SparseSkOp) to 4: one constrained definition per shape with if-constexpr dispatch and a static_assert with a readable message on unsupported operator kinds. Previously an unsupported SketchingOperator type selected a declared-but-undefined primary template and died with an undefined-symbol link error. The FULL(S) overloads forward to the SUBMAT(S) definitions with zero offsets. lsksy3/rsksy3 share one transpose_copy_to_layout helper for the layout-mismatch path.

… SkOp path; reuse library utilities New negative tests assert RandBLAS::Error on the contracts the previous commits added: leading-dimension lower bounds (spsymm and the sparse sketch_symmetric branch), submatrix windows exceeding the operator, and one-based indices. The suite previously had zero error-path coverage for these entry points (the old symmetry-check test was removed with the API change and nothing replaced it). CaseD_Int32_Indices pins the two-SparseMatrix overload for a non-int64 index type against a hand-computed reference; it exercises whichever branch the build selects (expand-A + spgemm when the index width matches MKL_INT, densify-B otherwise). sparse_skop_unmaterialized hands sketch_symmetric an unfilled operator, exercising the submatrix_as_coo window-sampling path; the reference is built from a materialized twin with the same distribution and seed. Reuse cleanups: zero_other_triangle delegates to RandBLAS::overwrite_triangle instead of hand-rolled ColMajor loops, and the SkOp densification loop in test_sparse_skop is replaced by coo_to_dense (layout-aware, already tested).

make_dense_symmetric mirrors via RandBLAS::symmetrize, and the two hand-rolled two-pass dense-to-CSR converters (~55 lines) become overwrite_triangle + dense_to_csr calls. Behaviour unchanged: abs_tol=0 in dense_to_csr drops exactly the entries the manual loops skipped.

DevNotes and the dispatcher-adjacent RTD text still described the provisional design: an MKL path that fell back for CSC and side=Right (both are handled; side=Right never reaches MKL at all now), a shared internal::apply_beta_scale helper that never existed in the tree, and two-axpy-per-nonzero kernel descriptions. Rewritten to match the code: side normalization at dispatch, the validation and single-beta contract, column-driven parallel kernels, the Case D expand-A design with densify-B as the non-MKL fallback, and the window-sampling behaviour of the sparse sketch wrappers.

clang (Apple and LLVM alike, including the tsan lane) rejects referencing structured bindings inside OpenMP regions ('capturing a structured binding is not yet supported in OpenMP'); gcc accepts it, which is why the Linux gcc lanes were green while every clang lane failed to compile. Unpack the stride_64t fields into named int64_t locals before the parallel loops. No behavioural change.

Per review: the Uplo rework must extend the API, not break it. The four pre-Uplo sketch_symmetric overloads are restored verbatim in behaviour (runtime symmetry check via require_symmetric, both triangles read, forwarding to sketch_general) alongside the Uplo overloads. Overload resolution cannot collide: blas::Uplo is a scoped enum with no conversion to or from the legacy signatures' argument types. The old symmetry-check-throws test returns, plus a test pinning legacy results bitwise-adjacent to the Uplo overloads on symmetric input; the web docs advertise the legacy variants again under a marked dropdown.

Post-#196, left_spmm passes the caller's beta straight to MKL (which fuses it) and applies lascl only before the hand kernels; the old pre-scale-then-beta=1 pattern was removed there. spsymm now follows the same contract in both Cases: MKL and mkl_spgemm_to_dense receive beta (in Case D this also restores the alpha=1/beta=0 direct-write fast path, which the pre-applied beta=1 was defeating), and the dispatcher's lascl moved to just before the pure-accumulator fallbacks. Also per the #196 contract: empty products (a zero dimension, alpha == 0, or a structurally empty operand) now leave beta*Y at the dispatcher, since MKL rejects some valid empty sparse matrices at handle creation; without the guards, spsymm threw on inputs left_spmm handles gracefully. The public dispatch docstring keeps the contract; MKL routing mechanics moved to a plain comment and DevNotes.

…c register Window bounds now go through validate_submat_dims (the overflow-safe validator #196 centralized; the addition-form checks it replaced could overflow and accept negative dims), and the manual layout-flip ternaries go through flipped_layout, including one site where a local variable shadowed that helper's name. The spsymm kernels and the *_to_dense templates use the SignedInteger concept like every neighboring kernel. Documentation register fixes: internal helpers (the fallback kernels, lsksy3/rsksy3/lsksys/rsksys) drop their Doxygen-visible /// blocks for plain comments; the public spsymm convenience wrapper's math renders through the \math alias; all references to a non-repository planning file are gone (the four-case design lives in sparse_data/DevNotes.md); require_symmetric's docstring describes the current legacy-overload relationship instead of a phase history. Symmetric<SpMat> and as_symmetric gain their own directives in the web API reference. sksy.hh includes <vector> and <type_traits> directly.

…se D pairings Test-case names follow the snake_case behavior-phrase convention (STYLE_GUIDE 'Other files'). New coverage: empty sparse operands leave beta*Y in both spsymm overloads (the #196 contract the previous commit adopted); the three previously untested sparse-times-sparse format pairings (CSR-COO, COO-CSC, CSC-COO) and a non-CSR side=Right cell, completing the 3x3 grid the docs claim. The remaining std::mt19937 use carries its stated reason per CONTRIBUTING.

… helpers Every timed method is now checked against a trusted reference before its timing rows print (PASS/FAIL column), matching the practice in saso_sampling_performance. The program gains --help, argument validation, a --threads flag with requested-vs-effective reporting, and the OpenMPSettingsGuard from RandBLAS/testing/benchmarking.hh. Timing uses int64_t with static_cast, the file carries the full license block, and comments describe methods by behavior (one-triangle vs both-triangles vs GEMM-forwarding) rather than by PR chronology.

mmelnich · 2026-08-24T17:52:25Z

@rileyjmurray how's this looking?

rileyjmurray

There are lots of design questions here. These are just some initial thoughts.

rileyjmurray · 2026-08-25T14:42:19Z

-Symmetric matrices have to be stored as general matrices.
-  This stems partly from a desire for sketch_symmetric work equally well with DenseSkOp and SparseSkOp.
-  Another reason is that BLAS' SYMM function doesn't allow transposes, which is a key tool we use
-  in sketch_general to resolve layout discrepancies between the various arguments.
+sketch_symmetric supports both DenseSkOp and SparseSkOp.
+  ``sketch_symmetric`` now takes a ``blas::Uplo`` and exploits symmetry of :math:`A`:
+  the ``DenseSkOp`` branch dispatches to ``blas::symm`` (via ``lsksy3`` / ``rsksy3``)
+  and the ``SparseSkOp`` branch dispatches to a column-driven accumulation kernel
+  (``lsksys`` / ``rsksys``) that reads only the triangle of :math:`A` named by
+  ``uplo``. See ``RandBLAS/sparse_data/DevNotes.md`` for the access-pattern detail.
+
+Layout-mismatched ``DenseSkOp`` in ``sketch_symmetric`` transpose-copies the operand.
+  When the ``DenseSkOp``'s storage layout differs from the caller's ``layout`` parameter, ``sketch_symmetric`` transpose-copies the operand into a tight buffer in the caller's layout and then calls ``blas::symm``. The copy is ``O(d * n)``; ``blas::symm`` has no on-the-fly transpose flag for the dense operand, so this is what it costs to keep the SYMM speedup over a ``blas::gemm`` fallback. The layout-matched case skips the copy.


This material is under a Limitations heading, but not written in a way that describes limitations. It's just a very detailed behavior explanation of a certain codepath.Please resolve. Have the robot read this file in its entirety so it keeps style consistent.

Note: sketching operator layout-mismatches are already resolved by copy-transpose, and we don't mention that on this page; see fill_dense_unpacked.

rileyjmurray · 2026-08-25T14:45:56Z

+These overloads accept a ``blas::Uplo`` parameter naming the triangle of
+:math:`\mtxA` that is structurally stored; the opposite triangle is implied
+by symmetry and is **not** read. Both ``DenseSkOp`` and ``SparseSkOp`` are
+supported: DenseSkOp dispatches to ``blas::symm`` via ``lsksy3`` / ``rsksy3``,
+SparseSkOp dispatches to a column-driven accumulation kernel via
+``lsksys`` / ``rsksys``. See ``RandBLAS/sparse_data/DevNotes.md`` for the
+SparseSkOp access-pattern detail.


This level of detail isn't present in similar pages.

rileyjmurray · 2026-08-25T14:49:57Z

+       The "dense-symm A times sparse SkOp" case (Case B) is also implemented:
+       it lives in ``sketch_symmetric`` (the SparseSkOp branch) and uses a
+       column-driven accumulation kernel that reads only the named
+       triangle of A. See ``RandBLAS/sparse_data/DevNotes.md`` for the
+       access pattern.
+
+       The "sparse-symm A times sparse RHS, dense output" case (Case D) is
+       implemented via a two-``SparseMatrix``-arg ``spsymm`` overload,
+       covering all 3 x 3 = 9 sparse-format pairings for ``(A, B)``.
+       ``mkl_sparse_sp2m`` rejects symmetric descriptors and
+       ``mkl_sparse_?_spmmd`` takes no descriptor, so the symmetric
+       expansion happens on the RandBLAS side: on MKL builds the stored
+       triangle of ``A`` is expanded to a general sparse matrix in
+       ``O(nnz)`` memory and the existing sparse-times-sparse routine runs
+       with a GENERAL descriptor, keeping ``B`` sparse. Non-MKL builds
+       (and index widths mismatched with ``MKL_INT``) fall back to
+       densifying ``B`` into a temporary ``m``-by-``n`` buffer and
+       composing through the Case-C dispatcher.


References to Cases B, C, D aren't appropriate here, since those aren't defined on this page.

Separately --- comments at the kernel-dispatch level don't belong on the website.

rileyjmurray · 2026-08-25T14:51:48Z

+.. dropdown:: :math:`\mtxY = \alpha \cdot \mtxA \cdot \mtxB + \beta \cdot \mtxY,` with sparse symmetric :math:`\mtxA` (spsymm)
+    :animate: fade-in-slide-down
+    :color: light
+
+    .. doxygenfunction:: RandBLAS::spsymm(blas::Layout layout, blas::Uplo uplo, int64_t m, int64_t n, T alpha, const SpMat &A, const T *B, int64_t ldb, T beta, T *Y, int64_t ldy)
+      :project: RandBLAS
+
+    .. doxygenfunction:: RandBLAS::spsymm(blas::Layout layout, int64_t m, int64_t n, T alpha, const Symmetric<SpMat> &A_sym, const T *B, int64_t ldb, T beta, T *Y, int64_t ldy)
+      :project: RandBLAS


Variable name comment: name the matrices (A, B, C) or (A, X, Y). The current approach (A, B, Y) is no good. Apply appropriate changes everywhere that makes sense -- source and docs.

rileyjmurray · 2026-08-25T14:55:10Z

 .. dropdown:: :math:`\mtxB = \alpha \cdot \op(\submat(\mtxS))\cdot \op(\mtxA) + \beta \cdot \mtxB`
    :animate: fade-in-slide-down
    :color: light

    .. doxygenfunction:: RandBLAS::sketch_sparse(blas::Layout layout, blas::Op opS, blas::Op opA, int64_t d, int64_t n, int64_t m, T alpha, const DenseSkOp &S, int64_t S_ro, int64_t S_co, const SpMat &A, T beta, T *B, int64_t ldb) 
      :project: RandBLAS

 .. dropdown:: :math:`\mtxB = \alpha \cdot \op(\mtxA)\cdot \op(\submat(\mtxS)) + \beta \cdot \mtxB`
    :animate: fade-in-slide-down
    :color: light

    .. doxygenfunction:: RandBLAS::sketch_sparse(blas::Layout layout, blas::Op opA, blas::Op opS, int64_t m, int64_t d, int64_t n, T alpha, const SpMat &A, const DenseSkOp &S, int64_t S_ro, int64_t S_co, T beta, T *B, int64_t ldb) 
      :project: RandBLAS


Don't restrict template to DenseSkOp.

rileyjmurray · 2026-08-25T14:57:38Z

+    .. doxygenfunction:: RandBLAS::spsymm(blas::Layout layout, blas::Uplo uplo, int64_t m, int64_t n, T alpha, const SpMat &A, const T *B, int64_t ldb, T beta, T *Y, int64_t ldy)
+      :project: RandBLAS


Does MKL have two dimensional parameters (m and n) in its SPMM function? I think it'd only have one.

rileyjmurray · 2026-08-25T14:59:03Z

+    .. doxygenfunction:: RandBLAS::spsymm(blas::Layout layout, int64_t m, int64_t n, T alpha, const Symmetric<SpMat> &A_sym, const T *B, int64_t ldb, T beta, T *Y, int64_t ldy)
+      :project: RandBLAS


If we're going use function overloading here then the typed Symmetric<SpMat> argument should probably be in an overload of spmm, not spsymm.

Rename output matrix (A,B,Y) -> (A,B,C) to match spmm/spgemm. Slim public spsymm to a single dimension arg (order of A comes from the matrix, matching MKL's mkl_sparse_?_mm convention). Move the Symmetric<SpMat> carrier overload from spsymm to spmm. Trim FAQ/sketch_dense/sketch_sparse docs to drop pseudo-limitations and dispatch-level detail that don't belong on the website.

…hanges FAQ: correct the sparse-symmetric answer (spmm+Symmetric is a deterministic product, not a sketch); restore the layout-mismatch transpose-copy note as an actual Limitations entry; fix the stale sksy naming-convention line. sketch_sparse.rst: restore Case D documentation. symmetric.hh: fix a docstring that contradicted the still-live bare-SpMat spsymm entry point; add n_rows/n_cols aliases matching the DenseSkOp/ SparseSkOp convention. spmm(Symmetric) now delegates to spsymm instead of duplicating its logic. spsymm_dispatch.hh now includes spmm_dispatch.hh so any TU seeing the new overload sees the full spmm overload set. Rename stragglers (Ysk_ref, a stale side=Left comment) cleaned up.

rileyjmurray reviewed May 14, 2026

View reviewed changes

Comment thread RandBLAS/sparse_data/csr_spsymm_impl.hh Outdated

rileyjmurray requested changes May 15, 2026

View reviewed changes

mmelnich force-pushed the add-symm-kernels branch from 513568a to 0af25b0 Compare May 15, 2026 20:10

mmelnich force-pushed the add-symm-kernels branch from e7986dc to c8b036d Compare August 13, 2026 19:29

mmelnich mentioned this pull request Aug 13, 2026

Adopt one-triangle symmetric sketching (BallisticLA/RandBLAS#163) in ExplicitSymLinOp BallisticLA/RandLAPACK#167

Open

rileyjmurray reviewed Aug 13, 2026

View reviewed changes

Comment thread RandBLAS/sparse_data/coo_sksys_impl.hh

Comment thread RandBLAS/sparse_data/mkl_spmm_impl.hh

Comment thread RandBLAS/sparse_data/mkl_spmm_impl.hh

Comment thread rtd/source/api_reference/sketch_dense.rst

mmelnich added 15 commits August 24, 2026 12:58

mmelnich added 18 commits August 24, 2026 13:01

sparse_data: make the *_to_dense helpers index-generic

fc2af90

coo/csr/csc_to_dense were bound to the default int64_t index type, so any code routing a non-int64-indexed matrix through them failed to compile. Deduce sint_t (default preserved); no caller changes needed.

mmelnich force-pushed the add-symm-kernels branch from e77e7c3 to bd6ea1a Compare August 24, 2026 17:36

rileyjmurray requested changes Aug 25, 2026

View reviewed changes

mmelnich added 2 commits August 26, 2026 10:52

		.. doxygenfunction:: RandBLAS::spsymm(blas::Layout layout, blas::Uplo uplo, int64_t m, int64_t n, T alpha, const SpMat &A, const T B, int64_t ldb, T beta, T Y, int64_t ldy)
		:project: RandBLAS

Uh oh!

Conversation

mmelnich commented May 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

The four cases

What's in the PR

API: extension, not a break

Test coverage

Uh oh!

mmelnich commented May 13, 2026

Uh oh!

Uh oh!

rileyjmurray commented May 14, 2026

Uh oh!

mmelnich commented May 14, 2026

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

mmelnich commented Aug 13, 2026

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

mmelnich commented Aug 23, 2026

Uh oh!

rileyjmurray commented Aug 24, 2026

Uh oh!

mmelnich commented Aug 24, 2026

Uh oh!

rileyjmurray left a comment

Choose a reason for hiding this comment

Uh oh!

rileyjmurray Aug 25, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

rileyjmurray Aug 25, 2026

Choose a reason for hiding this comment

Uh oh!

rileyjmurray Aug 25, 2026

Choose a reason for hiding this comment

Uh oh!

rileyjmurray Aug 25, 2026

Choose a reason for hiding this comment

Uh oh!

rileyjmurray Aug 25, 2026

Choose a reason for hiding this comment

Uh oh!

rileyjmurray Aug 25, 2026

Choose a reason for hiding this comment

Uh oh!

rileyjmurray Aug 25, 2026

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

mmelnich commented May 13, 2026 •

edited

Loading