Add SYMM kernels for sparse symmetric matrices; rework sketch_symmetric to exploit symmetry#163
Open
mmelnich wants to merge 19 commits into
Open
Add SYMM kernels for sparse symmetric matrices; rework sketch_symmetric to exploit symmetry#163mmelnich wants to merge 19 commits into
mmelnich wants to merge 19 commits into
Conversation
Contributor
Author
|
@rileyjmurray, this is still a WIP, but feel free to look. |
Contributor
|
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 See also RandBLAS/RandBLAS/sparse_data/mkl_spmm_impl.hh Lines 379 to 383 in d8a153f |
mmelnich
added a commit
that referenced
this pull request
May 14, 2026
…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).
Contributor
Author
|
@rileyjmurray lmk your thoughts |
rileyjmurray
requested changes
May 15, 2026
Comment on lines
+137
to
+138
| prevents a symmetric sparse matrix from being accidentally passed to | ||
| the general ``spmm`` / ``spgemm`` routines. |
Contributor
There was a problem hiding this comment.
the "symmetric sparse matrix" mentioned here means "a matrix that the user interprets as symmetric while only storing one of the two triangles", right?
Contributor
Author
There was a problem hiding this comment.
Yes, correct. That would be the only dangerous case for spmm/spgemm application path.
| | A | dense-symm × dense | dense, one triangle | dense | Implemented via ``blas::symm`` in ``sksy.hh``. | | ||
| | B | dense-symm × sparse | dense, one triangle | sparse | Implemented via hand-rolled ``lsksys`` / ``rsksys`` in ``sksy.hh`` (two-axpy scatter per stored nonzero of S, reading only the named triangle of A). | | ||
| | C | sparse-symm × dense (→ dense) | sparse, one triangle | dense | Implemented in ``spsymm_dispatch.hh`` (MKL fast path + per-format fallbacks). | | ||
| | D | sparse-symm × sparse → dense | sparse, one triangle | sparse | Implemented via densify-B + Case-C composition in ``spsymm_dispatch.hh``. | |
Contributor
There was a problem hiding this comment.
This would be the only place where RandBLAS supports sparse-times-sparse. That kind of feature needs to be introduced in a separate PR.
Contributor
Author
There was a problem hiding this comment.
But we already have sparse * sparse through MKL.
Comment on lines
+417
to
+419
| /// * Defines :math:`\submat(\mtxS).` SparseSkOp is currently not supported and | ||
| /// calls with a SparseSkOp will throw a :math:`\ttt{RandBLAS::Error}` --- | ||
| /// this is the Case-B kernel from `project-plans/randblas-symm-plan.md`. |
Contributor
There was a problem hiding this comment.
So this actually reduces RandBLAS' functionality. Before this PR we'd be able to sketch a dense structurally-symmetric matrix with a SparseSkOp. After this PR we'd throw an error. That's a "move to RandBLAS 2.0" kind of change. I don't want to make that move in this PR. At the very least this PR has a bunch of stuff that's perfectly useful for RandBLAS 1.x.
Contributor
Author
There was a problem hiding this comment.
Stale docstring, sorry. Case B is actually implemented in this PR (commit 77a3b45): SparseSkOp dispatches to a hand-rolled lsksys / rsksys two-AXPY scatter. No throw, no functionality removed. Fixing the docstring in the next commit.
Comment on lines
+475
to
486
| template <typename T, typename RNG, RandBLAS::SignedInteger sint_t> | ||
| inline void sketch_symmetric( | ||
| blas::Layout layout, blas::Uplo uplo, | ||
| int64_t d, int64_t n, | ||
| T alpha, | ||
| const SparseSkOp<T, RNG, sint_t> &S, int64_t ro_s, int64_t co_s, | ||
| const T* A, int64_t lda, | ||
| T beta, | ||
| T* B, int64_t ldb | ||
| ) { | ||
| RandBLAS::sparse::lsksys(layout, uplo, d, n, alpha, S, ro_s, co_s, A, lda, beta, B, ldb); | ||
| } |
Contributor
There was a problem hiding this comment.
If this PR changes things so that LSKSYS errors on a SparseSkOp S, then we should just have a static_assert(false) here.
Contributor
Author
There was a problem hiding this comment.
This isn't a stub. Case B is implemented in this PR (commit 77a3b45). lsksys is a real kernel doing the per-stored-nonzero two-AXPY scatter against dense symm A; the wrapper just forwards into it. After the refactor on this branch it'll call into RandBLAS::sparse_data::coo_lsksys instead, but it's still a real kernel with no error path.
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.
mmelnich
force-pushed
the
add-symm-kernels
branch
from
513568a to
0af25b0
Compare
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.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Adds SYMM-shaped kernels to RandBLAS so that symmetric matrices can finally be sketched with their symmetry exploited. Implements all four cases (A, B, C, D), and folds in repo-wide cleanups uncovered along the way (a
lasclutility, helper extractions, a CSC fast path).Drives the funNyström++ work on
BallisticLA/RandLAPACK#132:sketch_symmetricwas silently forwarding tosketch_general→ GEMM, so A's symmetry wasn't exploited and a sparse SkOp had to be densified before the call. This PR closes the dense side and lays the foundation for the sparse-symmetric side.Design plan:
randblas-symm-plan.md. Condensed 4-case table + MKL availability matrix lives inRandBLAS/sparse_data/DevNotes.md.The four cases
blas::symminsksy.hhlsksys/rsksyswrappers insksy.hhhandle the SparseSkOp materialization-and-recurse; the COO walk plus two-axpy-per-stored-nonzero scatter lives inRandBLAS::sparse_data::coo_lsksys/coo_rsksysinsparse_data/coo_sksys_impl.hhand reads only A's named trianglesparse_data/spsymm_dispatch.hh(MKL fast path covers all of{CSR, CSC, COO} × {side=Left, side=Right}; per-format hand kernels remain as the non-MKL fallback)sparse_data/spsymm_dispatch.hh(newspsymmoverload taking twoSparseMatrixargs)What's in the PR
Cases A, B, C, all real implementations.
sketch_symmetricdispatches by SkOp type. DenseSkOp goes throughdense::lsksy3/rsksy3(SYMM-backed, with transpose-copy into matching layout on layout-mismatch so the SYMM speedup is retained). SparseSkOp goes throughsparse::lsksys/rsksyswrappers insksy.hh, which handle the SparseSkOp materialization-and-recurse and forward intoRandBLAS::sparse_data::coo_lsksys/coo_rsksysinsparse_data/coo_sksys_impl.hh; the kernel does a per-stored-nonzero two-AXPY scatter that reads only A's named triangle, with submatrix offsets filtered inline on the COO view.spsymmdispatches to an MKL fast path covering both Side values: side=Left direct, side=Right via opposite-layout reinterpretation of B/Y (no-copy trick exploitingA == A^T), and CSC viaA.transpose()lightweight CSR view plus flippeduplo. The hand-rolled per-format fallback (csr/csc/coo_spsymm_impl.hh) remains as the non-MKL path; its per-stored-entry scatter is shared viainternal::spsymm_scatter_{left,right}helpers.Case D, densify-B plus Case-C composition. The new two-
SparseMatrix-argspsymmoverload allocates anm-by-ntemporarystd::vector<T>in the caller's layout, fills it via the format-specificcoo_to_dense/csr_to_dense/csc_to_densehelper picked byif constexpr, and then calls the existing Case-C dispatcher on the densified buffer. Covers all 3 × 3 = 9 sparse-format pairings for(A, B). MKL'smkl_sparse_sp2mreturnsSPARSE_STATUS_NOT_SUPPORTEDwhendescrA.type == SPARSE_MATRIX_TYPE_SYMMETRIC(onlyGENERALis accepted there), andmkl_sparse_d_spmmdtakes no descriptor at all, so the symmetric expansion has to happen on the RandBLAS side either way. Cost:O(m * n)temporary buffer, small for the typical RandNLA workload whereBis a sketching operator withnnz(B) << m*n.Symmetric<SpMat>wrapper (sparse_data/symmetric.hh): non-owning carrier holdingconst SpMat&andblas::Uplo. Square-matrix invariant at construction. Type-safety guard: aSymmetric<SpMat>argument fails to bind to generalspmm/spgemm. PublicRandBLAS::spsymmhas both a raw (SpMat+uplo) overload and a wrapper-taking overload.Cleanups along the way:
RandBLAS::util::lascl(util.hh):A := α·Afor a dense m×n matrix, LAPACK-style naming. Replaces hand-rolled "loop-over-blas::scal" beta-scale patterns at 7 sites acrosstrsm_dispatch.hh,spmm_dispatch.hh,mkl_spmm_impl.hh(mkl_spgemm_to_dense), the three new*_spsymm_impl.hhfiles, and two example programs.internal::mkl_sparse_mm_call(mkl_spmm_impl.hh): type-dispatched wrapper aroundmkl_sparse_d_mm/mkl_sparse_s_mm, shared bymkl_left_spmm(GENERAL descriptor) andmkl_spsymm(SYMMETRIC descriptor).<numeric>include insparse_data/base.hh: latent gap (usedstd::iotabut had relied on transitive inclusion viacoo_matrix.hh); the newsymmetric.hhexposed it.RandBLAS.hhumbrella: now pulls insparse_data/spsymm_dispatch.hhso downstream#include <RandBLAS.hh>consumers seeRandBLAS::spsymmand theSymmetric<SpMat>wrapper automatically (mirrors howspmm/spgemm/sketch_symmetricare already exposed).Perf benchmark: new
examples/simple-kernel-benchmarks/spsymm_performance.cccompares (a) sparse:RandBLAS::spsymm(one-triangle storage + MKL fast path) vs the pre-PR workaround (both triangles +RandBLAS::spmm) vs a denseblas::symmreference, and (b) dense: the newsketch_symmetric(SYMM-backed) vs the equivalentsketch_generalcall (the pre-PR GEMM-forwarding behaviour). Single-config CLI or default 3-point sweep.API break worth flagging
sketch_symmetricnow takesblas::Uplo uplo(previously:T sym_check_tol = 0). All four overloads were affected. Justification:sketch_symmetricperformed an O(n²) runtime symmetry check viarequire_symmetric, then forwarded to GEMM. The check was a workaround for the fact that GEMM reads both triangles.blas::symm, only the named triangle is read. The check is meaningless and pure overhead.uplois the BLAS-standard way to name which triangle is read.In-tree callers patched in this PR:
test/linops/test_sketch_symmetric.cc. The only known external consumer is RandLAPACK'sExplicitSymLinOp(via the funNyström++ branch); follow-up RandLAPACK PR will threaduplothrough.require_symmetricremains available as a standalone validator.What's still missing
Deliberate scope cuts; not blockers for this PR.
B, but the dispatch grid is messier (mixed A/B sparse formats × side × layout × uplo) and no portable reference exists in MKL / Ginkgo / SparseBLAS reference / MAGMA-sparse (surveyed 2026-05). The densify-B + Case-C composition shipped here is correct, covers all 9 format pairings, and has a small constant-factor overhead in the regime that matters for RandNLA workloads. If the dense-temp cost becomes a bottleneck, a hand-rolled kernel can be a targeted follow-up.ExplicitSymLinOpand any othersketch_symmetricconsumers needuplothreaded through. Small follow-up; will open after this PR merges.Test coverage
24-cell
{COO, CSR, CSC} × {ColMajor, RowMajor} × {Upper, Lower} × {Left, Right}forspsymmCase C (intest/linops/test_spsymm.cc), plus float/beta=0/alpha=0/wrapper-routing edge cases, plus 7 Case D tests ({CSR-CSR, CSC-CSC, COO-COO, Mixed-format} × Left, plusCSR-CSR-Right,Float, andAlphaZero-BetaScale), plus 5 cases for theSymmetric<SpMat>wrapper, plus 14 cases intest_sketch_symmetric.cc(8 DenseSkOp paths + 6 new SparseSkOp paths exercising Case B across both sides, both layouts, both uplo, and a lift configuration).Tolerance for
spsymmtests:atol = 100*eps,rtol = 10*eps. The dense reference (blas::symmon a fully-symmetrized A) and the sparse path accumulate FMAs in a different order (off-diagonal entries contribute via two AXPYs in the sparse path vs. one fused dotprod in dense SYMM), giving a few-ULP divergence. The relaxed tolerance documents this is expected, not a bug.Verified locally: 477/477 ctest pass on Linux + GCC 13.3 + CUDA-aware blaspp + MKL sparse (includes 17/17 TestSpsymm = 10 Case C + 7 Case D, plus 14/14 TestSketchSymmetric = 8 DenseSkOp + 6 SparseSkOp).