From 168198e0a6502b38832800827bb9fc656330fde1 Mon Sep 17 00:00:00 2001 From: mmelnich Date: Wed, 13 May 2026 18:03:51 -0400 Subject: [PATCH 01/19] Add SYMM-backed dispatch to sketch_symmetric; stub sparse-SkOp branch 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. --- RandBLAS/sksy.hh | 746 +++++++++++++-------------- RandBLAS/util.hh | 6 +- test/linops/test_sketch_symmetric.cc | 63 +-- 3 files changed, 388 insertions(+), 427 deletions(-) diff --git a/RandBLAS/sksy.hh b/RandBLAS/sksy.hh index e49208b6..4d26da13 100644 --- a/RandBLAS/sksy.hh +++ b/RandBLAS/sksy.hh @@ -33,161 +33,167 @@ #include "RandBLAS/base.hh" #include "RandBLAS/skge.hh" -namespace RandBLAS { -using namespace RandBLAS::dense; -using namespace RandBLAS::sparse; +// ============================================================================= +// Symmetric sketching helpers (SYMM-backed). See project-plans/randblas-symm-plan.md +// for the four-case API design and the implementation status of each case. +// - lsksy3, rsksy3: Case A (dense-symm A x dense Omega), via blas::symm. +// - SparseSkOp branches of sketch_symmetric: Case B stubs (not implemented). +// ============================================================================= +namespace RandBLAS::dense { + +// ============================================================================= +/// LSKSY3: SYMM-backed left-sketch with a symmetric matrix A. +/// +/// Computes B = alpha * submat(S) * mat(A) + beta * B, where: +/// - mat(A) is n-by-n symmetric. Only the triangle named by `uplo` is read. +/// - submat(S) is the d-by-n view of S at (ro_s, co_s). +/// - mat(B) is d-by-n. +/// +/// When S has no materialized buffer, the submatrix is realized via +/// `submatrix_as_blackbox` (same pattern as `lskge3`). When the buffered S's +/// storage layout matches the caller's `layout`, the final call is +/// `blas::symm` with `side = Right` (since A is on the right of S in the +/// operation). When layouts mismatch, SYMM cannot transpose S on-the-fly; +/// this first cut falls back to `blas::gemm` with `opS = Trans`, sacrificing +/// the SYMM speedup on that path. Optimization target: transpose-copy of S +/// into matching layout, then SYMM. See project-plans/randblas-symm-plan.md §7. +template +void lsksy3( + blas::Layout layout, + blas::Uplo uplo, + int64_t d, // B is d-by-n + int64_t n, // A is n-by-n, S is d-by-n (after submat view) + T alpha, + const DenseSkOp &S, + int64_t ro_s, + int64_t co_s, + const T *A, + int64_t lda, + T beta, + T *B, + int64_t ldb +){ + constexpr bool maybe_denseskop = !std::is_same_v, BLASFriendlyOperator>; + if constexpr (maybe_denseskop) { + if (!S.buff) { + auto submat_S = submatrix_as_blackbox>(S, d, n, ro_s, co_s); + lsksy3(layout, uplo, d, n, alpha, submat_S, 0, 0, A, lda, beta, B, ldb); + return; + } + } + randblas_require( S.buff != nullptr ); + randblas_require( S.n_rows >= d + ro_s ); + randblas_require( S.n_cols >= n + co_s ); + if (layout == blas::Layout::ColMajor) { + randblas_require(lda >= n); + randblas_require(ldb >= d); + } else { + randblas_require(lda >= n); + randblas_require(ldb >= n); + } + + auto [pos, lds] = offset_and_ldim(S.layout, S.n_rows, S.n_cols, ro_s, co_s); + T* S_ptr = &S.buff[pos]; + + if (S.layout == layout) { + // Fast path: SYMM directly. + blas::symm(layout, blas::Side::Right, uplo, d, n, + alpha, A, lda, S_ptr, lds, beta, B, ldb); + } else { + // Layout mismatch: SYMM has no opS flag. Fall back to GEMM with + // opS = Trans. Loses the SYMM 1.3-1.8x speedup on this path but is + // correct. Optimization target tracked in randblas-symm-plan.md §7. + blas::gemm(layout, blas::Op::Trans, blas::Op::NoTrans, d, n, n, + alpha, S_ptr, lds, A, lda, beta, B, ldb); + } + return; +} -// MARK: SUBMAT(S) // ============================================================================= -/// \fn sketch_symmetric(blas::Layout layout, int64_t n, -/// int64_t d, T alpha, const T *A, int64_t lda, -/// const SKOP &S, int64_t ro_s, int64_t co_s, -/// T beta, T *B, int64_t ldb, T sym_check_tol = 0 -/// ) -/// @verbatim embed:rst:leading-slashes -/// Check that :math:`\mat(A)` is symmetric up to tolerance :math:`\texttt{sym_check_tol}`, then sketch from the right in a SYMM-like operation -/// -/// .. math:: -/// \mat(B) = \alpha \cdot \underbrace{\mat(A)}_{n \times n} \cdot \underbrace{\submat(\mtxS)}_{n \times d} + \beta \cdot \underbrace{\mat(B)}_{n \times d}, \tag{$\star$} -/// -/// where :math:`\alpha` and :math:`\beta` are real scalars and :math:`\mtxS` is a sketching operator. -/// -/// .. dropdown:: FAQ -/// :animate: fade-in-slide-down -/// -/// **What's** :math:`\mat(A)?` -/// -/// It's a symmetric matrix of order :math:`n`. Its precise contents depend on :math:`(A, \lda)`, -/// according to -/// -/// .. math:: -/// \mat(A)_{ij} = A[i + j \cdot \lda] = A[i \cdot \lda + j]. -/// -/// Note that the the "layout" parameter passed to this function is not used here. -/// That's because this function requires :math:`\mat(A)` to be stored in the format -/// of a general matrix (with both upper and lower triangles). -/// -/// This function's default behavior is to check that :math:`\mat(A)` is symmetric before -/// attempting sketching. That check can be skipped (at your own peril!) by calling this -/// function with sym_check_tol < 0. -/// -/// **What's** :math:`\mat(B)?` -/// -/// It's an :math:`n \times d` matrix. Its precise contents depend on :math:`(B,\ldb)` and "layout." -/// -/// If layout == ColMajor, then -/// -/// .. math:: -/// \mat(B)_{ij} = B[i + j \cdot \ldb]. -/// -/// In this case, :math:`\ldb` must be :math:`\geq n.` -/// -/// If layout == RowMajor, then -/// -/// .. math:: -/// \mat(B)_{ij} = B[i \cdot \ldb + j]. -/// -/// In this case, :math:`\ldb` must be :math:`\geq d.` -/// -/// **What is** :math:`\submat(\mtxS)` **?** -/// -/// It's the :math:`n \times d` submatrix of :math:`{\mtxS}` whose upper-left corner appears -/// at index :math:`(\texttt{ro_s}, \texttt{co_s})` of :math:`{\mtxS}.` -/// -/// .. dropdown:: Full parameter descriptions -/// :animate: fade-in-slide-down -/// -/// layout - [in] -/// * Either Layout::ColMajor or Layout::RowMajor -/// * Matrix storage for :math:`\mat(B).` -/// -/// n - [in] -/// * A nonnegative integer. -/// * The number of rows in :math:`\mat(B).` -/// * The number of rows and columns in :math:`\mat(A).` -/// -/// d - [in] -/// * A nonnegative integer. -/// * The number of columns in :math:`\mat(B)` and :math:`\submat(\mtxS).` -/// -/// alpha - [in] -/// * A real scalar. -/// * If zero, then :math:`A` is not accessed. -/// -/// A - [in] -/// * Pointer to a 1D array of real scalars. -/// * Defines :math:`\mat(A).` -/// -/// lda - [in] -/// * A nonnegative integer. -/// * Leading dimension of :math:`\mat(A)` when reading from :math:`A.` -/// -/// S - [in] -/// * A DenseSkOp or SparseSkOp object. -/// * Defines :math:`\submat(\mtxS).` -/// -/// ro_s - [in] -/// * A nonnegative integer. -/// * The rows of :math:`\submat(\mtxS)` are a contiguous subset of rows of :math:`S.` -/// * The rows of :math:`\submat(\mtxS)` start at :math:`S[\texttt{ro_s}, :].` -/// -/// co_s - [in] -/// * A nonnnegative integer. -/// * The columns of :math:`\submat(\mtxS)` are a contiguous subset of columns of :math:`S.` -/// * The columns of :math:`\submat(\mtxS)` start at :math:`S[:,\texttt{co_s}].` -/// -/// beta - [in] -/// * A real scalar. -/// * If zero, then :math:`B` need not be set on input. -/// -/// B - [in,out] -/// * Pointer to 1D array of real scalars. -/// * On entry, defines :math:`\mat(B)` -/// on the RIGHT-hand side of :math:`(\star).` -/// * On exit, defines :math:`\mat(B)` -/// on the LEFT-hand side of :math:`(\star).` +/// RSKSY3: SYMM-backed right-sketch with a symmetric matrix A. /// -/// ldb - [in] -/// * A nonnegative integer. -/// * Leading dimension of :math:`\mat(B)` when reading from :math:`B.` +/// Computes B = alpha * mat(A) * submat(S) + beta * B, where: +/// - mat(A) is n-by-n symmetric. Only the triangle named by `uplo` is read. +/// - submat(S) is the n-by-d view of S at (ro_s, co_s). +/// - mat(B) is n-by-d. /// -/// @endverbatim -template -inline void sketch_symmetric( - // B = alpha*A*S + beta*B, where A is a symmetric matrix stored in the format of a general matrix. +/// Same materialization and layout-mismatch fallback semantics as `lsksy3`. +/// Final call (matching layout): `blas::symm` with `side = Left`. +template +void rsksy3( blas::Layout layout, - int64_t n, // number of rows in B - int64_t d, // number of columns in B + blas::Uplo uplo, + int64_t n, // A is n-by-n, S is n-by-d (after submat view), B is n-by-d + int64_t d, T alpha, - const T* A, + const T *A, int64_t lda, - const SKOP &S, + const DenseSkOp &S, int64_t ro_s, int64_t co_s, T beta, - T* B, - int64_t ldb, - T sym_check_tol = 0 -) { - RandBLAS::util::require_symmetric(layout, A, n, lda, sym_check_tol); - sketch_general(layout, blas::Op::NoTrans, blas::Op::NoTrans, n, d, n, alpha, A, lda, S, ro_s, co_s, beta, B, ldb); + T *B, + int64_t ldb +){ + constexpr bool maybe_denseskop = !std::is_same_v, BLASFriendlyOperator>; + if constexpr (maybe_denseskop) { + if (!S.buff) { + auto submat_S = submatrix_as_blackbox>(S, n, d, ro_s, co_s); + rsksy3(layout, uplo, n, d, alpha, A, lda, submat_S, 0, 0, beta, B, ldb); + return; + } + } + randblas_require( S.buff != nullptr ); + randblas_require( S.n_rows >= n + ro_s ); + randblas_require( S.n_cols >= d + co_s ); + if (layout == blas::Layout::ColMajor) { + randblas_require(lda >= n); + randblas_require(ldb >= n); + } else { + randblas_require(lda >= n); + randblas_require(ldb >= d); + } + + auto [pos, lds] = offset_and_ldim(S.layout, S.n_rows, S.n_cols, ro_s, co_s); + T* S_ptr = &S.buff[pos]; + + if (S.layout == layout) { + blas::symm(layout, blas::Side::Left, uplo, n, d, + alpha, A, lda, S_ptr, lds, beta, B, ldb); + } else { + // Layout-mismatch fallback (see lsksy3). + blas::gemm(layout, blas::Op::NoTrans, blas::Op::Trans, n, d, n, + alpha, A, lda, S_ptr, lds, beta, B, ldb); + } + return; } +} // end namespace RandBLAS::dense + + +namespace RandBLAS { + +using namespace RandBLAS::dense; +using namespace RandBLAS::sparse; + + +// MARK: SUBMAT(S) // ============================================================================= -/// \fn sketch_symmetric(blas::Layout layout, int64_t d, -/// int64_t n, T alpha, const SKOP &S, int64_t ro_s, int64_t co_s, -/// const T *A, int64_t lda, T beta, T *B, int64_t ldb, T sym_check_tol = 0 -/// ) +/// \fn sketch_symmetric(blas::Layout layout, blas::Uplo uplo, +/// int64_t d, int64_t n, T alpha, +/// const SKOP &S, int64_t ro_s, int64_t co_s, +/// const T *A, int64_t lda, T beta, T *B, int64_t ldb +/// ) /// @verbatim embed:rst:leading-slashes -/// Check that :math:`\mat(A)` is symmetric up to tolerance :math:`\texttt{sym_check_tol}`, then sketch from the left in a SYMM-like operation -/// +/// Sketch from the left in a SYMM-like operation +/// /// .. math:: /// \mat(B) = \alpha \cdot \underbrace{\submat(\mtxS)}_{d \times n} \cdot \underbrace{\mat(A)}_{n \times n} + \beta \cdot \underbrace{\mat(B)}_{d \times n}, \tag{$\star$} -/// +/// /// where :math:`\alpha` and :math:`\beta` are real scalars and :math:`\mtxS` is a sketching operator. /// /// .. dropdown:: FAQ @@ -195,49 +201,28 @@ inline void sketch_symmetric( /// /// **What's** :math:`\mat(A)?` /// -/// It's a symmetric matrix of order :math:`n`. Its precise contents depend on :math:`(A, \lda)`, -/// according to -/// -/// .. math:: -/// \mat(A)_{ij} = A[i + j \cdot \lda] = A[i \cdot \lda + j]. -/// -/// Note that the the "layout" parameter passed to this function is not used here. -/// That's because this function requires :math:`\mat(A)` to be stored in the format -/// of a general matrix (with both upper and lower triangles). -/// -/// This function's default behavior is to check that :math:`\mat(A)` is symmetric before -/// attempting sketching. That check can be skipped (at your own peril!) by calling this -/// function with sym_check_tol < 0. -/// -/// **What's** :math:`\mat(B)?` +/// It's a symmetric matrix of order :math:`n`, stored with only the triangle named by +/// :math:`\texttt{uplo}` populated. The opposite triangle is not read. /// -/// It's a :math:`d \times n` matrix. Its precise contents depend on :math:`(B,\ldb)` and "layout." -/// -/// If layout == ColMajor, then -/// -/// .. math:: -/// \mat(B)_{ij} = B[i + j \cdot \ldb]. -/// -/// In this case, :math:`\ldb` must be :math:`\geq d.` -/// -/// If layout == RowMajor, then +/// The :math:`\texttt{layout}` parameter governs the indexing convention into :math:`A`: /// /// .. math:: -/// \mat(B)_{ij} = B[i \cdot \ldb + j]. +/// \mat(A)_{ij} = A[i + j \cdot \lda] \quad (\text{ColMajor}) +/// = A[i \cdot \lda + j] \quad (\text{RowMajor}). /// -/// In this case, :math:`\ldb` must be :math:`\geq n.` -/// -/// **What is** :math:`\submat(\mtxS)` **?** -/// -/// It's the :math:`d \times n` submatrix of :math:`{\mtxS}` whose upper-left corner appears -/// at index :math:`(\texttt{ro_s}, \texttt{co_s})` of :math:`{\mtxS}.` +/// Unlike the pre-SYMM API, both triangles need not match: only the stored triangle is +/// read by :math:`\ttt{blas::symm}` underneath. /// /// .. dropdown:: Full parameter descriptions /// :animate: fade-in-slide-down /// /// layout - [in] /// * Either Layout::ColMajor or Layout::RowMajor -/// * Matrix storage for :math:`\mat(B).` +/// * Matrix storage for :math:`\mat(A)` and :math:`\mat(B).` +/// +/// uplo - [in] +/// * Either Uplo::Upper or Uplo::Lower +/// * Names the triangle of :math:`\mat(A)` that is stored and read. /// /// d - [in] /// * A nonnegative integer. @@ -245,26 +230,26 @@ inline void sketch_symmetric( /// /// n - [in] /// * A nonnegative integer. -/// * The number of columns in :math:`\mat(B).` +/// * The number of columns in :math:`\mat(B)` and :math:`\submat(\mtxS).` /// * The number of rows and columns in :math:`\mat(A).` /// /// alpha - [in] /// * A real scalar. /// * If zero, then :math:`A` is not accessed. /// -/// S - [in] -/// * A DenseSkOp or SparseSkOp object. -/// * Defines :math:`\submat(\mtxS).` +/// S - [in] +/// * A SketchingOperator object (DenseSkOp or SparseSkOp). +/// * 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`. /// /// ro_s - [in] /// * A nonnegative integer. -/// * The rows of :math:`\submat(\mtxS)` are a contiguous subset of rows of :math:`S.` -/// * The rows of :math:`\submat(\mtxS)` start at :math:`S[\texttt{ro_s}, :].` +/// * The rows of :math:`\submat(\mtxS)` start at :math:`\mtxS[\texttt{ro_s}, :].` /// /// co_s - [in] -/// * A nonnnegative integer. -/// * The columns of :math:`\submat(\mtxS)` are a contiguous subset of columns of :math:`S.` -/// * The columns of :math:`\submat(\mtxS)` start at :math:`S[:,\texttt{co_s}].` +/// * A nonnegative integer. +/// * The columns of :math:`\submat(\mtxS)` start at :math:`\mtxS[:, \texttt{co_s}].` /// /// A - [in] /// * Pointer to a 1D array of real scalars. @@ -278,255 +263,254 @@ inline void sketch_symmetric( /// * A real scalar. /// * If zero, then :math:`B` need not be set on input. /// -/// B - [in,out] +/// B - [in, out] /// * Pointer to 1D array of real scalars. -/// * On entry, defines :math:`\mat(B)` -/// on the RIGHT-hand side of :math:`(\star).` -/// * On exit, defines :math:`\mat(B)` -/// on the LEFT-hand side of :math:`(\star).` +/// * On entry, defines :math:`\mat(B)` on the RIGHT-hand side of :math:`(\star).` +/// * On exit, defines :math:`\mat(B)` on the LEFT-hand side of :math:`(\star).` /// /// ldb - [in] /// * A nonnegative integer. /// * Leading dimension of :math:`\mat(B)` when reading from :math:`B.` /// /// @endverbatim -template +template inline void sketch_symmetric( - // B = alpha*S*A + beta*B - blas::Layout layout, - int64_t d, // number of rows in B - int64_t n, // number of columns in B + blas::Layout layout, blas::Uplo uplo, + int64_t d, int64_t n, T alpha, - const SKOP &S, - int64_t ro_s, - int64_t co_s, - const T* A, - int64_t lda, + const SKOP &S, int64_t ro_s, int64_t co_s, + const T* A, int64_t lda, T beta, - T* B, - int64_t ldb, - T sym_check_tol = 0 + T* B, int64_t ldb +); + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + int64_t d, int64_t n, + T alpha, + const DenseSkOp &S, int64_t ro_s, int64_t co_s, + const T* A, int64_t lda, + T beta, + T* B, int64_t ldb ) { - RandBLAS::util::require_symmetric(layout, A, n, lda, sym_check_tol); - sketch_general(layout, blas::Op::NoTrans, blas::Op::NoTrans, d, n, n, alpha, S, ro_s, co_s, A, lda, beta, B, ldb); + RandBLAS::dense::lsksy3(layout, uplo, d, n, alpha, S, ro_s, co_s, A, lda, beta, B, ldb); +} + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + int64_t d, int64_t n, + T alpha, + const SparseSkOp &S, int64_t ro_s, int64_t co_s, + const T* A, int64_t lda, + T beta, + T* B, int64_t ldb +) { + (void) layout; (void) uplo; (void) d; (void) n; (void) alpha; + (void) S; (void) ro_s; (void) co_s; (void) A; (void) lda; + (void) beta; (void) B; (void) ldb; + randblas_require( + false && + "RandBLAS::sketch_symmetric with a sparse sketching operator is the " + "Case-B kernel from project-plans/randblas-symm-plan.md. Not implemented " + "in this PR --- the API signature is reserved so future PRs can fill in " + "the body without breaking source compatibility. Composition fallback: " + "densify the SkOp then call sketch_symmetric on the dense version." + ); } -// MARK: FULL(S) // ============================================================================= -/// \fn sketch_symmetric(blas::Layout layout, T alpha, -/// const T *A, int64_t lda, const SKOP &S, -/// T beta, T *B, int64_t ldb, T sym_check_tol = 0 -/// ) +/// \fn sketch_symmetric(blas::Layout layout, blas::Uplo uplo, +/// int64_t n, int64_t d, T alpha, +/// const T *A, int64_t lda, +/// const SKOP &S, int64_t ro_s, int64_t co_s, +/// T beta, T *B, int64_t ldb +/// ) /// @verbatim embed:rst:leading-slashes -/// Check that :math:`\mat(A)` is symmetric up to tolerance :math:`\texttt{sym_check_tol}`, then sketch from the right in a SYMM-like operation -/// -/// .. math:: -/// \mat(B) = \alpha \cdot \underbrace{\mat(A)}_{n \times n} \cdot \mtxS + \beta \cdot \underbrace{\mat(B)}_{n \times d}, \tag{$\star$} -/// -/// where :math:`\alpha` and :math:`\beta` are real scalars and :math:`\mtxS` is an :math:`n \times d` sketching operator. -/// -/// .. dropdown:: FAQ -/// :animate: fade-in-slide-down -/// -/// **What's** :math:`\mat(A)?` -/// -/// It's a symmetric matrix of order :math:`n`, where :math:`n = \texttt{S.dist.n_cols}`. -/// Its precise contents depend on :math:`(A, \lda)`, according to +/// Sketch from the right in a SYMM-like operation /// -/// .. math:: -/// \mat(A)_{ij} = A[i + j \cdot \lda] = A[i \cdot \lda + j]. -/// -/// Note that the the "layout" parameter passed to this function is not used here. -/// That's because this function requires :math:`\mat(A)` to be stored in the format -/// of a general matrix (with both upper and lower triangles). -/// -/// This function's default behavior is to check that :math:`\mat(A)` is symmetric before -/// attempting sketching. That check can be skipped (at your own peril!) by calling this -/// function with sym_check_tol < 0. -/// -/// **What's** :math:`\mat(B)?` -/// -/// It's an :math:`n \times d` matrix, where :math:`n = \texttt{S.dist.n_cols}` -/// and :math:`d = \texttt{S.dist.n_rows}`. -/// Its precise contents depend on :math:`(B,\ldb)` and "layout." -/// -/// If layout == ColMajor, then -/// -/// .. math:: -/// \mat(B)_{ij} = B[i + j \cdot \ldb]. -/// -/// In this case, :math:`\ldb` must be :math:`\geq n.` -/// -/// If layout == RowMajor, then -/// -/// .. math:: -/// \mat(B)_{ij} = B[i \cdot \ldb + j]. -/// -/// In this case, :math:`\ldb` must be :math:`\geq d.` -/// -/// .. dropdown:: Full parameter descriptions -/// :animate: fade-in-slide-down -/// -/// layout - [in] -/// * Either Layout::ColMajor or Layout::RowMajor -/// * Matrix storage for :math:`\mat(B).` -/// -/// alpha - [in] -/// * A real scalar. -/// * If zero, then :math:`A` is not accessed. +/// .. math:: +/// \mat(B) = \alpha \cdot \underbrace{\mat(A)}_{n \times n} \cdot \underbrace{\submat(\mtxS)}_{n \times d} + \beta \cdot \underbrace{\mat(B)}_{n \times d}, \tag{$\star$} /// -/// A - [in] -/// * Pointer to a 1D array of real scalars. -/// * Defines :math:`\mat(A).` +/// where :math:`\alpha` and :math:`\beta` are real scalars and :math:`\mtxS` is a sketching operator. /// -/// lda - [in] -/// * A nonnegative integer. -/// * Leading dimension of :math:`\mat(A)` when reading from :math:`A.` +/// See the left-side overload above for the meaning of :math:`\mat(A)` and its +/// :math:`\texttt{uplo}` storage convention. The roles of :math:`d` and :math:`n` are mirrored: +/// :math:`d` is the embedding dimension (cols of :math:`\submat(\mtxS)`) and :math:`n` is the +/// order of :math:`\mat(A).` /// -/// S - [in] -/// * A DenseSkOp or SparseSkOp object. +/// @endverbatim +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + int64_t n, int64_t d, + T alpha, + const T* A, int64_t lda, + const SKOP &S, int64_t ro_s, int64_t co_s, + T beta, + T* B, int64_t ldb +); + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + int64_t n, int64_t d, + T alpha, + const T* A, int64_t lda, + const DenseSkOp &S, int64_t ro_s, int64_t co_s, + T beta, + T* B, int64_t ldb +) { + RandBLAS::dense::rsksy3(layout, uplo, n, d, alpha, A, lda, S, ro_s, co_s, beta, B, ldb); +} + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + int64_t n, int64_t d, + T alpha, + const T* A, int64_t lda, + const SparseSkOp &S, int64_t ro_s, int64_t co_s, + T beta, + T* B, int64_t ldb +) { + (void) layout; (void) uplo; (void) n; (void) d; (void) alpha; + (void) S; (void) ro_s; (void) co_s; (void) A; (void) lda; + (void) beta; (void) B; (void) ldb; + randblas_require( + false && + "RandBLAS::sketch_symmetric with a sparse sketching operator is the " + "Case-B kernel from project-plans/randblas-symm-plan.md. Not implemented." + ); +} + + +// MARK: FULL(S) + +// ============================================================================= +/// \fn sketch_symmetric(blas::Layout layout, blas::Uplo uplo, T alpha, +/// const SKOP &S, const T *A, int64_t lda, +/// T beta, T *B, int64_t ldb +/// ) +/// @verbatim embed:rst:leading-slashes +/// Sketch from the left in a SYMM-like operation, with :math:`\mtxS` used in full +/// (no submatrix offsets). /// -/// beta - [in] -/// * A real scalar. -/// * If zero, then :math:`B` need not be set on input. +/// .. math:: +/// \mat(B) = \alpha \cdot \underbrace{\mtxS}_{d \times n} \cdot \underbrace{\mat(A)}_{n \times n} + \beta \cdot \underbrace{\mat(B)}_{d \times n}, \tag{$\star$} /// -/// B - [in,out] -/// * Pointer to 1D array of real scalars. -/// * On entry, defines :math:`\mat(B)` -/// on the RIGHT-hand side of :math:`(\star).` -/// * On exit, defines :math:`\mat(B)` -/// on the LEFT-hand side of :math:`(\star).` +/// The dimensions :math:`d` and :math:`n` are taken from :math:`\mtxS` directly +/// (:math:`d = \mtxS.\ttt{dist.n\_rows}`, :math:`n = \mtxS.\ttt{dist.n\_cols}`). /// -/// ldb - [in] -/// * A nonnegative integer. -/// * Leading dimension of :math:`\mat(B)` when reading from :math:`B.` +/// See the submatrix overload above for the meaning of :math:`\mat(A)` and +/// :math:`\texttt{uplo}`. /// /// @endverbatim -template +template inline void sketch_symmetric( - // B = alpha*A*S + beta*B, where A is a symmetric matrix stored in the format of a general matrix. - blas::Layout layout, + blas::Layout layout, blas::Uplo uplo, T alpha, - const T* A, - int64_t lda, const SKOP &S, + const T* A, int64_t lda, T beta, - T* B, - int64_t ldb, - T sym_check_tol = 0 + T* B, int64_t ldb +); + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + T alpha, + const DenseSkOp &S, + const T* A, int64_t lda, + T beta, + T* B, int64_t ldb ) { - int64_t n = S.dist.n_rows; - int64_t d = S.dist.n_cols; - RandBLAS::util::require_symmetric(layout, A, n, lda, sym_check_tol); - sketch_general(layout, blas::Op::NoTrans, blas::Op::NoTrans, n, d, n, alpha, A, lda, S, 0, 0, beta, B, ldb); + int64_t d = S.dist.n_rows; + int64_t n = S.dist.n_cols; + RandBLAS::dense::lsksy3(layout, uplo, d, n, alpha, S, 0, 0, A, lda, beta, B, ldb); +} + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + T alpha, + const SparseSkOp &S, + const T* A, int64_t lda, + T beta, + T* B, int64_t ldb +) { + (void) layout; (void) uplo; (void) alpha; + (void) S; (void) A; (void) lda; + (void) beta; (void) B; (void) ldb; + randblas_require( + false && + "RandBLAS::sketch_symmetric with a sparse sketching operator is the " + "Case-B kernel from project-plans/randblas-symm-plan.md. Not implemented." + ); } // ============================================================================= -/// \fn sketch_symmetric(blas::Layout layout, T alpha, const SKOP &S, -/// const T *A, int64_t lda, T beta, T *B, int64_t ldb, T sym_check_tol = 0 -/// ) +/// \fn sketch_symmetric(blas::Layout layout, blas::Uplo uplo, T alpha, +/// const T *A, int64_t lda, const SKOP &S, +/// T beta, T *B, int64_t ldb +/// ) /// @verbatim embed:rst:leading-slashes -/// Check that :math:`\mat(A)` is symmetric up to tolerance :math:`\texttt{sym_check_tol}`, then sketch from the left in a SYMM-like operation -/// -/// .. math:: -/// \mat(B) = \alpha \cdot \mtxS \cdot \underbrace{\mat(A)}_{n \times n} + \beta \cdot \underbrace{\mat(B)}_{d \times n}, \tag{$\star$} -/// -/// where :math:`\alpha` and :math:`\beta` are real scalars and :math:`\mtxS` is a :math:`d \times n` sketching operator. -/// -/// .. dropdown:: FAQ -/// :animate: fade-in-slide-down -/// -/// **What's** :math:`\mat(A)?` -/// -/// It's a symmetric matrix of order :math:`n`. Its precise contents depend on :math:`(A, \lda)`, -/// according to -/// -/// .. math:: -/// \mat(A)_{ij} = A[i + j \cdot \lda] = A[i \cdot \lda + j]. -/// -/// Note that the the "layout" parameter passed to this function is not used here. -/// That's because this function requires :math:`\mat(A)` to be stored in the format -/// of a general matrix (with both upper and lower triangles). -/// -/// This function's default behavior is to check that :math:`\mat(A)` is symmetric before -/// attempting sketching. That check can be skipped (at your own peril!) by calling this -/// function with sym_check_tol < 0. +/// Sketch from the right in a SYMM-like operation, with :math:`\mtxS` used in full. /// -/// **What's** :math:`\mat(B)?` -/// -/// It's a :math:`d \times n` matrix. Its precise contents depend on :math:`(B,\ldb)` and "layout." -/// -/// If layout == ColMajor, then -/// -/// .. math:: -/// \mat(B)_{ij} = B[i + j \cdot \ldb]. -/// -/// In this case, :math:`\ldb` must be :math:`\geq d.` -/// -/// If layout == RowMajor, then -/// -/// .. math:: -/// \mat(B)_{ij} = B[i \cdot \ldb + j]. -/// -/// In this case, :math:`\ldb` must be :math:`\geq n.` -/// -/// .. dropdown:: Full parameter descriptions -/// :animate: fade-in-slide-down -/// -/// layout - [in] -/// * Either Layout::ColMajor or Layout::RowMajor -/// * Matrix storage for :math:`\mat(B).` -/// -/// alpha - [in] -/// * A real scalar. -/// * If zero, then :math:`A` is not accessed. -/// -/// S - [in] -/// * A DenseSkOp or SparseSkOp object. -/// -/// A - [in] -/// * Pointer to a 1D array of real scalars. -/// * Defines :math:`\mat(A).` -/// -/// lda - [in] -/// * A nonnegative integer. -/// * Leading dimension of :math:`\mat(A)` when reading from :math:`A.` +/// .. math:: +/// \mat(B) = \alpha \cdot \underbrace{\mat(A)}_{n \times n} \cdot \underbrace{\mtxS}_{n \times d} + \beta \cdot \underbrace{\mat(B)}_{n \times d}, \tag{$\star$} /// -/// beta - [in] -/// * A real scalar. -/// * If zero, then :math:`B` need not be set on input. +/// The dimensions :math:`n` and :math:`d` are taken from :math:`\mtxS` directly +/// (:math:`n = \mtxS.\ttt{dist.n\_rows}`, :math:`d = \mtxS.\ttt{dist.n\_cols}`). /// -/// B - [in,out] -/// * Pointer to 1D array of real scalars. -/// * On entry, defines :math:`\mat(B)` -/// on the RIGHT-hand side of :math:`(\star).` -/// * On exit, defines :math:`\mat(B)` -/// on the LEFT-hand side of :math:`(\star).` -/// -/// ldb - [in] -/// * A nonnegative integer. -/// * Leading dimension of :math:`\mat(B)` when reading from :math:`B.` +/// See the submatrix overload above for the meaning of :math:`\mat(A)` and +/// :math:`\texttt{uplo}`. /// /// @endverbatim -template +template inline void sketch_symmetric( - // B = alpha*S*A + beta*B - blas::Layout layout, + blas::Layout layout, blas::Uplo uplo, T alpha, + const T* A, int64_t lda, const SKOP &S, - const T* A, - int64_t lda, T beta, - T* B, - int64_t ldb, - T sym_check_tol = 0 + T* B, int64_t ldb +); + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + T alpha, + const T* A, int64_t lda, + const DenseSkOp &S, + T beta, + T* B, int64_t ldb ) { - int64_t d = S.dist.n_rows; - int64_t n = S.dist.n_cols; - RandBLAS::util::require_symmetric(layout, A, n, lda, sym_check_tol); - sketch_general(layout, blas::Op::NoTrans, blas::Op::NoTrans, d, n, n, alpha, S, 0, 0, A, lda, beta, B, ldb); + int64_t n = S.dist.n_rows; + int64_t d = S.dist.n_cols; + RandBLAS::dense::rsksy3(layout, uplo, n, d, alpha, A, lda, S, 0, 0, beta, B, ldb); +} + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + T alpha, + const T* A, int64_t lda, + const SparseSkOp &S, + T beta, + T* B, int64_t ldb +) { + (void) layout; (void) uplo; (void) alpha; + (void) S; (void) A; (void) lda; + (void) beta; (void) B; (void) ldb; + randblas_require( + false && + "RandBLAS::sketch_symmetric with a sparse sketching operator is the " + "Case-B kernel from project-plans/randblas-symm-plan.md. Not implemented." + ); } } // end namespace RandBLAS diff --git a/RandBLAS/util.hh b/RandBLAS/util.hh index 597b1e55..5ee923b2 100644 --- a/RandBLAS/util.hh +++ b/RandBLAS/util.hh @@ -122,7 +122,11 @@ void flip_layout(blas::Layout layout_in, int64_t m, int64_t n, std::vector &A /// for all :math:`i,j \in \\{0,\ldots,n-1\\}.` An error is raised if any such check fails. /// This function returns immediately without performing any checks if :math:`\ttt{tol} < 0.` /// @endverbatim -/// sketch_symmetric calls this function with \math{\ttt{tol} = 0} by default. +/// (Historical note: pre-Phase-1 of project-plans/randblas-symm-plan.md, +/// sketch_symmetric called this function with \math{\ttt{tol} = 0} by default. +/// With the SYMM-based dispatch, only the triangle named by \math{\ttt{uplo}} +/// is read, so the symmetry check is no longer invoked from sketch_symmetric. +/// require_symmetric remains available as a standalone validator.) /// template void require_symmetric(blas::Layout layout, const T* A, int64_t n, int64_t lda, T tol) { diff --git a/test/linops/test_sketch_symmetric.cc b/test/linops/test_sketch_symmetric.cc index 5ee3269f..94383cb0 100644 --- a/test/linops/test_sketch_symmetric.cc +++ b/test/linops/test_sketch_symmetric.cc @@ -60,14 +60,15 @@ void random_symmetric_mat(int64_t n, T* A, int64_t lda, STATE s) { template blas::Side sketch_symmetric_side( - blas::Side side_skop, blas::Layout layout, int64_t rows_out, int64_t cols_out, + blas::Side side_skop, blas::Layout layout, blas::Uplo uplo, + int64_t rows_out, int64_t cols_out, T alpha, const T* A, int64_t lda, SKOP &S, int64_t ro_s, int64_t co_s, T beta, T* B, int64_t ldb ) { if (side_skop == blas::Side::Left) { - RandBLAS::sketch_symmetric(layout, rows_out, cols_out, alpha, S, ro_s, co_s, A, lda, beta, B, ldb); + RandBLAS::sketch_symmetric(layout, uplo, rows_out, cols_out, alpha, S, ro_s, co_s, A, lda, beta, B, ldb); return blas::Side::Right; } else { - RandBLAS::sketch_symmetric(layout, rows_out, cols_out, alpha, A, lda, S, ro_s, co_s, beta, B, ldb); + RandBLAS::sketch_symmetric(layout, uplo, rows_out, cols_out, alpha, A, lda, S, ro_s, co_s, beta, B, ldb); return blas::Side::Left; } } @@ -88,7 +89,8 @@ class TestSketchSymmetric : public ::testing::Test { template static void test_same_layouts( - uint32_t seed_a, uint32_t seed_skop, Axis major_axis, T alpha, int64_t d, int64_t n, int64_t lda, T beta, blas::Side side_skop + uint32_t seed_a, uint32_t seed_skop, Axis major_axis, T alpha, int64_t d, int64_t n, int64_t lda, T beta, blas::Side side_skop, + blas::Uplo uplo = blas::Uplo::Upper ) { auto [rows_out, cols_out] = dims_of_sketch_symmetric_output(d, n, side_skop); std::vector A(lda*lda, 0.0); @@ -104,9 +106,9 @@ class TestSketchSymmetric : public ::testing::Test { std::vector B_expect(B_actual); // Compute the actual output - auto side_a = sketch_symmetric_side(side_skop, S.layout, rows_out, cols_out, alpha, A.data(), lda, S, 0, 0, beta, B_actual.data(), ldb); + auto side_a = sketch_symmetric_side(side_skop, S.layout, uplo, rows_out, cols_out, alpha, A.data(), lda, S, 0, 0, beta, B_actual.data(), ldb); // Compute the expected output - blas::symm(S.layout, side_a, Uplo::Upper, rows_out, cols_out, alpha, A.data(), lda, S.buff, lds, beta, B_expect.data(), ldb); + blas::symm(S.layout, side_a, uplo, rows_out, cols_out, alpha, A.data(), lda, S.buff, lds, beta, B_expect.data(), ldb); auto msg = RandBLAS::testing::matrices_approx_equal( S.layout, blas::Op::NoTrans, rows_out, cols_out, B_actual.data(), ldb, B_expect.data(), ldb, @@ -120,7 +122,8 @@ class TestSketchSymmetric : public ::testing::Test { template static void test_opposing_layouts( - uint32_t seed_a, uint32_t seed_skop, Axis major_axis, T alpha, int64_t d, int64_t n, int64_t lda, T beta, blas::Side side_skop + uint32_t seed_a, uint32_t seed_skop, Axis major_axis, T alpha, int64_t d, int64_t n, int64_t lda, T beta, blas::Side side_skop, + blas::Uplo uplo = blas::Uplo::Upper ) { auto [rows_out, cols_out] = dims_of_sketch_symmetric_output(d, n, side_skop); std::vector A(lda*lda, 0.0); @@ -144,11 +147,11 @@ class TestSketchSymmetric : public ::testing::Test { RandBLAS::fill_dense(D, B_actual.data(), RNGState(seed_b)); std::vector B_expect(B_actual); // Compute the actual output - auto side_a = sketch_symmetric_side(side_skop, layout_B, rows_out, cols_out, alpha, A.data(), lda, S, 0, 0, beta, B_actual.data(), ldb); + auto side_a = sketch_symmetric_side(side_skop, layout_B, uplo, rows_out, cols_out, alpha, A.data(), lda, S, 0, 0, beta, B_actual.data(), ldb); // Compute the expected output std::vector S_flipped(S.buff, S.buff + d*n); RandBLAS::util::flip_layout(S.layout, rows_out, cols_out, S_flipped, lds_init, ldb); - blas::symm(layout_B, side_a, Uplo::Upper, rows_out, cols_out, alpha, A.data(), lda, S_flipped.data(), ldb, beta, B_expect.data(), ldb); + blas::symm(layout_B, side_a, uplo, rows_out, cols_out, alpha, A.data(), lda, S_flipped.data(), ldb, beta, B_expect.data(), ldb); auto msg = RandBLAS::testing::matrices_approx_equal( layout_B, blas::Op::NoTrans, rows_out, cols_out, B_actual.data(), ldb, B_expect.data(), ldb, @@ -160,36 +163,12 @@ class TestSketchSymmetric : public ::testing::Test { return; } - template - static void test_error_on_asymmetric() { - // Build a 3x3 non-symmetric matrix (A[0,1] != A[1,0]) and verify - // that sketch_symmetric raises RandBLAS::Error due to the symmetry check. - int64_t n = 3; - int64_t d = 2; - // Column-major 3x3 matrix: symmetric except A(0,1)=5 vs A(1,0)=0. - // col0 col1 col2 - // 1.0 5.0 0.0 - // 0.0 2.0 0.0 - // 0.0 0.0 3.0 - std::vector A = { - 1.0, 0.0, 0.0, - 5.0, 2.0, 0.0, - 0.0, 0.0, 3.0 - }; - DenseDist D(d, n, ScalarDist::Uniform, Axis::Short); - DenseSkOp S(D, 42); - RandBLAS::fill_dense(S); - std::vector B(d * n, 0.0); - try { - RandBLAS::sketch_symmetric(Layout::ColMajor, d, n, (T)1.0, S, 0, 0, A.data(), n, (T)0.0, B.data(), d); - FAIL() << "Expected RandBLAS::Error for asymmetric matrix"; - } catch (const RandBLAS::Error& e) { - std::string msg = e.what(); - EXPECT_NE(msg.find("Symmetry check failed"), std::string::npos) - << "Error message did not mention symmetry check: " << msg; - } - return; - } + // Note: an earlier `test_error_on_asymmetric` test verified that + // sketch_symmetric threw on asymmetric input via require_symmetric. With + // the SYMM-based dispatch (project-plans/randblas-symm-plan.md, Phase 1), + // only the triangle named by `uplo` is read --- the opposite triangle + // contributing wrong values is undefined behavior at the caller, not a + // RandBLAS check. The test was removed accordingly. }; @@ -391,9 +370,3 @@ TEST_F(TestSketchSymmetric, right_lift_opposing_layouts) { test_opposing_layouts(31, 33, Axis::Long, 0.5, 50, 10, 19, -1.0, blas::Side::Right); } -// MARK: validation errors - -TEST_F(TestSketchSymmetric, symmetry_check_fails_for_asymmetric_matrix) { - test_error_on_asymmetric(); - test_error_on_asymmetric(); -} From c73bd19cc0504d1f9e4be95d91537605bfe922c3 Mon Sep 17 00:00:00 2001 From: mmelnich Date: Wed, 13 May 2026 18:33:53 -0400 Subject: [PATCH 02/19] Add Symmetric wrapper for sparse symmetric matrices 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 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 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 . 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.) --- RandBLAS/sparse_data/base.hh | 1 + RandBLAS/sparse_data/symmetric.hh | 97 +++++++++++++++++++ test/CMakeLists.txt | 1 + test/datastructures/test_symmetric_wrapper.cc | 94 ++++++++++++++++++ 4 files changed, 193 insertions(+) create mode 100644 RandBLAS/sparse_data/symmetric.hh create mode 100644 test/datastructures/test_symmetric_wrapper.cc diff --git a/RandBLAS/sparse_data/base.hh b/RandBLAS/sparse_data/base.hh index 7eb4df2a..8926b81d 100644 --- a/RandBLAS/sparse_data/base.hh +++ b/RandBLAS/sparse_data/base.hh @@ -31,6 +31,7 @@ #include "RandBLAS/config.h" #include "RandBLAS/base.hh" #include +#include // std::iota, used below #ifdef __cpp_concepts #include diff --git a/RandBLAS/sparse_data/symmetric.hh b/RandBLAS/sparse_data/symmetric.hh new file mode 100644 index 00000000..dee1a52a --- /dev/null +++ b/RandBLAS/sparse_data/symmetric.hh @@ -0,0 +1,97 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include "RandBLAS/sparse_data/base.hh" +#include "RandBLAS/exceptions.hh" + + +namespace RandBLAS::sparse_data { + +// ============================================================================= +/// Lightweight non-owning wrapper marking a SparseMatrix as symmetric with a +/// stored triangle. +/// +/// @verbatim embed:rst:leading-slashes +/// Holds: +/// - A const reference to the underlying SparseMatrix :math:`A.` +/// - :math:`\ttt{blas::Uplo uplo}`: names the triangle of :math:`A` that is +/// structurally populated. The opposite triangle is implied by symmetry. +/// +/// The wrapper performs **no validation** of the matrix contents --- it is the +/// caller's responsibility to guarantee that the named triangle is correctly +/// populated and that the opposite triangle is either structurally absent or +/// will be ignored by the SYMM-aware consumer (kernels in this library that +/// accept a ``Symmetric`` agree to consult ``uplo`` and respect this +/// contract). Construction does enforce that :math:`A` is square +/// (``A.n_rows == A.n_cols``) via :math:`\ttt{randblas\_require}`. +/// +/// This wrapper is the carrier for symmetric sparse matrices into the +/// :math:`\ttt{spsymm}`-family kernels (project-plans/randblas-symm-plan.md +/// Case C and beyond). It is intentionally separate from the underlying +/// :math:`\ttt{SparseMatrix}` concept so that calling +/// :math:`\ttt{spmm}` / :math:`\ttt{spgemm}` with a ``Symmetric`` +/// argument fails to compile rather than silently treating the matrix as +/// general --- a type-system guard against accidental loss of symmetry. +/// +/// The wrapper is non-owning: it holds a reference to :math:`A`, not a copy. +/// The caller must keep :math:`A` alive for the lifetime of the wrapper. +/// @endverbatim +template +struct Symmetric { + const SpMat& A; + const blas::Uplo uplo; + + using scalar_t = typename SpMat::scalar_t; + using index_t = typename SpMat::index_t; + + Symmetric(const SpMat& A_in, blas::Uplo uplo_in) : A(A_in), uplo(uplo_in) { + randblas_require(A_in.n_rows == A_in.n_cols); + } +}; + + +// ============================================================================= +/// Construct a ``Symmetric`` wrapper around ``A`` with the named triangle. +/// Syntactic sugar for the constructor; lets callers write +/// ``as_symmetric(A, blas::Uplo::Upper)`` and rely on template argument +/// deduction. +template +inline Symmetric as_symmetric(const SpMat& A, blas::Uplo uplo) { + return Symmetric(A, uplo); +} + +} // end namespace RandBLAS::sparse_data + + +namespace RandBLAS { + using RandBLAS::sparse_data::Symmetric; + using RandBLAS::sparse_data::as_symmetric; +} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index cdf202a7..60de660c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -38,6 +38,7 @@ if (GTest_FOUND) datastructures/test_csr_matrix.cc datastructures/test_coo_matrix.cc datastructures/test_spmat_conversions.cc + datastructures/test_symmetric_wrapper.cc linops/test_spgemm.cc linops/test_sparse_trsm.cc diff --git a/test/datastructures/test_symmetric_wrapper.cc b/test/datastructures/test_symmetric_wrapper.cc new file mode 100644 index 00000000..bce1d076 --- /dev/null +++ b/test/datastructures/test_symmetric_wrapper.cc @@ -0,0 +1,94 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#include "RandBLAS/sparse_data/coo_matrix.hh" +#include "RandBLAS/sparse_data/csr_matrix.hh" +#include "RandBLAS/sparse_data/csc_matrix.hh" +#include "RandBLAS/sparse_data/symmetric.hh" +#include "RandBLAS/exceptions.hh" + +#include +#include + +using namespace RandBLAS::sparse_data; + + +class TestSymmetricWrapper : public ::testing::Test {}; + + +// Construct from each of the three sparse formats and confirm field access + +// trait aliases. The wrapper does not require any actual data, since it only +// stores a reference + uplo. + +TEST_F(TestSymmetricWrapper, ConstructFromCOO) { + COOMatrix A(4, 4); + auto wrapped = as_symmetric(A, blas::Uplo::Upper); + EXPECT_EQ(&wrapped.A, &A); + EXPECT_EQ(wrapped.uplo, blas::Uplo::Upper); + EXPECT_EQ(wrapped.A.n_rows, 4); + EXPECT_EQ(wrapped.A.n_cols, 4); + static_assert(std::is_same_v, + "Symmetric>::scalar_t must be double"); +} + +TEST_F(TestSymmetricWrapper, ConstructFromCSR) { + CSRMatrix A(5, 5); + auto wrapped = as_symmetric(A, blas::Uplo::Lower); + EXPECT_EQ(&wrapped.A, &A); + EXPECT_EQ(wrapped.uplo, blas::Uplo::Lower); + EXPECT_EQ(wrapped.A.n_rows, 5); + static_assert(std::is_same_v, + "Symmetric>::scalar_t must be float"); +} + +TEST_F(TestSymmetricWrapper, ConstructFromCSC) { + CSCMatrix A(3, 3); + auto wrapped = as_symmetric(A, blas::Uplo::Upper); + EXPECT_EQ(&wrapped.A, &A); + EXPECT_EQ(wrapped.uplo, blas::Uplo::Upper); + EXPECT_EQ(wrapped.A.n_rows, 3); +} + +TEST_F(TestSymmetricWrapper, NamespacingExportsAtRandBLASScope) { + // The wrapper and the helper must be accessible via the top-level + // RandBLAS:: namespace (re-exported from sparse_data::). + COOMatrix A(2, 2); + RandBLAS::Symmetric> wrapped_typed(A, blas::Uplo::Upper); + EXPECT_EQ(&wrapped_typed.A, &A); + auto wrapped_sugar = RandBLAS::as_symmetric(A, blas::Uplo::Lower); + EXPECT_EQ(wrapped_sugar.uplo, blas::Uplo::Lower); +} + +TEST_F(TestSymmetricWrapper, NonSquareRejectedAtConstruction) { + COOMatrix A(3, 4); // not square + EXPECT_THROW({ + auto wrapped = as_symmetric(A, blas::Uplo::Upper); + (void) wrapped; + }, RandBLAS::Error); +} From a26f1db855f2c1556a913ceef8d28989bed9c166 Mon Sep 17 00:00:00 2001 From: mmelnich Date: Wed, 13 May 2026 19:28:51 -0400 Subject: [PATCH 03/19] Add spsymm: symmetric sparse x dense matmul with MKL fast path and per-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 A_sym, B, ldb, beta, Y, ldy) Wrapper-routing overload. Extracts uplo from the Symmetric 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 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. --- RandBLAS/sparse_data/coo_spsymm_impl.hh | 99 ++++++++++ RandBLAS/sparse_data/csc_spsymm_impl.hh | 109 +++++++++++ RandBLAS/sparse_data/csr_spsymm_impl.hh | 151 +++++++++++++++ RandBLAS/sparse_data/mkl_spsymm_impl.hh | 131 +++++++++++++ RandBLAS/sparse_data/spsymm_dispatch.hh | 174 ++++++++++++++++++ test/CMakeLists.txt | 1 + test/linops/test_spsymm.cc | 232 ++++++++++++++++++++++++ 7 files changed, 897 insertions(+) create mode 100644 RandBLAS/sparse_data/coo_spsymm_impl.hh create mode 100644 RandBLAS/sparse_data/csc_spsymm_impl.hh create mode 100644 RandBLAS/sparse_data/csr_spsymm_impl.hh create mode 100644 RandBLAS/sparse_data/mkl_spsymm_impl.hh create mode 100644 RandBLAS/sparse_data/spsymm_dispatch.hh create mode 100644 test/linops/test_spsymm.cc diff --git a/RandBLAS/sparse_data/coo_spsymm_impl.hh b/RandBLAS/sparse_data/coo_spsymm_impl.hh new file mode 100644 index 00000000..981262a8 --- /dev/null +++ b/RandBLAS/sparse_data/coo_spsymm_impl.hh @@ -0,0 +1,99 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include "RandBLAS/exceptions.hh" +#include "RandBLAS/sparse_data/base.hh" +#include "RandBLAS/sparse_data/coo_matrix.hh" +#include "RandBLAS/sparse_data/csr_spsymm_impl.hh" // for internal::apply_beta_scale +#include + +namespace RandBLAS::sparse_data { + +// ============================================================================= +/// COO fallback for symmetric sparse-times-dense. +/// +/// Iterates over the (row, col, val) triples directly. For uplo=Upper, +/// structurally populated entries satisfy row <= col; for uplo=Lower, +/// row >= col. Entries outside the named triangle are silently skipped. +/// No assumption on the order of the COO entries. +template +void coo_spsymm( + blas::Layout layout, + blas::Side side, + blas::Uplo uplo, + int64_t m, int64_t n, + T alpha, + const COOMatrix& A, + const T* B, int64_t ldb, + T beta, + T* Y, int64_t ldy +) { + randblas_require(A.n_rows == A.n_cols); + int64_t k = (side == blas::Side::Left) ? m : n; + randblas_require(A.n_rows == k); + + internal::apply_beta_scale(layout, m, n, beta, Y, ldy); + if (alpha == T(0)) return; + + const bool col_major = (layout == blas::Layout::ColMajor); + + for (int64_t p = 0; p < A.nnz; ++p) { + sint_t i = A.rows[p]; + sint_t j = A.cols[p]; + if (uplo == blas::Uplo::Upper && j < i) continue; + if (uplo == blas::Uplo::Lower && j > i) continue; + T av = alpha * A.vals[p]; + + if (side == blas::Side::Left) { + if (col_major) { + blas::axpy(n, av, &B[j], ldb, &Y[i], ldy); + if (j != i) + blas::axpy(n, av, &B[i], ldb, &Y[j], ldy); + } else { + blas::axpy(n, av, &B[j*ldb], 1, &Y[i*ldy], 1); + if (j != i) + blas::axpy(n, av, &B[i*ldb], 1, &Y[j*ldy], 1); + } + } else { + if (col_major) { + blas::axpy(m, av, &B[i*ldb], 1, &Y[j*ldy], 1); + if (j != i) + blas::axpy(m, av, &B[j*ldb], 1, &Y[i*ldy], 1); + } else { + blas::axpy(m, av, &B[i], ldb, &Y[j], ldy); + if (j != i) + blas::axpy(m, av, &B[j], ldb, &Y[i], ldy); + } + } + } +} + +} // namespace RandBLAS::sparse_data diff --git a/RandBLAS/sparse_data/csc_spsymm_impl.hh b/RandBLAS/sparse_data/csc_spsymm_impl.hh new file mode 100644 index 00000000..e860df6a --- /dev/null +++ b/RandBLAS/sparse_data/csc_spsymm_impl.hh @@ -0,0 +1,109 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include "RandBLAS/exceptions.hh" +#include "RandBLAS/sparse_data/base.hh" +#include "RandBLAS/sparse_data/csc_matrix.hh" +#include "RandBLAS/sparse_data/csr_spsymm_impl.hh" // for internal::apply_beta_scale +#include + +namespace RandBLAS::sparse_data { + +// ============================================================================= +/// CSC fallback for symmetric sparse-times-dense. +/// +/// Same semantics as the CSR variant; iteration is column-driven via +/// colptr/rowidxs/vals. For CSC stored with uplo=Upper, structurally +/// populated entries satisfy row <= col; for uplo=Lower, row >= col. +/// Entries outside the named triangle are silently skipped. +template +void csc_spsymm( + blas::Layout layout, + blas::Side side, + blas::Uplo uplo, + int64_t m, int64_t n, + T alpha, + const CSCMatrix& A, + const T* B, int64_t ldb, + T beta, + T* Y, int64_t ldy +) { + randblas_require(A.n_rows == A.n_cols); + int64_t k = (side == blas::Side::Left) ? m : n; + randblas_require(A.n_rows == k); + + internal::apply_beta_scale(layout, m, n, beta, Y, ldy); + if (alpha == T(0)) return; + + const bool col_major = (layout == blas::Layout::ColMajor); + + if (side == blas::Side::Left) { + // Y = alpha * A * B + ... (A is m-by-m) + for (int64_t j = 0; j < m; ++j) { + for (int64_t p = A.colptr[j]; p < A.colptr[j+1]; ++p) { + sint_t i = A.rowidxs[p]; + if (uplo == blas::Uplo::Upper && i > j) continue; + if (uplo == blas::Uplo::Lower && i < j) continue; + T av = alpha * A.vals[p]; + if (col_major) { + blas::axpy(n, av, &B[j], ldb, &Y[i], ldy); + if (i != j) + blas::axpy(n, av, &B[i], ldb, &Y[j], ldy); + } else { + blas::axpy(n, av, &B[j*ldb], 1, &Y[i*ldy], 1); + if (i != j) + blas::axpy(n, av, &B[i*ldb], 1, &Y[j*ldy], 1); + } + } + } + } else { + // Y = alpha * B * A + ... (A is n-by-n) + for (int64_t j = 0; j < n; ++j) { + for (int64_t p = A.colptr[j]; p < A.colptr[j+1]; ++p) { + sint_t i = A.rowidxs[p]; + if (uplo == blas::Uplo::Upper && i > j) continue; + if (uplo == blas::Uplo::Lower && i < j) continue; + T av = alpha * A.vals[p]; + if (col_major) { + blas::axpy(m, av, &B[i*ldb], 1, &Y[j*ldy], 1); + if (i != j) + blas::axpy(m, av, &B[j*ldb], 1, &Y[i*ldy], 1); + } else { + blas::axpy(m, av, &B[i], ldb, &Y[j], ldy); + if (i != j) + blas::axpy(m, av, &B[j], ldb, &Y[i], ldy); + } + } + } + } +} + +} // namespace RandBLAS::sparse_data diff --git a/RandBLAS/sparse_data/csr_spsymm_impl.hh b/RandBLAS/sparse_data/csr_spsymm_impl.hh new file mode 100644 index 00000000..67fb2562 --- /dev/null +++ b/RandBLAS/sparse_data/csr_spsymm_impl.hh @@ -0,0 +1,151 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include "RandBLAS/exceptions.hh" +#include "RandBLAS/sparse_data/base.hh" +#include "RandBLAS/sparse_data/csr_matrix.hh" +#include +#include + +namespace RandBLAS::sparse_data { + +namespace internal { + +// Apply Y <- beta * Y for an m-by-n dense matrix with leading dimension ldy. +// Used by the spsymm fallback kernels to fold beta-scaling out of the +// nonzero-traversal loop. +template +inline void apply_beta_scale(blas::Layout layout, int64_t m, int64_t n, T beta, T* Y, int64_t ldy) { + if (beta == T(1)) return; + if (beta == T(0)) { + if (layout == blas::Layout::ColMajor) { + for (int64_t j = 0; j < n; ++j) + std::fill(Y + j*ldy, Y + j*ldy + m, T(0)); + } else { + for (int64_t i = 0; i < m; ++i) + std::fill(Y + i*ldy, Y + i*ldy + n, T(0)); + } + return; + } + if (layout == blas::Layout::ColMajor) { + for (int64_t j = 0; j < n; ++j) + blas::scal(m, beta, Y + j*ldy, 1); + } else { + for (int64_t i = 0; i < m; ++i) + blas::scal(n, beta, Y + i*ldy, 1); + } +} + +} // namespace internal + + +// ============================================================================= +/// CSR fallback for the symmetric sparse-times-dense kernel. +/// +/// Computes Y = alpha * op(A) * B + beta * Y (side=Left) or +/// Y = alpha * B * op(A) + beta * Y (side=Right), +/// where A is symmetric and only the triangle named by uplo is structurally +/// stored. The off-diagonal entries contribute twice (the (i,j) and the +/// implied (j,i)); diagonal entries contribute once. Entries that fall +/// outside the named triangle are silently skipped, so the kernel is robust +/// against callers who store both triangles by mistake (it just behaves like +/// the "correctly stored" case). +/// +/// Inner-loop updates use blas::axpy on slices of B and Y. +template +void csr_spsymm( + blas::Layout layout, + blas::Side side, + blas::Uplo uplo, + int64_t m, int64_t n, + T alpha, + const CSRMatrix& A, + const T* B, int64_t ldb, + T beta, + T* Y, int64_t ldy +) { + randblas_require(A.n_rows == A.n_cols); + int64_t k = (side == blas::Side::Left) ? m : n; + randblas_require(A.n_rows == k); + + internal::apply_beta_scale(layout, m, n, beta, Y, ldy); + if (alpha == T(0)) return; + + const bool col_major = (layout == blas::Layout::ColMajor); + + if (side == blas::Side::Left) { + // Y = alpha * A * B + ... (A is m-by-m) + // For stored (i, j, v): + // Y[i, :] += alpha * v * B[j, :] + // if i != j: Y[j, :] += alpha * v * B[i, :] + for (int64_t i = 0; i < m; ++i) { + for (int64_t p = A.rowptr[i]; p < A.rowptr[i+1]; ++p) { + sint_t j = A.colidxs[p]; + if (uplo == blas::Uplo::Upper && j < i) continue; + if (uplo == blas::Uplo::Lower && j > i) continue; + T av = alpha * A.vals[p]; + if (col_major) { + blas::axpy(n, av, &B[j], ldb, &Y[i], ldy); + if (j != i) + blas::axpy(n, av, &B[i], ldb, &Y[j], ldy); + } else { + blas::axpy(n, av, &B[j*ldb], 1, &Y[i*ldy], 1); + if (j != i) + blas::axpy(n, av, &B[i*ldb], 1, &Y[j*ldy], 1); + } + } + } + } else { + // Y = alpha * B * A + ... (A is n-by-n) + // For stored (i, j, v): + // Y[:, j] += alpha * v * B[:, i] + // if i != j: Y[:, i] += alpha * v * B[:, j] + for (int64_t i = 0; i < n; ++i) { + for (int64_t p = A.rowptr[i]; p < A.rowptr[i+1]; ++p) { + sint_t j = A.colidxs[p]; + if (uplo == blas::Uplo::Upper && j < i) continue; + if (uplo == blas::Uplo::Lower && j > i) continue; + T av = alpha * A.vals[p]; + if (col_major) { + blas::axpy(m, av, &B[i*ldb], 1, &Y[j*ldy], 1); + if (j != i) + blas::axpy(m, av, &B[j*ldb], 1, &Y[i*ldy], 1); + } else { + blas::axpy(m, av, &B[i], ldb, &Y[j], ldy); + if (j != i) + blas::axpy(m, av, &B[j], ldb, &Y[i], ldy); + } + } + } + } +} + +} // namespace RandBLAS::sparse_data diff --git a/RandBLAS/sparse_data/mkl_spsymm_impl.hh b/RandBLAS/sparse_data/mkl_spsymm_impl.hh new file mode 100644 index 00000000..260514b6 --- /dev/null +++ b/RandBLAS/sparse_data/mkl_spsymm_impl.hh @@ -0,0 +1,131 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include "RandBLAS/config.h" + +#if defined(RandBLAS_HAS_MKL) + +#if defined(BLAS_ILP64) && !defined(MKL_ILP64) +#define MKL_ILP64 +#endif + +#include +#include + +#include "RandBLAS/exceptions.hh" +#include "RandBLAS/sparse_data/base.hh" +#include "RandBLAS/sparse_data/coo_matrix.hh" +#include "RandBLAS/sparse_data/csr_matrix.hh" +#include "RandBLAS/sparse_data/csc_matrix.hh" +#include "RandBLAS/sparse_data/mkl_spmm_impl.hh" // reuse make_mkl_handle, to_mkl_layout, check_mkl_status + +namespace RandBLAS::sparse_data::mkl { + +// ============================================================================ +// MKL-accelerated symmetric SpMM: Y = alpha * A * B + beta * Y (side=Left only). +// A is symmetric sparse (one triangle stored, named by uplo); B and Y dense. +// +// Returns true if MKL handled the call, false to signal fallback to the +// hand-rolled per-format kernel. MKL applies alpha and beta internally; the +// caller therefore passes them through unchanged. +// +// Known limitations that trigger a fallback: +// - side == Right: MKL's mkl_sparse_d_mm has no side parameter; the +// symmetric matrix is always on the left of the dense block. The +// transpose-trick to express side=Right via layout flips depends on +// leading-dim assumptions that don't always hold; safer to fall back. +// - CSC format: MKL's mkl_sparse_d_mm returns NOT_SUPPORTED for CSC even +// with a symmetric descriptor (mirrors the behavior in mkl_left_spmm). +// - Index type mismatched with MKL_INT. +// ============================================================================ +template +bool mkl_spsymm( + blas::Layout layout, + blas::Side side, + 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 +) { + using sint_t = typename SpMat::index_t; + constexpr bool is_csc = std::is_same_v>; + + if (side != blas::Side::Left) + return false; + if constexpr (is_csc) + return false; + + auto h = make_mkl_handle(A); + + struct matrix_descr descr; + descr.type = SPARSE_MATRIX_TYPE_SYMMETRIC; + descr.mode = (uplo == blas::Uplo::Upper) + ? SPARSE_FILL_MODE_UPPER + : SPARSE_FILL_MODE_LOWER; + descr.diag = SPARSE_DIAG_NON_UNIT; + + sparse_status_t status; + if constexpr (std::is_same_v) { + status = mkl_sparse_d_mm( + SPARSE_OPERATION_NON_TRANSPOSE, alpha, h.handle, descr, + to_mkl_layout(layout), + B, (MKL_INT)n, (MKL_INT)ldb, + beta, Y, (MKL_INT)ldy + ); + } else if constexpr (std::is_same_v) { + status = mkl_sparse_s_mm( + SPARSE_OPERATION_NON_TRANSPOSE, alpha, h.handle, descr, + to_mkl_layout(layout), + B, (MKL_INT)n, (MKL_INT)ldb, + beta, Y, (MKL_INT)ldy + ); + } else { + static_assert(sizeof(T) == 0, "MKL sparse BLAS only supports float and double."); + } + + // Some MKL versions return NOT_SUPPORTED for combinations we couldn't + // predict. Don't throw -- signal fallback. + if (status == SPARSE_STATUS_NOT_SUPPORTED) + return false; + + check_mkl_status(status, "mkl_sparse_mm (symmetric)"); + (void) m; // m is implied by A.n_rows; kept for signature symmetry + return true; +} + +} // namespace RandBLAS::sparse_data::mkl + +#endif // RandBLAS_HAS_MKL diff --git a/RandBLAS/sparse_data/spsymm_dispatch.hh b/RandBLAS/sparse_data/spsymm_dispatch.hh new file mode 100644 index 00000000..b12b3f65 --- /dev/null +++ b/RandBLAS/sparse_data/spsymm_dispatch.hh @@ -0,0 +1,174 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include "RandBLAS/base.hh" +#include "RandBLAS/exceptions.hh" +#include "RandBLAS/sparse_data/base.hh" +#include "RandBLAS/sparse_data/coo_matrix.hh" +#include "RandBLAS/sparse_data/csr_matrix.hh" +#include "RandBLAS/sparse_data/csc_matrix.hh" +#include "RandBLAS/sparse_data/coo_spsymm_impl.hh" +#include "RandBLAS/sparse_data/csr_spsymm_impl.hh" +#include "RandBLAS/sparse_data/csc_spsymm_impl.hh" +#include "RandBLAS/sparse_data/symmetric.hh" +#include "RandBLAS/config.h" + +#if defined(RandBLAS_HAS_MKL) +#include "RandBLAS/sparse_data/mkl_spsymm_impl.hh" +#endif + +#include + + +namespace RandBLAS::sparse_data { + +// ============================================================================= +/// Dispatched symmetric sparse-times-dense kernel (Case C in +/// project-plans/randblas-symm-plan.md). +/// +/// @verbatim embed:rst:leading-slashes +/// Computes +/// +/// .. math:: +/// \mat(Y) = \alpha \cdot \op_{\ttt{side}}(\mat(A), \mat(B)) + \beta \cdot \mat(Y), +/// +/// where: +/// +/// - :math:`\mat(A)` is a sparse symmetric matrix stored in COO, CSR, or +/// CSC format. Only the triangle named by :math:`\ttt{uplo}` is read. +/// The opposite triangle is implied by symmetry. The matrix is required +/// to be square. +/// - :math:`\mat(B)` and :math:`\mat(Y)` are dense, both m-by-n. +/// - For :math:`\ttt{side} = \ttt{Left}`, :math:`\mat(A)` is m-by-m and the +/// operation is :math:`\mat(Y) = \alpha A B + \beta Y`. +/// - For :math:`\ttt{side} = \ttt{Right}`, :math:`\mat(A)` is n-by-n and +/// the operation is :math:`\mat(Y) = \alpha B A + \beta Y`. +/// +/// Dispatch: +/// - If RandBLAS was built with MKL and the index type matches MKL_INT, +/// try the MKL fast path (mkl_sparse_d_mm with +/// ``SPARSE_MATRIX_TYPE_SYMMETRIC`` descriptor). The MKL path covers +/// side=Left for CSR and COO; it returns false (fall back) for CSC and +/// for side=Right. +/// - Otherwise (or on MKL fall back), call the format-specific fallback +/// kernel (``csr_spsymm`` / ``csc_spsymm`` / ``coo_spsymm``). +/// @endverbatim +template +void spsymm( + blas::Layout layout, + blas::Side side, + 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 +) { + using sint_t = typename SpMat::index_t; + constexpr bool is_coo = std::is_same_v>; + constexpr bool is_csr = std::is_same_v>; + constexpr bool is_csc = std::is_same_v>; + static_assert(is_coo || is_csr || is_csc, + "RandBLAS::sparse_data::spsymm requires COO, CSR, or CSC."); + + randblas_require(A.n_rows == A.n_cols); + int64_t k = (side == blas::Side::Left) ? m : n; + randblas_require(A.n_rows == k); + +#if defined(RandBLAS_HAS_MKL) + if constexpr (sizeof(sint_t) == sizeof(MKL_INT)) { + bool handled = mkl::mkl_spsymm( + layout, side, uplo, m, n, + alpha, A, B, ldb, beta, Y, ldy + ); + if (handled) return; + } +#endif + + // Fallback path + if constexpr (is_csr) { + csr_spsymm(layout, side, uplo, m, n, alpha, A, B, ldb, beta, Y, ldy); + } else if constexpr (is_csc) { + csc_spsymm(layout, side, uplo, m, n, alpha, A, B, ldb, beta, Y, ldy); + } else if constexpr (is_coo) { + coo_spsymm(layout, side, uplo, m, n, alpha, A, B, ldb, beta, Y, ldy); + } +} + +} // end namespace RandBLAS::sparse_data + + +namespace RandBLAS { + +// ============================================================================= +/// Convenience wrapper for symmetric sparse-times-dense matmul. +/// +/// Computes :math:`\ttt{Y} = \alpha A B + \beta Y` (side=Left default), where +/// :math:`A` is the symmetric sparse matrix and only the triangle named by +/// :math:`\ttt{uplo}` is read. +template +inline void 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 +) { + RandBLAS::sparse_data::spsymm( + layout, blas::Side::Left, uplo, m, n, + alpha, A, B, ldb, beta, Y, ldy + ); +} + +// ============================================================================= +/// Convenience overload taking a Symmetric wrapper. The uplo is read +/// from the wrapper. Defaults to side=Left. +template +inline void spsymm( + blas::Layout layout, + int64_t m, int64_t n, + T alpha, + const Symmetric& A_sym, + const T* B, int64_t ldb, + T beta, + T* Y, int64_t ldy +) { + RandBLAS::sparse_data::spsymm( + layout, blas::Side::Left, A_sym.uplo, m, n, + alpha, A_sym.A, B, ldb, beta, Y, ldy + ); +} + +} // end namespace RandBLAS diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 60de660c..6d966522 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -45,6 +45,7 @@ if (GTest_FOUND) linops/test_spmm/test_spmm_csc.cc linops/test_spmm/test_spmm_csr.cc linops/test_spmm/test_spmm_coo.cc + linops/test_spsymm.cc linops/test_sketch_sparse.cc ) diff --git a/test/linops/test_spsymm.cc b/test/linops/test_spsymm.cc new file mode 100644 index 00000000..fff66fd2 --- /dev/null +++ b/test/linops/test_spsymm.cc @@ -0,0 +1,232 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#include "RandBLAS/dense_skops.hh" +#include "RandBLAS/sparse_data/coo_matrix.hh" +#include "RandBLAS/sparse_data/csr_matrix.hh" +#include "RandBLAS/sparse_data/csc_matrix.hh" +#include "RandBLAS/sparse_data/spsymm_dispatch.hh" +#include "RandBLAS/util.hh" +#include "RandBLAS/testing/comparison.hh" + +#include +#include +#include + +using namespace RandBLAS::sparse_data; +using namespace RandBLAS::sparse_data::coo; +using namespace RandBLAS::sparse_data::csr; +using namespace RandBLAS::sparse_data::csc; +using blas::Layout; +using blas::Side; +using blas::Uplo; +using RandBLAS::Axis; +using RandBLAS::DenseDist; +using RandBLAS::RNGState; +using RandBLAS::ScalarDist; + + +class TestSpsymm : public ::testing::Test { +protected: + template + static void fill_sym_dense(int64_t n, T* A, int64_t lda, uint32_t seed) { + DenseDist D(lda, lda, ScalarDist::Uniform); + RandBLAS::fill_dense_unpacked(Layout::ColMajor, D, n, n, 0, 0, A, RNGState(seed)); + RandBLAS::symmetrize(Layout::ColMajor, Uplo::Upper, n, A, lda); + } + + template + static void zero_other_triangle(int64_t n, T* A, int64_t lda, Uplo uplo) { + // Zero out the strict triangle NOT named by uplo, leaving only the + // structurally-stored side + diagonal. + if (uplo == Uplo::Upper) { + for (int64_t j = 0; j < n; ++j) + for (int64_t i = j + 1; i < n; ++i) + A[i + j * lda] = T(0); + } else { + for (int64_t j = 0; j < n; ++j) + for (int64_t i = 0; i < j; ++i) + A[i + j * lda] = T(0); + } + } + + template + static void dense_to_sparse_format(Layout layout, T* dense_buf, T abs_tol, SpMat& sp) { + using sint_t = typename SpMat::index_t; + if constexpr (std::is_same_v>) { + dense_to_coo(layout, dense_buf, abs_tol, sp); + } else if constexpr (std::is_same_v>) { + dense_to_csr(layout, dense_buf, abs_tol, sp); + } else if constexpr (std::is_same_v>) { + dense_to_csc(layout, dense_buf, abs_tol, sp); + } else { + static_assert(sizeof(SpMat) == 0, "Unsupported sparse format."); + } + } + + template + static void run_case( + Layout layout, Side side, Uplo uplo, + int64_t n_A, int64_t d, + T alpha, T beta, + uint32_t seed_A, uint32_t seed_B + ) { + // For side=Left: Y = alpha*A*B + beta*Y, A is n_A x n_A, B and Y are n_A x d. + // For side=Right: Y = alpha*B*A + beta*Y, B and Y are d x n_A, A is n_A x n_A. + int64_t m_BY, n_BY; + if (side == Side::Left) { m_BY = n_A; n_BY = d; } + else { m_BY = d; n_BY = n_A; } + + // Build dense symmetric A. + int64_t lda = n_A; + std::vector A_full(lda * n_A, T(0)); + fill_sym_dense(n_A, A_full.data(), lda, seed_A); + + // Build sparse A: take A_full, zero out the non-named triangle, then convert. + std::vector A_triangle(A_full); + zero_other_triangle(n_A, A_triangle.data(), lda, uplo); + SpMat A_sparse(n_A, n_A); + dense_to_sparse_format(Layout::ColMajor, A_triangle.data(), T(0), A_sparse); + + // Build random B and an initial Y (for beta != 0 to be non-trivial). + int64_t ldb = (layout == Layout::ColMajor) ? m_BY : n_BY; + int64_t ldy = ldb; + std::vector B(m_BY * n_BY); + DenseDist DB(m_BY, n_BY, ScalarDist::Uniform); + RandBLAS::fill_dense_unpacked(layout, DB, m_BY, n_BY, 0, 0, B.data(), RNGState(seed_B)); + + std::vector Y_actual(m_BY * n_BY); + RandBLAS::fill_dense_unpacked(layout, DB, m_BY, n_BY, 0, 0, Y_actual.data(), RNGState(seed_B + 7)); + std::vector Y_expect = Y_actual; + + // Reference using dense blas::symm on the fully-populated A_full + // (both triangles match because A_full is symmetrized; choice of uplo + // for the reference doesn't matter, but we pass `uplo` for consistency). + blas::symm(layout, side, uplo, m_BY, n_BY, + alpha, A_full.data(), lda, B.data(), ldb, + beta, Y_expect.data(), ldy); + + // Under test: spsymm on the one-triangle sparse A. + RandBLAS::sparse_data::spsymm(layout, side, uplo, m_BY, n_BY, + alpha, A_sparse, B.data(), ldb, + beta, Y_actual.data(), ldy); + + // Tolerance: the dense reference (blas::symm) and the sparse path + // accumulate in different orders, so we get a few ULPs of accumulation + // divergence. Use 100*eps absolute tolerance to absorb that. + T atol = T(100) * std::numeric_limits::epsilon(); + T rtol = T(10) * std::numeric_limits::epsilon(); + auto msg = RandBLAS::testing::matrices_approx_equal( + layout, blas::Op::NoTrans, m_BY, n_BY, + Y_actual.data(), ldy, Y_expect.data(), ldy, + __PRETTY_FUNCTION__, __FILE__, __LINE__, atol, rtol + ); + if (!msg.empty()) FAIL() << msg; + } + + template + static void sweep_layout_uplo(Side side, int64_t n_A, int64_t d, double alpha, double beta) { + using T = typename SpMat::scalar_t; + for (auto layout : {Layout::ColMajor, Layout::RowMajor}) { + for (auto uplo : {Uplo::Upper, Uplo::Lower}) { + SCOPED_TRACE(testing::Message() << + "layout=" << (layout == Layout::ColMajor ? "Col" : "Row") + << " uplo=" << (uplo == Uplo::Upper ? "U" : "L")); + run_case(layout, side, uplo, n_A, d, T(alpha), T(beta), 0, 1); + } + } + } +}; + + +// 24-cell coverage: {COO, CSR, CSC} x {ColMajor, RowMajor} x {Upper, Lower} x {Left, Right}. +// Each TEST_F sweeps the (layout, uplo) plane internally for one (format, side) combination. + +TEST_F(TestSpsymm, CSR_Left) { sweep_layout_uplo>(Side::Left, 10, 4, 1.5, -0.5); } +TEST_F(TestSpsymm, CSR_Right) { sweep_layout_uplo>(Side::Right, 10, 4, 1.5, -0.5); } +TEST_F(TestSpsymm, CSC_Left) { sweep_layout_uplo>(Side::Left, 10, 4, 1.5, -0.5); } +TEST_F(TestSpsymm, CSC_Right) { sweep_layout_uplo>(Side::Right, 10, 4, 1.5, -0.5); } +TEST_F(TestSpsymm, COO_Left) { sweep_layout_uplo>(Side::Left, 10, 4, 1.5, -0.5); } +TEST_F(TestSpsymm, COO_Right) { sweep_layout_uplo>(Side::Right, 10, 4, 1.5, -0.5); } + +// Float coverage for one representative format +TEST_F(TestSpsymm, CSR_Left_Float) { sweep_layout_uplo>(Side::Left, 10, 4, 1.5, -0.5); } + +// beta=0 (init-from-zero) edge case +TEST_F(TestSpsymm, CSR_Left_Beta0) { + for (auto layout : {Layout::ColMajor, Layout::RowMajor}) + for (auto uplo : {Uplo::Upper, Uplo::Lower}) + run_case>(layout, Side::Left, uplo, 10, 4, 1.0, 0.0, 0, 2); +} + +// alpha=0 (only beta-scaling) edge case +TEST_F(TestSpsymm, CSR_Left_Alpha0) { + run_case>(Layout::ColMajor, Side::Left, Uplo::Upper, 10, 4, 0.0, 0.5, 0, 3); +} + +// Symmetric wrapper routing: covers the public RandBLAS::spsymm(layout, Symmetric, ...) overload +TEST_F(TestSpsymm, SymmetricWrapper) { + using SpMat = CSRMatrix; + int64_t n_A = 8; + int64_t d = 3; + int64_t lda = n_A; + std::vector A_full(lda * n_A, 0.0); + fill_sym_dense(n_A, A_full.data(), lda, 0); + std::vector A_tri(A_full); + zero_other_triangle(n_A, A_tri.data(), lda, Uplo::Upper); + SpMat A_sparse(n_A, n_A); + dense_to_csr(Layout::ColMajor, A_tri.data(), 0.0, A_sparse); + + int64_t m_BY = n_A, n_BY = d; + int64_t ldb = m_BY, ldy = m_BY; + std::vector B(m_BY * n_BY); + DenseDist DB(m_BY, n_BY, ScalarDist::Uniform); + RandBLAS::fill_dense_unpacked(Layout::ColMajor, DB, m_BY, n_BY, 0, 0, B.data(), RNGState(11)); + std::vector Y_actual(m_BY * n_BY, 0.0); + std::vector Y_expect = Y_actual; + + double alpha = 1.0, beta = 0.0; + blas::symm(Layout::ColMajor, Side::Left, Uplo::Upper, m_BY, n_BY, + alpha, A_full.data(), lda, B.data(), ldb, + beta, Y_expect.data(), ldy); + + // Route via wrapper overload at the public RandBLAS:: namespace + auto A_sym = RandBLAS::as_symmetric(A_sparse, Uplo::Upper); + RandBLAS::spsymm(Layout::ColMajor, m_BY, n_BY, alpha, A_sym, + B.data(), ldb, beta, Y_actual.data(), ldy); + + double atol = 100.0 * std::numeric_limits::epsilon(); + double rtol = 10.0 * std::numeric_limits::epsilon(); + auto msg = RandBLAS::testing::matrices_approx_equal( + Layout::ColMajor, blas::Op::NoTrans, m_BY, n_BY, + Y_actual.data(), ldy, Y_expect.data(), ldy, + __PRETTY_FUNCTION__, __FILE__, __LINE__, atol, rtol + ); + if (!msg.empty()) FAIL() << msg; +} From fce1e5a63a7f0c24908f9f14933672287d64010d Mon Sep 17 00:00:00 2001 From: mmelnich Date: Wed, 13 May 2026 19:35:24 -0400 Subject: [PATCH 04/19] Add Case D stub (sparse-symm x sparse -> dense) and document the 4-case 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 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). --- RandBLAS/sparse_data/DevNotes.md | 67 +++++++++++++++++++++++++ RandBLAS/sparse_data/spsymm_dispatch.hh | 56 +++++++++++++++++++++ test/linops/test_spsymm.cc | 16 ++++++ 3 files changed, 139 insertions(+) diff --git a/RandBLAS/sparse_data/DevNotes.md b/RandBLAS/sparse_data/DevNotes.md index 24d52336..2c3b017f 100644 --- a/RandBLAS/sparse_data/DevNotes.md +++ b/RandBLAS/sparse_data/DevNotes.md @@ -54,6 +54,73 @@ either ``lskges`` or ``rskges``. Here's what would happen after we entered one o (if inside ``lskges``) or ``right_spmm`` (if inside ``rskges``). +## SYMM-shaped kernels (spsymm) + +RandBLAS exposes a SYMM-style API for sparse symmetric matrices via the +``spsymm`` family. The design covers four cases based on the structure of the +two operands (``A`` symmetric vs. the second factor ``B``): + +| Tag | Operation | A storage | B storage | Status this PR | +|-----|---------------------------------|---------------------|-----------|-------------------------------------------------| +| A | dense-symm × dense | dense, one triangle | dense | Implemented via ``blas::symm`` in ``sksy.hh`` | +| B | dense-symm × sparse | dense, one triangle | sparse | Stub in ``sksy.hh`` (``sketch_symmetric`` on SparseSkOp). Throws ``RandBLAS::Error``. | +| 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 | Stub in ``spsymm_dispatch.hh`` (overload taking two ``SparseMatrix`` args). Throws ``RandBLAS::Error``. | + +### MKL availability + +| Tag | MKL native? | Notes | +|-----|----------------|-----------------------------------------------------------------------------------------------------------------| +| A | No (BLAS++) | ``blas::symm`` directly. | +| B | No | The transpose trick puts the sparse op on the left of ``mkl_sparse_d_mm``, but the dense A has no ``matrix_descr``, so MKL can't be told A is symmetric. | +| C | Yes (mostly) | ``mkl_sparse_d_mm`` with ``descr.type = SPARSE_MATRIX_TYPE_SYMMETRIC``. RandBLAS falls back to a hand kernel for side=Right (MKL has no Side parameter) and for CSC (``mkl_sparse_d_mm`` returns NOT_SUPPORTED on CSC). | +| D | Partial | ``mkl_sparse_sp2m`` accepts a symmetric descriptor on A but writes sparse output; densifying afterward is a workable composition fallback but not a single-call kernel. | + +### Case C dispatch (``spsymm_dispatch.hh``) + +``RandBLAS::sparse_data::spsymm(layout, side, uplo, m, n, alpha, A, B, ldb, beta, Y, ldy)`` +dispatches as follows: + +1. Validate: A is square (``A.n_rows == A.n_cols``), and matches the side + convention (``A.n_rows == m`` for side=Left, ``A.n_rows == n`` for side=Right). +2. If RandBLAS was built with MKL and the index width matches ``MKL_INT``, + try ``mkl::mkl_spsymm``. It applies the ``SPARSE_MATRIX_TYPE_SYMMETRIC`` + descriptor and calls ``mkl_sparse_d_mm`` directly. Returns false for + side=Right and CSC; control falls through to step 3 in either case. +3. Format-specific fallback: ``csr_spsymm`` / ``csc_spsymm`` / ``coo_spsymm``. + Each iterates the named triangle once. For each stored entry ``A(i,j) = v``, + it emits one ``blas::axpy`` for the structural location and a second one + for the implied symmetric counterpart (when ``i != j``). Diagonal entries + contribute once. Entries outside the named triangle are silently skipped, + so a caller that mistakenly stored both triangles still gets the correct + answer (the kernel just behaves as if the "extra" entries were absent). + +A shared ``internal::apply_beta_scale`` helper (defined in +``csr_spsymm_impl.hh``, re-included by the other two format files) handles +the ``Y <- beta * Y`` pass on entry. + +The public-facing wrappers in the top-level ``RandBLAS::`` namespace are: + + - ``spsymm(layout, uplo, m, n, alpha, A, B, ldb, beta, Y, ldy)`` -- + convenience for side=Left. + - ``spsymm(layout, m, n, alpha, Symmetric A_sym, B, ldb, beta, Y, ldy)`` + -- routes via the ``Symmetric`` carrier so the uplo annotation + travels with the matrix. + +### Cases B and D: why stub-only + +No portable kernel for "dense-symm × sparse" or "sparse-symm × sparse → dense" +exists in MKL, Ginkgo, the SparseBLAS reference implementation, the C++ +Sparse BLAS proposal (arXiv:2411.13259), or MAGMA-sparse (surveyed 2026-05). +Hand-rolling is feasible but non-trivial: case B's access pattern is "for +each nonzero of B, gather a column of A from the stored triangle", which is +awkward to vectorize; case D layers a symmetric-triangle filter on top of +an spgemm-into-dense pattern. + +The stubs exist so the API surface is locked: future PRs can fill in the +bodies without breaking source compatibility for downstream callers. The +stub message points back to ``project-plans/randblas-symm-plan.md``. + ## Sketching sparse data with dense operators If we call ``sketch_sparse`` with a DenseSkOp, ``S``, and a sparse matrix, ``A``, then we'll get routed to either diff --git a/RandBLAS/sparse_data/spsymm_dispatch.hh b/RandBLAS/sparse_data/spsymm_dispatch.hh index b12b3f65..4c5012f2 100644 --- a/RandBLAS/sparse_data/spsymm_dispatch.hh +++ b/RandBLAS/sparse_data/spsymm_dispatch.hh @@ -124,6 +124,62 @@ void spsymm( } } +// ============================================================================= +/// Case D stub: sparse-symmetric A times sparse B, dense output. +/// +/// @verbatim embed:rst:leading-slashes +/// Not implemented in this PR. The function signature is reserved here so that +/// future PRs can fill in the body without breaking source compatibility for +/// downstream users. Calling this overload triggers ``RandBLAS::Error`` via +/// ``randblas_require(false, ...)``. +/// +/// Why no implementation? No portable kernel for "sparse-symmetric A times +/// sparse B into a dense result" exists in MKL, Ginkgo, the SparseBLAS +/// reference (``SparseBLAS/spblas-reference``), or MAGMA-sparse (all surveyed +/// 2026-05). The C++ Sparse BLAS standardization proposal (arXiv:2411.13259) +/// does not specify a SYMM-shaped sparse op either. Hand-rolling the kernel is +/// feasible but ~1-2 weeks of work, deferred outside the current PR. +/// +/// Composition fallbacks callers can use today: +/// +/// - Densify B and call the Case-C spsymm (above). Exploits A's symmetry +/// but loses B's sparsity benefit. +/// - With Intel MKL: ``mkl_sparse_sp2m`` accepts a symmetric ``matrix_descr`` +/// on A and produces a sparse C, which can be densified afterward. +/// Exploits both structures but pays an intermediate sparse-C plus a +/// dense-fill step. +/// +/// See project-plans/randblas-symm-plan.md for the full design context. +/// @endverbatim +template +void spsymm( + blas::Layout layout, + blas::Side side, + blas::Uplo uplo, + int64_t m, int64_t n, + T alpha, + const SpMatA& A, + const SpMatB& B, + T beta, + T* Y, int64_t ldy +) { + (void) layout; (void) side; (void) uplo; + (void) m; (void) n; (void) alpha; + (void) A; (void) B; (void) beta; + (void) Y; (void) ldy; + randblas_require( + false && + "RandBLAS::sparse_data::spsymm(..., const SpMatA& A, const SpMatB& B, ...) " + "(Case D: sparse-symmetric A times sparse B into dense Y) is not " + "implemented. No portable reference kernel exists in MKL, Ginkgo, " + "spblas-reference, or MAGMA-sparse (as of 2026-05). Composition " + "fallbacks: (a) densify B and call the Case-C spsymm, or (b) call " + "mkl_sparse_sp2m with a symmetric descriptor on A and densify the " + "resulting sparse C. See project-plans/randblas-symm-plan.md." + ); +} + } // end namespace RandBLAS::sparse_data diff --git a/test/linops/test_spsymm.cc b/test/linops/test_spsymm.cc index fff66fd2..d3df6ad4 100644 --- a/test/linops/test_spsymm.cc +++ b/test/linops/test_spsymm.cc @@ -191,6 +191,22 @@ TEST_F(TestSpsymm, CSR_Left_Alpha0) { } // Symmetric wrapper routing: covers the public RandBLAS::spsymm(layout, Symmetric, ...) overload +// Case D stub: sparse-symmetric A times sparse B must throw RandBLAS::Error. +// The API surface is reserved; the body is "not implemented" pending a future PR. +TEST_F(TestSpsymm, CaseD_SparseSparseThrows) { + CSRMatrix A_sparse(4, 4); + CSRMatrix B_sparse(4, 3); + std::vector Y(4 * 3, 0.0); + EXPECT_THROW({ + RandBLAS::sparse_data::spsymm( + Layout::ColMajor, Side::Left, Uplo::Upper, + 4, 3, 1.0, + A_sparse, B_sparse, + 0.0, Y.data(), 4 + ); + }, RandBLAS::Error); +} + TEST_F(TestSpsymm, SymmetricWrapper) { using SpMat = CSRMatrix; int64_t n_A = 8; From 512fa11c14a538e42fd250ca7142235a894a7b4e Mon Sep 17 00:00:00 2001 From: mmelnich Date: Wed, 13 May 2026 19:39:00 -0400 Subject: [PATCH 05/19] RTD: document the new sketch_symmetric uplo signature and add spsymm 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-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. --- rtd/source/FAQ.rst | 11 +++++--- rtd/source/api_reference/sketch_dense.rst | 14 ++++++++--- rtd/source/api_reference/sketch_sparse.rst | 29 ++++++++++++++++++++++ 3 files changed, 46 insertions(+), 8 deletions(-) diff --git a/rtd/source/FAQ.rst b/rtd/source/FAQ.rst index ea9a0e1c..3172781a 100644 --- a/rtd/source/FAQ.rst +++ b/rtd/source/FAQ.rst @@ -110,10 +110,13 @@ No support for negative values of "incx" or "incy" in sketch_vector. to extend our SPMV kernels to support that, then we'd happily accept such a contribution. (It shouldn't be hard! We just haven't gotten around to this.) -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. +SparseSkOp is not yet supported in sketch_symmetric. + ``sketch_symmetric`` now takes a ``blas::Uplo`` and dispatches to ``blas::symm`` for ``DenseSkOp``, exploiting symmetry of :math:`A`. The ``SparseSkOp`` branch throws ``RandBLAS::Error`` --- this is Case B in the SYMM-kernels plan + (``RandBLAS/sparse_data/DevNotes.md``). No portable kernel for "dense-symm matrix times sparse operator" exists in MKL, Ginkgo, the SparseBLAS reference, or MAGMA-sparse, so we deferred the hand-roll to a future PR. + As a composition fallback, ``DenseSkOp(densify(sparse_op))`` can be passed to ``sketch_symmetric`` instead --- it loses the operator-side sparsity benefit but exploits A's symmetry. + +Layout-mismatched ``DenseSkOp`` in ``sketch_symmetric`` falls back to GEMM. + When the ``DenseSkOp``'s storage layout differs from the caller's ``layout`` parameter, ``sketch_symmetric`` falls back to ``blas::gemm`` with the transpose flag --- ``blas::symm`` has no on-the-fly transpose flag for the dense operand. The layout-matched case still gets the SYMM speedup. Language interoperability diff --git a/rtd/source/api_reference/sketch_dense.rst b/rtd/source/api_reference/sketch_dense.rst index 7dd0a571..88ae0611 100644 --- a/rtd/source/api_reference/sketch_dense.rst +++ b/rtd/source/api_reference/sketch_dense.rst @@ -67,18 +67,24 @@ Analogs to GEMM Analogs to SYMM --------------- +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. Currently only ``DenseSkOp`` is supported; +calling these with a ``SparseSkOp`` throws ``RandBLAS::Error`` (Case B of the +SYMM-kernels plan; see ``RandBLAS/sparse_data/DevNotes.md``). + .. dropdown:: :math:`\mtxB = \alpha \cdot \mtxS \cdot \mtxA + \beta \cdot \mtxB` :animate: fade-in-slide-down :color: light - .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, T alpha, const SKOP &S, const T *A, int64_t lda, T beta, T *B, int64_t ldb, T sym_check_tol = 0) + .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, blas::Uplo uplo, T alpha, const SKOP &S, const T *A, int64_t lda, T beta, T *B, int64_t ldb) :project: RandBLAS .. dropdown:: :math:`\mtxB = \alpha \cdot \mtxA \cdot \mtxS + \beta \cdot \mtxB` :animate: fade-in-slide-down :color: light - .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, T alpha, const T *A, int64_t lda, const SKOP &S, T beta, T *B, int64_t ldb, T sym_check_tol = 0) + .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, blas::Uplo uplo, T alpha, const T *A, int64_t lda, const SKOP &S, T beta, T *B, int64_t ldb) :project: RandBLAS @@ -86,10 +92,10 @@ Analogs to SYMM :animate: fade-in-slide-down :color: light - .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, int64_t d, int64_t n, T alpha, const SKOP &S, int64_t ro_s, int64_t co_s, const T *A, int64_t lda, T beta, T *B, int64_t ldb, T sym_check_tol = 0) + .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, blas::Uplo uplo, int64_t d, int64_t n, T alpha, const SKOP &S, int64_t ro_s, int64_t co_s, const T *A, int64_t lda, T beta, T *B, int64_t ldb) :project: RandBLAS - .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, int64_t n, int64_t d, T alpha, const T *A, int64_t lda, const SKOP &S, int64_t ro_s, int64_t co_s, T beta, T *B, int64_t ldb, T sym_check_tol = 0) + .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, blas::Uplo uplo, int64_t n, int64_t d, T alpha, const T *A, int64_t lda, const SKOP &S, int64_t ro_s, int64_t co_s, T beta, T *B, int64_t ldb) :project: RandBLAS diff --git a/rtd/source/api_reference/sketch_sparse.rst b/rtd/source/api_reference/sketch_sparse.rst index 98a14df9..1b56f55a 100644 --- a/rtd/source/api_reference/sketch_sparse.rst +++ b/rtd/source/api_reference/sketch_sparse.rst @@ -114,6 +114,35 @@ Deterministic operations This function requires Intel MKL and only supports single and double precision (``float`` and ``double``), in contrast to other RandBLAS kernels that work with any scalar type. +.. 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 &A_sym, const T *B, int64_t ldb, T beta, T *Y, int64_t ldy) + :project: RandBLAS + + Only the triangle named by ``uplo`` is read from :math:`\mtxA`. With Intel MKL, + the call dispatches to ``mkl_sparse_d_mm`` with a symmetric ``matrix_descr`` + for the fast path; otherwise it uses a hand-rolled per-format kernel + (CSR / CSC / COO). See ``RandBLAS/sparse_data/DevNotes.md`` for the dispatch + flow and the broader 4-case design that ``spsymm`` belongs to. + + .. note:: + + The ``Symmetric`` carrier (a non-owning wrapper holding a + ``const SpMat &`` and a ``blas::Uplo``) is the type-safety hook that + prevents a symmetric sparse matrix from being accidentally passed to + the general ``spmm`` / ``spgemm`` routines. + + Companion stubs exist for the cases where the second factor is also + sparse (Case D) or where the symmetric factor is dense and the second + is sparse (Case B, via ``sketch_symmetric``). They throw + ``RandBLAS::Error`` with a pointer to the design plan; no portable + kernel for those shapes exists in current sparse-BLAS libraries. + .. dropdown:: :math:`\mtxB = \alpha \cdot \op(\mtxA)^{-1} \cdot \mtxB,` with sparse triangular :math:`\mtxA` :animate: fade-in-slide-down :color: light From c6638ce557c6e30b6fef4810b98630d71c034cb6 Mon Sep 17 00:00:00 2001 From: mmelnich Date: Thu, 14 May 2026 12:01:55 -0400 Subject: [PATCH 06/19] Extract lascl helper to RandBLAS::util; replace hand-rolled beta-scale 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 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. --- RandBLAS/sparse_data/coo_spsymm_impl.hh | 4 +-- RandBLAS/sparse_data/csc_spsymm_impl.hh | 4 +-- RandBLAS/sparse_data/csr_spsymm_impl.hh | 33 ++------------------ RandBLAS/sparse_data/mkl_spmm_impl.hh | 28 +++++------------ RandBLAS/util.hh | 40 +++++++++++++++++++++++++ 5 files changed, 53 insertions(+), 56 deletions(-) diff --git a/RandBLAS/sparse_data/coo_spsymm_impl.hh b/RandBLAS/sparse_data/coo_spsymm_impl.hh index 981262a8..4664f402 100644 --- a/RandBLAS/sparse_data/coo_spsymm_impl.hh +++ b/RandBLAS/sparse_data/coo_spsymm_impl.hh @@ -30,9 +30,9 @@ #pragma once #include "RandBLAS/exceptions.hh" +#include "RandBLAS/util.hh" #include "RandBLAS/sparse_data/base.hh" #include "RandBLAS/sparse_data/coo_matrix.hh" -#include "RandBLAS/sparse_data/csr_spsymm_impl.hh" // for internal::apply_beta_scale #include namespace RandBLAS::sparse_data { @@ -60,7 +60,7 @@ void coo_spsymm( int64_t k = (side == blas::Side::Left) ? m : n; randblas_require(A.n_rows == k); - internal::apply_beta_scale(layout, m, n, beta, Y, ldy); + RandBLAS::util::lascl(layout, m, n, beta, Y, ldy); if (alpha == T(0)) return; const bool col_major = (layout == blas::Layout::ColMajor); diff --git a/RandBLAS/sparse_data/csc_spsymm_impl.hh b/RandBLAS/sparse_data/csc_spsymm_impl.hh index e860df6a..e41a2acf 100644 --- a/RandBLAS/sparse_data/csc_spsymm_impl.hh +++ b/RandBLAS/sparse_data/csc_spsymm_impl.hh @@ -30,9 +30,9 @@ #pragma once #include "RandBLAS/exceptions.hh" +#include "RandBLAS/util.hh" #include "RandBLAS/sparse_data/base.hh" #include "RandBLAS/sparse_data/csc_matrix.hh" -#include "RandBLAS/sparse_data/csr_spsymm_impl.hh" // for internal::apply_beta_scale #include namespace RandBLAS::sparse_data { @@ -60,7 +60,7 @@ void csc_spsymm( int64_t k = (side == blas::Side::Left) ? m : n; randblas_require(A.n_rows == k); - internal::apply_beta_scale(layout, m, n, beta, Y, ldy); + RandBLAS::util::lascl(layout, m, n, beta, Y, ldy); if (alpha == T(0)) return; const bool col_major = (layout == blas::Layout::ColMajor); diff --git a/RandBLAS/sparse_data/csr_spsymm_impl.hh b/RandBLAS/sparse_data/csr_spsymm_impl.hh index 67fb2562..86835ce2 100644 --- a/RandBLAS/sparse_data/csr_spsymm_impl.hh +++ b/RandBLAS/sparse_data/csr_spsymm_impl.hh @@ -30,42 +30,13 @@ #pragma once #include "RandBLAS/exceptions.hh" +#include "RandBLAS/util.hh" #include "RandBLAS/sparse_data/base.hh" #include "RandBLAS/sparse_data/csr_matrix.hh" #include -#include namespace RandBLAS::sparse_data { -namespace internal { - -// Apply Y <- beta * Y for an m-by-n dense matrix with leading dimension ldy. -// Used by the spsymm fallback kernels to fold beta-scaling out of the -// nonzero-traversal loop. -template -inline void apply_beta_scale(blas::Layout layout, int64_t m, int64_t n, T beta, T* Y, int64_t ldy) { - if (beta == T(1)) return; - if (beta == T(0)) { - if (layout == blas::Layout::ColMajor) { - for (int64_t j = 0; j < n; ++j) - std::fill(Y + j*ldy, Y + j*ldy + m, T(0)); - } else { - for (int64_t i = 0; i < m; ++i) - std::fill(Y + i*ldy, Y + i*ldy + n, T(0)); - } - return; - } - if (layout == blas::Layout::ColMajor) { - for (int64_t j = 0; j < n; ++j) - blas::scal(m, beta, Y + j*ldy, 1); - } else { - for (int64_t i = 0; i < m; ++i) - blas::scal(n, beta, Y + i*ldy, 1); - } -} - -} // namespace internal - // ============================================================================= /// CSR fallback for the symmetric sparse-times-dense kernel. @@ -96,7 +67,7 @@ void csr_spsymm( int64_t k = (side == blas::Side::Left) ? m : n; randblas_require(A.n_rows == k); - internal::apply_beta_scale(layout, m, n, beta, Y, ldy); + RandBLAS::util::lascl(layout, m, n, beta, Y, ldy); if (alpha == T(0)) return; const bool col_major = (layout == blas::Layout::ColMajor); diff --git a/RandBLAS/sparse_data/mkl_spmm_impl.hh b/RandBLAS/sparse_data/mkl_spmm_impl.hh index 6d2af501..f027b633 100644 --- a/RandBLAS/sparse_data/mkl_spmm_impl.hh +++ b/RandBLAS/sparse_data/mkl_spmm_impl.hh @@ -42,6 +42,7 @@ #include #include "RandBLAS/exceptions.hh" +#include "RandBLAS/util.hh" #include "RandBLAS/sparse_data/base.hh" #include "RandBLAS/sparse_data/coo_matrix.hh" #include "RandBLAS/sparse_data/csr_matrix.hh" @@ -369,20 +370,8 @@ void mkl_spgemm_to_dense( // 3. C = alpha * temp + C (if needed) if (alpha == (T)0) { - // Just scale C by beta - if (beta == (T)0) { - int64_t total = (layout == blas::Layout::ColMajor) ? ldc * n : ldc * m; - std::fill(C, C + total, (T)0); - } else if (beta != (T)1) { - // Scale each column/row of C - if (layout == blas::Layout::ColMajor) { - 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); - } - } + // Just scale C by beta. + RandBLAS::util::lascl(layout, m, n, beta, C, ldc); return; } @@ -412,17 +401,14 @@ void mkl_spgemm_to_dense( check_mkl_status(status, "mkl_sparse_spmmd"); if (!direct_write) { - // C = alpha * target + beta * C + // C = alpha * target + beta * C: hoist the scale, then per-vector axpy. + RandBLAS::util::lascl(layout, m, n, beta, C, ldc); if (layout == blas::Layout::ColMajor) { - for (int64_t j = 0; j < n; ++j) { - blas::scal(m, beta, &C[j * ldc], 1); + for (int64_t j = 0; j < n; ++j) blas::axpy(m, alpha, &target[j * ldc], 1, &C[j * ldc], 1); - } } else { - for (int64_t i = 0; i < m; ++i) { - blas::scal(n, beta, &C[i * ldc], 1); + for (int64_t i = 0; i < m; ++i) blas::axpy(n, alpha, &target[i * ldc], 1, &C[i * ldc], 1); - } } } } diff --git a/RandBLAS/util.hh b/RandBLAS/util.hh index 5ee923b2..3834018e 100644 --- a/RandBLAS/util.hh +++ b/RandBLAS/util.hh @@ -61,6 +61,46 @@ void safe_scal(int64_t n, T a, T* x) { } } +// ============================================================================= +/// \fn lascl(blas::Layout layout, int64_t m, int64_t n, T alpha, T* A, int64_t lda) +/// @verbatim embed:rst:leading-slashes +/// In-place scale a dense :math:`m \times n` matrix: +/// +/// .. math:: +/// A \leftarrow \alpha \cdot A. +/// +/// Named after LAPACK's ``?lascl`` for the convention; the signature is +/// simpler than LAPACK's (no ``CFROM`` / ``CTO`` overflow protection, no +/// matrix-type flag --- general dense only). +/// +/// Fast paths: +/// +/// - :math:`\alpha = 1`: returns immediately. +/// - :math:`\alpha = 0`: ``std::fill``-based zero out. +/// - otherwise: per-column (ColMajor) or per-row (RowMajor) ``blas::scal``. +/// @endverbatim +template +inline void lascl(blas::Layout layout, int64_t m, int64_t n, T alpha, T* A, int64_t lda) { + if (alpha == T(1)) return; + if (alpha == T(0)) { + if (layout == blas::Layout::ColMajor) { + for (int64_t j = 0; j < n; ++j) + std::fill(A + j * lda, A + j * lda + m, T(0)); + } else { + for (int64_t i = 0; i < m; ++i) + std::fill(A + i * lda, A + i * lda + n, T(0)); + } + return; + } + if (layout == blas::Layout::ColMajor) { + for (int64_t j = 0; j < n; ++j) + blas::scal(m, alpha, A + j * lda, 1); + } else { + for (int64_t i = 0; i < m; ++i) + blas::scal(n, alpha, A + i * lda, 1); + } +} + template void omatcopy(int64_t m, int64_t n, const T* A, int64_t irs_a, int64_t ics_a, T* B, int64_t irs_b, int64_t ics_b) { // TODO: From 09b1dc1851a1437175041c6a2c64b4fe83270ffd Mon Sep 17 00:00:00 2001 From: mmelnich Date: Thu, 14 May 2026 12:11:19 -0400 Subject: [PATCH 07/19] Extract spsymm_scatter helpers shared across CSR / CSC / COO fallbacks 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 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 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). --- RandBLAS/sparse_data/coo_spsymm_impl.hh | 29 ++------- RandBLAS/sparse_data/csc_spsymm_impl.hh | 23 +------ RandBLAS/sparse_data/csr_spsymm_impl.hh | 29 +-------- RandBLAS/sparse_data/spsymm_internal.hh | 87 +++++++++++++++++++++++++ 4 files changed, 98 insertions(+), 70 deletions(-) create mode 100644 RandBLAS/sparse_data/spsymm_internal.hh diff --git a/RandBLAS/sparse_data/coo_spsymm_impl.hh b/RandBLAS/sparse_data/coo_spsymm_impl.hh index 4664f402..31f5a7a5 100644 --- a/RandBLAS/sparse_data/coo_spsymm_impl.hh +++ b/RandBLAS/sparse_data/coo_spsymm_impl.hh @@ -33,6 +33,7 @@ #include "RandBLAS/util.hh" #include "RandBLAS/sparse_data/base.hh" #include "RandBLAS/sparse_data/coo_matrix.hh" +#include "RandBLAS/sparse_data/spsymm_internal.hh" #include namespace RandBLAS::sparse_data { @@ -63,36 +64,16 @@ void coo_spsymm( RandBLAS::util::lascl(layout, m, n, beta, Y, ldy); if (alpha == T(0)) return; - const bool col_major = (layout == blas::Layout::ColMajor); - for (int64_t p = 0; p < A.nnz; ++p) { sint_t i = A.rows[p]; sint_t j = A.cols[p]; if (uplo == blas::Uplo::Upper && j < i) continue; if (uplo == blas::Uplo::Lower && j > i) continue; T av = alpha * A.vals[p]; - - if (side == blas::Side::Left) { - if (col_major) { - blas::axpy(n, av, &B[j], ldb, &Y[i], ldy); - if (j != i) - blas::axpy(n, av, &B[i], ldb, &Y[j], ldy); - } else { - blas::axpy(n, av, &B[j*ldb], 1, &Y[i*ldy], 1); - if (j != i) - blas::axpy(n, av, &B[i*ldb], 1, &Y[j*ldy], 1); - } - } else { - if (col_major) { - blas::axpy(m, av, &B[i*ldb], 1, &Y[j*ldy], 1); - if (j != i) - blas::axpy(m, av, &B[j*ldb], 1, &Y[i*ldy], 1); - } else { - blas::axpy(m, av, &B[i], ldb, &Y[j], ldy); - if (j != i) - blas::axpy(m, av, &B[j], ldb, &Y[i], ldy); - } - } + if (side == blas::Side::Left) + internal::spsymm_scatter_left(layout, n, av, i, j, B, ldb, Y, ldy); + else + internal::spsymm_scatter_right(layout, m, av, i, j, B, ldb, Y, ldy); } } diff --git a/RandBLAS/sparse_data/csc_spsymm_impl.hh b/RandBLAS/sparse_data/csc_spsymm_impl.hh index e41a2acf..f5fef594 100644 --- a/RandBLAS/sparse_data/csc_spsymm_impl.hh +++ b/RandBLAS/sparse_data/csc_spsymm_impl.hh @@ -33,6 +33,7 @@ #include "RandBLAS/util.hh" #include "RandBLAS/sparse_data/base.hh" #include "RandBLAS/sparse_data/csc_matrix.hh" +#include "RandBLAS/sparse_data/spsymm_internal.hh" #include namespace RandBLAS::sparse_data { @@ -63,8 +64,6 @@ void csc_spsymm( RandBLAS::util::lascl(layout, m, n, beta, Y, ldy); if (alpha == T(0)) return; - const bool col_major = (layout == blas::Layout::ColMajor); - if (side == blas::Side::Left) { // Y = alpha * A * B + ... (A is m-by-m) for (int64_t j = 0; j < m; ++j) { @@ -73,15 +72,7 @@ void csc_spsymm( if (uplo == blas::Uplo::Upper && i > j) continue; if (uplo == blas::Uplo::Lower && i < j) continue; T av = alpha * A.vals[p]; - if (col_major) { - blas::axpy(n, av, &B[j], ldb, &Y[i], ldy); - if (i != j) - blas::axpy(n, av, &B[i], ldb, &Y[j], ldy); - } else { - blas::axpy(n, av, &B[j*ldb], 1, &Y[i*ldy], 1); - if (i != j) - blas::axpy(n, av, &B[i*ldb], 1, &Y[j*ldy], 1); - } + internal::spsymm_scatter_left(layout, n, av, i, j, B, ldb, Y, ldy); } } } else { @@ -92,15 +83,7 @@ void csc_spsymm( if (uplo == blas::Uplo::Upper && i > j) continue; if (uplo == blas::Uplo::Lower && i < j) continue; T av = alpha * A.vals[p]; - if (col_major) { - blas::axpy(m, av, &B[i*ldb], 1, &Y[j*ldy], 1); - if (i != j) - blas::axpy(m, av, &B[j*ldb], 1, &Y[i*ldy], 1); - } else { - blas::axpy(m, av, &B[i], ldb, &Y[j], ldy); - if (i != j) - blas::axpy(m, av, &B[j], ldb, &Y[i], ldy); - } + internal::spsymm_scatter_right(layout, m, av, i, j, B, ldb, Y, ldy); } } } diff --git a/RandBLAS/sparse_data/csr_spsymm_impl.hh b/RandBLAS/sparse_data/csr_spsymm_impl.hh index 86835ce2..f088da96 100644 --- a/RandBLAS/sparse_data/csr_spsymm_impl.hh +++ b/RandBLAS/sparse_data/csr_spsymm_impl.hh @@ -33,6 +33,7 @@ #include "RandBLAS/util.hh" #include "RandBLAS/sparse_data/base.hh" #include "RandBLAS/sparse_data/csr_matrix.hh" +#include "RandBLAS/sparse_data/spsymm_internal.hh" #include namespace RandBLAS::sparse_data { @@ -70,50 +71,26 @@ void csr_spsymm( RandBLAS::util::lascl(layout, m, n, beta, Y, ldy); if (alpha == T(0)) return; - const bool col_major = (layout == blas::Layout::ColMajor); - if (side == blas::Side::Left) { // Y = alpha * A * B + ... (A is m-by-m) - // For stored (i, j, v): - // Y[i, :] += alpha * v * B[j, :] - // if i != j: Y[j, :] += alpha * v * B[i, :] for (int64_t i = 0; i < m; ++i) { for (int64_t p = A.rowptr[i]; p < A.rowptr[i+1]; ++p) { sint_t j = A.colidxs[p]; if (uplo == blas::Uplo::Upper && j < i) continue; if (uplo == blas::Uplo::Lower && j > i) continue; T av = alpha * A.vals[p]; - if (col_major) { - blas::axpy(n, av, &B[j], ldb, &Y[i], ldy); - if (j != i) - blas::axpy(n, av, &B[i], ldb, &Y[j], ldy); - } else { - blas::axpy(n, av, &B[j*ldb], 1, &Y[i*ldy], 1); - if (j != i) - blas::axpy(n, av, &B[i*ldb], 1, &Y[j*ldy], 1); - } + internal::spsymm_scatter_left(layout, n, av, i, j, B, ldb, Y, ldy); } } } else { // Y = alpha * B * A + ... (A is n-by-n) - // For stored (i, j, v): - // Y[:, j] += alpha * v * B[:, i] - // if i != j: Y[:, i] += alpha * v * B[:, j] for (int64_t i = 0; i < n; ++i) { for (int64_t p = A.rowptr[i]; p < A.rowptr[i+1]; ++p) { sint_t j = A.colidxs[p]; if (uplo == blas::Uplo::Upper && j < i) continue; if (uplo == blas::Uplo::Lower && j > i) continue; T av = alpha * A.vals[p]; - if (col_major) { - blas::axpy(m, av, &B[i*ldb], 1, &Y[j*ldy], 1); - if (j != i) - blas::axpy(m, av, &B[j*ldb], 1, &Y[i*ldy], 1); - } else { - blas::axpy(m, av, &B[i], ldb, &Y[j], ldy); - if (j != i) - blas::axpy(m, av, &B[j], ldb, &Y[i], ldy); - } + internal::spsymm_scatter_right(layout, m, av, i, j, B, ldb, Y, ldy); } } } diff --git a/RandBLAS/sparse_data/spsymm_internal.hh b/RandBLAS/sparse_data/spsymm_internal.hh new file mode 100644 index 00000000..4aea5678 --- /dev/null +++ b/RandBLAS/sparse_data/spsymm_internal.hh @@ -0,0 +1,87 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include +#include + +namespace RandBLAS::sparse_data::internal { + +// ============================================================================= +/// Add the contribution of one stored A(i, j) = v entry (and the implied +/// symmetric A(j, i) = v when i != j) to a side=Left spsymm output: +/// +/// Y[i, :] += av * B[j, :] +/// if (i != j) Y[j, :] += av * B[i, :] +/// +/// where av = alpha * v has already been folded by the caller. The inner +/// row updates use blas::axpy with strides determined by `layout`. +/// +/// Shared by all three spsymm fallback kernels (CSR / CSC / COO) -- they +/// differ only in their outer iteration order; the per-stored-entry scatter +/// is identical. +template +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) { + if (layout == blas::Layout::ColMajor) { + blas::axpy(n, av, &B[j], ldb, &Y[i], ldy); + if (i != j) + blas::axpy(n, av, &B[i], ldb, &Y[j], ldy); + } else { + blas::axpy(n, av, &B[j*ldb], 1, &Y[i*ldy], 1); + if (i != j) + blas::axpy(n, av, &B[i*ldb], 1, &Y[j*ldy], 1); + } +} + +// ============================================================================= +/// Side=Right counterpart of spsymm_scatter_left. For stored A(i, j) = v: +/// +/// Y[:, j] += av * B[:, i] +/// if (i != j) Y[:, i] += av * B[:, j] +template +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) { + if (layout == blas::Layout::ColMajor) { + blas::axpy(m, av, &B[i*ldb], 1, &Y[j*ldy], 1); + if (i != j) + blas::axpy(m, av, &B[j*ldb], 1, &Y[i*ldy], 1); + } else { + blas::axpy(m, av, &B[i], ldb, &Y[j], ldy); + if (i != j) + blas::axpy(m, av, &B[j], ldb, &Y[i], ldy); + } +} + +} // namespace RandBLAS::sparse_data::internal From 6db925253118c983b167beba0b2f445967797c69 Mon Sep 17 00:00:00 2001 From: mmelnich Date: Thu, 14 May 2026 12:13:57 -0400 Subject: [PATCH 08/19] mkl_spsymm: handle CSC via CSC.transpose() CSR-view + flipped uplo 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. --- RandBLAS/sparse_data/mkl_spsymm_impl.hh | 26 +++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/RandBLAS/sparse_data/mkl_spsymm_impl.hh b/RandBLAS/sparse_data/mkl_spsymm_impl.hh index 260514b6..43bc504f 100644 --- a/RandBLAS/sparse_data/mkl_spsymm_impl.hh +++ b/RandBLAS/sparse_data/mkl_spsymm_impl.hh @@ -62,9 +62,16 @@ namespace RandBLAS::sparse_data::mkl { // symmetric matrix is always on the left of the dense block. The // transpose-trick to express side=Right via layout flips depends on // leading-dim assumptions that don't always hold; safer to fall back. -// - CSC format: MKL's mkl_sparse_d_mm returns NOT_SUPPORTED for CSC even -// with a symmetric descriptor (mirrors the behavior in mkl_left_spmm). // - Index type mismatched with MKL_INT. +// +// CSC handling: MKL's mkl_sparse_d_mm returns NOT_SUPPORTED for CSC even +// with a symmetric descriptor. We work around this by taking the +// CSC.transpose() view (a lightweight CSR view over the same buffers, +// since A is symmetric so A == A^T) and recursing. The triangle the user +// named in the CSC is in the *opposite* triangle of the CSR view (a CSC +// Upper entry at (i, j) with i <= j becomes a CSR-view entry at (j, i) +// with j >= i, i.e., Lower in the CSR view), so the recursive call flips +// uplo. // ============================================================================ template bool mkl_spsymm( @@ -85,8 +92,19 @@ bool mkl_spsymm( if (side != blas::Side::Left) return false; - if constexpr (is_csc) - return false; + + if constexpr (is_csc) { + // Symmetric A: A == A^T. The CSC->CSR view is lightweight (same + // buffers, reinterpreted) and gives us a CSR matrix MKL accepts; + // the structurally stored triangle moves from {Upper,Lower} to + // the opposite side in the CSR view. + auto At = A.transpose(); + blas::Uplo uplo_flipped = (uplo == blas::Uplo::Upper) + ? blas::Uplo::Lower + : blas::Uplo::Upper; + return mkl_spsymm(layout, side, uplo_flipped, m, n, + alpha, At, B, ldb, beta, Y, ldy); + } auto h = make_mkl_handle(A); From c764d31047aec8c23d3d5f8c13ee832639b5c89c Mon Sep 17 00:00:00 2001 From: mmelnich Date: Thu, 14 May 2026 12:19:40 -0400 Subject: [PATCH 09/19] sksy: extract Case-B stub-throw helper, deduplicate four stub bodies 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 [[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. --- RandBLAS/sksy.hh | 62 +++++++++++++++++++++--------------------------- 1 file changed, 27 insertions(+), 35 deletions(-) diff --git a/RandBLAS/sksy.hh b/RandBLAS/sksy.hh index 4d26da13..c26914a2 100644 --- a/RandBLAS/sksy.hh +++ b/RandBLAS/sksy.hh @@ -180,6 +180,29 @@ using namespace RandBLAS::dense; using namespace RandBLAS::sparse; +namespace detail { + +// ============================================================================= +/// Shared Case-B stub-throw helper for the four SparseSkOp-taking +/// sketch_symmetric specializations. The variadic parameter pack consumes +/// the caller's arguments so the compiler doesn't warn about unused +/// parameters in the (genuinely-unused, by design) stub bodies. +template +[[noreturn]] inline void throw_sketch_symmetric_case_b(Args&&...) { + randblas_require( + false && + "RandBLAS::sketch_symmetric with a sparse sketching operator is the " + "Case-B kernel from project-plans/randblas-symm-plan.md. Not implemented " + "in this PR --- the API signature is reserved so future PRs can fill in " + "the body without breaking source compatibility. Composition fallback: " + "densify the SkOp then call sketch_symmetric on the dense version." + ); + __builtin_unreachable(); +} + +} // namespace detail + + // MARK: SUBMAT(S) // ============================================================================= @@ -307,17 +330,7 @@ inline void sketch_symmetric( T beta, T* B, int64_t ldb ) { - (void) layout; (void) uplo; (void) d; (void) n; (void) alpha; - (void) S; (void) ro_s; (void) co_s; (void) A; (void) lda; - (void) beta; (void) B; (void) ldb; - randblas_require( - false && - "RandBLAS::sketch_symmetric with a sparse sketching operator is the " - "Case-B kernel from project-plans/randblas-symm-plan.md. Not implemented " - "in this PR --- the API signature is reserved so future PRs can fill in " - "the body without breaking source compatibility. Composition fallback: " - "densify the SkOp then call sketch_symmetric on the dense version." - ); + detail::throw_sketch_symmetric_case_b(layout, uplo, d, n, alpha, S, ro_s, co_s, A, lda, beta, B, ldb); } @@ -376,14 +389,7 @@ inline void sketch_symmetric( T beta, T* B, int64_t ldb ) { - (void) layout; (void) uplo; (void) n; (void) d; (void) alpha; - (void) S; (void) ro_s; (void) co_s; (void) A; (void) lda; - (void) beta; (void) B; (void) ldb; - randblas_require( - false && - "RandBLAS::sketch_symmetric with a sparse sketching operator is the " - "Case-B kernel from project-plans/randblas-symm-plan.md. Not implemented." - ); + detail::throw_sketch_symmetric_case_b(layout, uplo, n, d, alpha, A, lda, S, ro_s, co_s, beta, B, ldb); } @@ -441,14 +447,7 @@ inline void sketch_symmetric( T beta, T* B, int64_t ldb ) { - (void) layout; (void) uplo; (void) alpha; - (void) S; (void) A; (void) lda; - (void) beta; (void) B; (void) ldb; - randblas_require( - false && - "RandBLAS::sketch_symmetric with a sparse sketching operator is the " - "Case-B kernel from project-plans/randblas-symm-plan.md. Not implemented." - ); + detail::throw_sketch_symmetric_case_b(layout, uplo, alpha, S, A, lda, beta, B, ldb); } @@ -503,14 +502,7 @@ inline void sketch_symmetric( T beta, T* B, int64_t ldb ) { - (void) layout; (void) uplo; (void) alpha; - (void) S; (void) A; (void) lda; - (void) beta; (void) B; (void) ldb; - randblas_require( - false && - "RandBLAS::sketch_symmetric with a sparse sketching operator is the " - "Case-B kernel from project-plans/randblas-symm-plan.md. Not implemented." - ); + detail::throw_sketch_symmetric_case_b(layout, uplo, alpha, S, A, lda, beta, B, ldb); } } // end namespace RandBLAS From 3071f5ebe20cca496b38cd9685d6023861a6d0fa Mon Sep 17 00:00:00 2001 From: mmelnich Date: Thu, 14 May 2026 12:23:42 -0400 Subject: [PATCH 10/19] test_spsymm: share run_case setup with the wrapper-routing test 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 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, ...) 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 overload at the public namespace and the as_symmetric sugar. Verified: 11 / 11 TestSpsymm tests pass, including SymmetricWrapper. --- test/linops/test_spsymm.cc | 69 ++++++++++++++++---------------------- 1 file changed, 28 insertions(+), 41 deletions(-) diff --git a/test/linops/test_spsymm.cc b/test/linops/test_spsymm.cc index d3df6ad4..17f33b17 100644 --- a/test/linops/test_spsymm.cc +++ b/test/linops/test_spsymm.cc @@ -95,7 +95,8 @@ class TestSpsymm : public ::testing::Test { Layout layout, Side side, Uplo uplo, int64_t n_A, int64_t d, T alpha, T beta, - uint32_t seed_A, uint32_t seed_B + uint32_t seed_A, uint32_t seed_B, + bool route_via_wrapper = false ) { // For side=Left: Y = alpha*A*B + beta*Y, A is n_A x n_A, B and Y are n_A x d. // For side=Right: Y = alpha*B*A + beta*Y, B and Y are d x n_A, A is n_A x n_A. @@ -132,10 +133,22 @@ class TestSpsymm : public ::testing::Test { alpha, A_full.data(), lda, B.data(), ldb, beta, Y_expect.data(), ldy); - // Under test: spsymm on the one-triangle sparse A. - RandBLAS::sparse_data::spsymm(layout, side, uplo, m_BY, n_BY, - alpha, A_sparse, B.data(), ldb, - beta, Y_actual.data(), ldy); + // Under test: spsymm on the one-triangle sparse A. By default we + // call the low-level dispatcher directly; if `route_via_wrapper` is + // set, we go through the public RandBLAS::spsymm(Symmetric) + // overload to exercise the wrapper-routing path. side=Left only + // for the wrapper path since the public wrapper defaults side=Left. + if (route_via_wrapper) { + randblas_require(side == Side::Left); + auto A_sym = RandBLAS::as_symmetric(A_sparse, uplo); + RandBLAS::spsymm(layout, m_BY, n_BY, + alpha, A_sym, B.data(), ldb, + beta, Y_actual.data(), ldy); + } else { + RandBLAS::sparse_data::spsymm(layout, side, uplo, m_BY, n_BY, + alpha, A_sparse, B.data(), ldb, + beta, Y_actual.data(), ldy); + } // Tolerance: the dense reference (blas::symm) and the sparse path // accumulate in different orders, so we get a few ULPs of accumulation @@ -207,42 +220,16 @@ TEST_F(TestSpsymm, CaseD_SparseSparseThrows) { }, RandBLAS::Error); } +// Routes through the public RandBLAS::spsymm(Symmetric) wrapper +// overload instead of the lower-level RandBLAS::sparse_data::spsymm. +// All other setup (dense reference, comparison tolerance) is identical +// to run_case; we set route_via_wrapper=true to flip the dispatch. TEST_F(TestSpsymm, SymmetricWrapper) { - using SpMat = CSRMatrix; - int64_t n_A = 8; - int64_t d = 3; - int64_t lda = n_A; - std::vector A_full(lda * n_A, 0.0); - fill_sym_dense(n_A, A_full.data(), lda, 0); - std::vector A_tri(A_full); - zero_other_triangle(n_A, A_tri.data(), lda, Uplo::Upper); - SpMat A_sparse(n_A, n_A); - dense_to_csr(Layout::ColMajor, A_tri.data(), 0.0, A_sparse); - - int64_t m_BY = n_A, n_BY = d; - int64_t ldb = m_BY, ldy = m_BY; - std::vector B(m_BY * n_BY); - DenseDist DB(m_BY, n_BY, ScalarDist::Uniform); - RandBLAS::fill_dense_unpacked(Layout::ColMajor, DB, m_BY, n_BY, 0, 0, B.data(), RNGState(11)); - std::vector Y_actual(m_BY * n_BY, 0.0); - std::vector Y_expect = Y_actual; - - double alpha = 1.0, beta = 0.0; - blas::symm(Layout::ColMajor, Side::Left, Uplo::Upper, m_BY, n_BY, - alpha, A_full.data(), lda, B.data(), ldb, - beta, Y_expect.data(), ldy); - - // Route via wrapper overload at the public RandBLAS:: namespace - auto A_sym = RandBLAS::as_symmetric(A_sparse, Uplo::Upper); - RandBLAS::spsymm(Layout::ColMajor, m_BY, n_BY, alpha, A_sym, - B.data(), ldb, beta, Y_actual.data(), ldy); - - double atol = 100.0 * std::numeric_limits::epsilon(); - double rtol = 10.0 * std::numeric_limits::epsilon(); - auto msg = RandBLAS::testing::matrices_approx_equal( - Layout::ColMajor, blas::Op::NoTrans, m_BY, n_BY, - Y_actual.data(), ldy, Y_expect.data(), ldy, - __PRETTY_FUNCTION__, __FILE__, __LINE__, atol, rtol + run_case>( + Layout::ColMajor, Side::Left, Uplo::Upper, + /*n_A=*/8, /*d=*/3, + /*alpha=*/1.0, /*beta=*/0.0, + /*seed_A=*/0, /*seed_B=*/11, + /*route_via_wrapper=*/true ); - if (!msg.empty()) FAIL() << msg; } From c0517eb7a9e3ecfef94ab43d259948f61c0a2961 Mon Sep 17 00:00:00 2001 From: mmelnich Date: Thu, 14 May 2026 12:27:26 -0400 Subject: [PATCH 11/19] Extract mkl_sparse_mm_call: shared type-dispatched wrapper around MKL'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 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). --- RandBLAS/sparse_data/mkl_spmm_impl.hh | 53 +++++++++++++++++-------- RandBLAS/sparse_data/mkl_spsymm_impl.hh | 23 +++-------- 2 files changed, 42 insertions(+), 34 deletions(-) diff --git a/RandBLAS/sparse_data/mkl_spmm_impl.hh b/RandBLAS/sparse_data/mkl_spmm_impl.hh index f027b633..1a3e0c27 100644 --- a/RandBLAS/sparse_data/mkl_spmm_impl.hh +++ b/RandBLAS/sparse_data/mkl_spmm_impl.hh @@ -257,6 +257,38 @@ MKLSparseHandle make_mkl_handle(const SpMat& A) { } } +// ============================================================================ +// Type-dispatched wrapper around mkl_sparse_d_mm / mkl_sparse_s_mm. +// Shared between mkl_left_spmm (general A) and mkl_spsymm (symmetric A); +// the caller controls the matrix_descr (general vs symmetric) and the +// `op` flag, plus the post-call status interpretation. +// ============================================================================ +template +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 +) { + if constexpr (std::is_same_v) { + return mkl_sparse_d_mm( + op, alpha, A_handle, descr, mkl_layout, + B, (MKL_INT)n_rhs, (MKL_INT)ldb, + beta, C, (MKL_INT)ldc + ); + } else if constexpr (std::is_same_v) { + return mkl_sparse_s_mm( + op, alpha, A_handle, descr, mkl_layout, + B, (MKL_INT)n_rhs, (MKL_INT)ldb, + beta, C, (MKL_INT)ldc + ); + } else { + static_assert(sizeof(T) == 0, "MKL sparse BLAS only supports float and double."); + // see GitHub PR #155 for why we don't use static_assert(false, ...). + } +} + // ============================================================================ // MKL-accelerated left_spmm: C = alpha * op(A) * op(B) + beta * C // where A is sparse, B and C are dense. @@ -313,22 +345,11 @@ bool mkl_left_spmm( struct matrix_descr descr; descr.type = SPARSE_MATRIX_TYPE_GENERAL; - sparse_status_t status; - if constexpr (std::is_same_v) { - status = mkl_sparse_d_mm( - to_mkl_op(opA), alpha, h.handle, descr, - to_mkl_layout(layout), - B, (MKL_INT)n, (MKL_INT)ldb, - beta, C, (MKL_INT)ldc - ); - } else if constexpr (std::is_same_v) { - status = mkl_sparse_s_mm( - to_mkl_op(opA), alpha, h.handle, descr, - to_mkl_layout(layout), - B, (MKL_INT)n, (MKL_INT)ldb, - beta, C, (MKL_INT)ldc - ); - } + sparse_status_t status = mkl_sparse_mm_call( + to_mkl_op(opA), alpha, h.handle, descr, + to_mkl_layout(layout), + B, n, ldb, beta, C, ldc + ); check_mkl_status(status, "mkl_sparse_mm"); return true; // signal: MKL handled it } diff --git a/RandBLAS/sparse_data/mkl_spsymm_impl.hh b/RandBLAS/sparse_data/mkl_spsymm_impl.hh index 43bc504f..614af92a 100644 --- a/RandBLAS/sparse_data/mkl_spsymm_impl.hh +++ b/RandBLAS/sparse_data/mkl_spsymm_impl.hh @@ -115,24 +115,11 @@ bool mkl_spsymm( : SPARSE_FILL_MODE_LOWER; descr.diag = SPARSE_DIAG_NON_UNIT; - sparse_status_t status; - if constexpr (std::is_same_v) { - status = mkl_sparse_d_mm( - SPARSE_OPERATION_NON_TRANSPOSE, alpha, h.handle, descr, - to_mkl_layout(layout), - B, (MKL_INT)n, (MKL_INT)ldb, - beta, Y, (MKL_INT)ldy - ); - } else if constexpr (std::is_same_v) { - status = mkl_sparse_s_mm( - SPARSE_OPERATION_NON_TRANSPOSE, alpha, h.handle, descr, - to_mkl_layout(layout), - B, (MKL_INT)n, (MKL_INT)ldb, - beta, Y, (MKL_INT)ldy - ); - } else { - static_assert(sizeof(T) == 0, "MKL sparse BLAS only supports float and double."); - } + sparse_status_t status = mkl_sparse_mm_call( + SPARSE_OPERATION_NON_TRANSPOSE, alpha, h.handle, descr, + to_mkl_layout(layout), + B, n, ldb, beta, Y, ldy + ); // Some MKL versions return NOT_SUPPORTED for combinations we couldn't // predict. Don't throw -- signal fallback. From 992a2f3808c7116f5c46ba9ec23eaac87b2eca0d Mon Sep 17 00:00:00 2001 From: mmelnich Date: Thu, 14 May 2026 12:35:30 -0400 Subject: [PATCH 12/19] Sweep pre-existing safe_scal-in-loop matrix-scale sites to use lascl 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. --- RandBLAS/sparse_data/spmm_dispatch.hh | 6 ++---- RandBLAS/sparse_data/trsm_dispatch.hh | 6 ++---- examples/simple-kernel-benchmarks/spmm_performance.cc | 8 +------- examples/sparse-low-rank-approx/qrcp_matrixmarket.cc | 5 +++-- 4 files changed, 8 insertions(+), 17 deletions(-) diff --git a/RandBLAS/sparse_data/spmm_dispatch.hh b/RandBLAS/sparse_data/spmm_dispatch.hh index c2edf112..192b6aee 100644 --- a/RandBLAS/sparse_data/spmm_dispatch.hh +++ b/RandBLAS/sparse_data/spmm_dispatch.hh @@ -31,6 +31,7 @@ #include "RandBLAS/base.hh" #include "RandBLAS/exceptions.hh" +#include "RandBLAS/util.hh" #include "RandBLAS/sparse_data/base.hh" #include "RandBLAS/sparse_data/coo_matrix.hh" #include "RandBLAS/sparse_data/csr_matrix.hh" @@ -125,14 +126,11 @@ void left_spmm( if (layout == Layout::ColMajor) { randblas_require(ldb >= rows_B); randblas_require(ldc >= d); - for (int64_t i = 0; i < n; ++i) - RandBLAS::util::safe_scal(d, beta, &C[i*ldc]); } else { randblas_require(ldc >= n); randblas_require(ldb >= cols_B); - for (int64_t i = 0; i < d; ++i) - RandBLAS::util::safe_scal(n, beta, &C[i*ldc]); } + RandBLAS::util::lascl(layout, d, n, beta, C, ldc); if (alpha == (T) 0) return; diff --git a/RandBLAS/sparse_data/trsm_dispatch.hh b/RandBLAS/sparse_data/trsm_dispatch.hh index 39ad2b9a..d2dcbb0e 100644 --- a/RandBLAS/sparse_data/trsm_dispatch.hh +++ b/RandBLAS/sparse_data/trsm_dispatch.hh @@ -31,6 +31,7 @@ #include "RandBLAS/base.hh" #include "RandBLAS/exceptions.hh" +#include "RandBLAS/util.hh" #include "RandBLAS/sparse_data/base.hh" #include "RandBLAS/sparse_data/coo_matrix.hh" #include "RandBLAS/sparse_data/csr_matrix.hh" @@ -199,13 +200,10 @@ void trsm( int64_t m = A.n_rows; if (layout == blas::Layout::ColMajor) { randblas_require(ldb >= m); - for (int64_t i = 0; i < n; ++i) - RandBLAS::util::safe_scal(m, alpha, &B[i*ldb]); } else { randblas_require(ldb >= n); - for (int64_t i = 0; i < m; ++i) - RandBLAS::util::safe_scal(n, alpha, &B[i*ldb]); } + RandBLAS::util::lascl(layout, m, n, alpha, B, ldb); if (alpha == static_cast(0)) return; diff --git a/examples/simple-kernel-benchmarks/spmm_performance.cc b/examples/simple-kernel-benchmarks/spmm_performance.cc index 03bb7c43..411fb055 100644 --- a/examples/simple-kernel-benchmarks/spmm_performance.cc +++ b/examples/simple-kernel-benchmarks/spmm_performance.cc @@ -88,13 +88,7 @@ void handrolled_left_spmm_csr( const T *B, int64_t ldb, T beta, T *C, int64_t ldc ) { // Apply beta to C (same as dispatch) - if (layout == Layout::ColMajor) { - for (int64_t i = 0; i < n; ++i) - RandBLAS::util::safe_scal(d, beta, &C[i*ldc]); - } else { - for (int64_t i = 0; i < d; ++i) - RandBLAS::util::safe_scal(n, beta, &C[i*ldc]); - } + RandBLAS::util::lascl(layout, d, n, beta, C, ldc); if (alpha == (T)0) return; // Call the hand-rolled CSR kernel directly (same as dispatch fallback) diff --git a/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc b/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc index 7c4f23af..49a56cbf 100644 --- a/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc +++ b/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc @@ -341,9 +341,10 @@ void sketch_to_tqrcp(SpMat &A, int64_t k, T* Q, int64_t ldq, T* Y, int64_t ldy, lapack::geqp3(k, n, Y, ldy, piv, tau), "GEQP3 : ") // ================================================================ - // Step 2: copy A(:, piv(0)-1), ..., A(:, piv(k)-1) into dense Q + // Step 2: copy A(:, piv(0)-1), ..., A(:, piv(k)-1) into dense Q. + // Zero the m-by-k destination block in one pass before scattering. + RandBLAS::util::lascl(blas::Layout::ColMajor, m, k, 0.0, Q, ldq); for (int64_t j = 0; j < k; ++j) { - RandBLAS::util::safe_scal(m, 0.0, Q + j*ldq); for (int64_t ell = A.colptr[piv[j]-1]; ell < A.colptr[piv[j]]; ++ell) { int64_t i = A.rowidxs[ell]; Q[i + ldq*j] = A.vals[ell]; From e7b9a517aba5b303ede5a31e37d5cfd35c2a6f5c Mon Sep 17 00:00:00 2001 From: mmelnich Date: Thu, 14 May 2026 12:54:03 -0400 Subject: [PATCH 13/19] mkl_spsymm: handle side=Right via opposite-layout reinterpretation 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. --- RandBLAS/sparse_data/mkl_spsymm_impl.hh | 36 +++++++++++++++++-------- 1 file changed, 25 insertions(+), 11 deletions(-) diff --git a/RandBLAS/sparse_data/mkl_spsymm_impl.hh b/RandBLAS/sparse_data/mkl_spsymm_impl.hh index 614af92a..9848f73c 100644 --- a/RandBLAS/sparse_data/mkl_spsymm_impl.hh +++ b/RandBLAS/sparse_data/mkl_spsymm_impl.hh @@ -57,13 +57,20 @@ namespace RandBLAS::sparse_data::mkl { // hand-rolled per-format kernel. MKL applies alpha and beta internally; the // caller therefore passes them through unchanged. // -// Known limitations that trigger a fallback: -// - side == Right: MKL's mkl_sparse_d_mm has no side parameter; the -// symmetric matrix is always on the left of the dense block. The -// transpose-trick to express side=Right via layout flips depends on -// leading-dim assumptions that don't always hold; safer to fall back. +// Known limitation that triggers a fallback: // - Index type mismatched with MKL_INT. // +// Both Side values are handled: +// - side=Left (Y = alpha*A*B + beta*Y): direct call to mkl_sparse_d_mm +// with the symmetric descriptor. +// - side=Right (Y = alpha*B*A + beta*Y): A is symmetric so A == A^T, +// and (B*A)^T = A^T * B^T = A * B^T. We tell MKL to compute A*B^T +// into Y^T by flipping the MKL layout flag. Concretely, the same +// user buffers for B and Y are reinterpreted in the opposite layout +// (ColMajor <-> RowMajor); ldb and ldy carry through unchanged +// because the reinterpretation has the same leading-dim semantics +// (rows-of-RowMajor have the same stride as cols-of-ColMajor). +// // CSC handling: MKL's mkl_sparse_d_mm returns NOT_SUPPORTED for CSC even // with a symmetric descriptor. We work around this by taking the // CSC.transpose() view (a lightweight CSR view over the same buffers, @@ -90,9 +97,6 @@ bool mkl_spsymm( using sint_t = typename SpMat::index_t; constexpr bool is_csc = std::is_same_v>; - if (side != blas::Side::Left) - return false; - if constexpr (is_csc) { // Symmetric A: A == A^T. The CSC->CSR view is lightweight (same // buffers, reinterpreted) and gives us a CSR matrix MKL accepts; @@ -115,10 +119,21 @@ bool mkl_spsymm( : SPARSE_FILL_MODE_LOWER; descr.diag = SPARSE_DIAG_NON_UNIT; + // For side=Right, reinterpret B and Y in the opposite layout to + // present them to MKL as B^T and Y^T; MKL then computes Y^T = A*B^T + // which is the transpose of Y = B*A. The number of MKL-side + // right-hand-side columns is m (the user's row count) rather than n + // (the user's col count, = A's order). + blas::Layout mkl_layout = (side == blas::Side::Left) + ? layout + : (layout == blas::Layout::ColMajor ? blas::Layout::RowMajor + : blas::Layout::ColMajor); + int64_t n_rhs = (side == blas::Side::Left) ? n : m; + sparse_status_t status = mkl_sparse_mm_call( SPARSE_OPERATION_NON_TRANSPOSE, alpha, h.handle, descr, - to_mkl_layout(layout), - B, n, ldb, beta, Y, ldy + to_mkl_layout(mkl_layout), + B, n_rhs, ldb, beta, Y, ldy ); // Some MKL versions return NOT_SUPPORTED for combinations we couldn't @@ -127,7 +142,6 @@ bool mkl_spsymm( return false; check_mkl_status(status, "mkl_sparse_mm (symmetric)"); - (void) m; // m is implied by A.n_rows; kept for signature symmetry return true; } From a7005f18d61da98b7178e0b523c03e3c786fe904 Mon Sep 17 00:00:00 2001 From: mmelnich Date: Thu, 14 May 2026 12:58:21 -0400 Subject: [PATCH 14/19] sksy: transpose-copy S into matching layout for SYMM on layout-mismatch 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 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. --- RandBLAS/sksy.hh | 37 +++++++++++++++++++++++++------------ 1 file changed, 25 insertions(+), 12 deletions(-) diff --git a/RandBLAS/sksy.hh b/RandBLAS/sksy.hh index c26914a2..c852e185 100644 --- a/RandBLAS/sksy.hh +++ b/RandBLAS/sksy.hh @@ -55,10 +55,10 @@ namespace RandBLAS::dense { /// `submatrix_as_blackbox` (same pattern as `lskge3`). When the buffered S's /// storage layout matches the caller's `layout`, the final call is /// `blas::symm` with `side = Right` (since A is on the right of S in the -/// operation). When layouts mismatch, SYMM cannot transpose S on-the-fly; -/// this first cut falls back to `blas::gemm` with `opS = Trans`, sacrificing -/// the SYMM speedup on that path. Optimization target: transpose-copy of S -/// into matching layout, then SYMM. See project-plans/randblas-symm-plan.md §7. +/// operation). When layouts mismatch, SYMM cannot transpose S on the fly, so +/// we transpose-copy S into a tight buffer matching the caller's layout (cost: +/// `O(d * n)` for the copy) and then call SYMM on the copy. This keeps the +/// SYMM speedup on the matvec, at the cost of the one-time copy. template void lsksy3( blas::Layout layout, @@ -102,11 +102,17 @@ void lsksy3( blas::symm(layout, blas::Side::Right, uplo, d, n, alpha, A, lda, S_ptr, lds, beta, B, ldb); } else { - // Layout mismatch: SYMM has no opS flag. Fall back to GEMM with - // opS = Trans. Loses the SYMM 1.3-1.8x speedup on this path but is - // correct. Optimization target tracked in randblas-symm-plan.md §7. - blas::gemm(layout, blas::Op::Trans, blas::Op::NoTrans, d, n, n, - alpha, S_ptr, lds, A, lda, beta, B, ldb); + // Layout mismatch: transpose-copy S into a tight buffer in the + // caller's layout, then SYMM. Costs O(d*n) for the copy, but keeps + // the SYMM speedup (~1.3-1.8x over GEMM) on the d*n*n_A matvec. + int64_t lds_new = (layout == blas::Layout::ColMajor) ? d : n; + std::vector S_copy(static_cast(d) * static_cast(n)); + auto [irs_in, ics_in] = layout_to_strides(S.layout, lds); + auto [irs_out, ics_out] = layout_to_strides(layout, lds_new); + util::omatcopy(d, n, S_ptr, irs_in, ics_in, + S_copy.data(), irs_out, ics_out); + blas::symm(layout, blas::Side::Right, uplo, d, n, + alpha, A, lda, S_copy.data(), lds_new, beta, B, ldb); } return; } @@ -164,9 +170,16 @@ void rsksy3( blas::symm(layout, blas::Side::Left, uplo, n, d, alpha, A, lda, S_ptr, lds, beta, B, ldb); } else { - // Layout-mismatch fallback (see lsksy3). - blas::gemm(layout, blas::Op::NoTrans, blas::Op::Trans, n, d, n, - alpha, A, lda, S_ptr, lds, beta, B, ldb); + // Layout mismatch: transpose-copy S into a tight matching-layout + // buffer, then SYMM. See lsksy3 for the trade-off discussion. + int64_t lds_new = (layout == blas::Layout::ColMajor) ? n : d; + std::vector S_copy(static_cast(n) * static_cast(d)); + auto [irs_in, ics_in] = layout_to_strides(S.layout, lds); + auto [irs_out, ics_out] = layout_to_strides(layout, lds_new); + util::omatcopy(n, d, S_ptr, irs_in, ics_in, + S_copy.data(), irs_out, ics_out); + blas::symm(layout, blas::Side::Left, uplo, n, d, + alpha, A, lda, S_copy.data(), lds_new, beta, B, ldb); } return; } From 096b06352d245f8571f8f4ed7e7925d58a1e0adb Mon Sep 17 00:00:00 2001 From: mmelnich Date: Thu, 14 May 2026 15:42:18 -0400 Subject: [PATCH 15/19] Add spsymm / sketch_symmetric perf benchmark; expose spsymm via RandBLAS.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 to the umbrella header so downstream code using #include gets RandBLAS::spsymm and the Symmetric 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). --- RandBLAS.hh | 1 + examples/CMakeLists.txt | 10 + .../spsymm_performance.cc | 296 ++++++++++++++++++ 3 files changed, 307 insertions(+) create mode 100644 examples/simple-kernel-benchmarks/spsymm_performance.cc diff --git a/RandBLAS.hh b/RandBLAS.hh index d789a753..4cd5f7b0 100644 --- a/RandBLAS.hh +++ b/RandBLAS.hh @@ -39,5 +39,6 @@ #include #include #include +#include #endif diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index a953bb35..4a3b4938 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -103,5 +103,15 @@ target_link_libraries( spmm_performance PUBLIC RandBLAS blaspp lapackpp ) +add_executable( + spsymm_performance simple-kernel-benchmarks/spsymm_performance.cc +) +target_include_directories( + spsymm_performance PUBLIC ${Random123_DIR} +) +target_link_libraries( + spsymm_performance PUBLIC RandBLAS blaspp lapackpp +) + diff --git a/examples/simple-kernel-benchmarks/spsymm_performance.cc b/examples/simple-kernel-benchmarks/spsymm_performance.cc new file mode 100644 index 00000000..bdebf44a --- /dev/null +++ b/examples/simple-kernel-benchmarks/spsymm_performance.cc @@ -0,0 +1,296 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// ============================================================================ +// SPSYMM / SKETCH_SYMMETRIC PERFORMANCE BENCHMARK +// ============================================================================ +// +// This benchmark answers two performance questions about the symm-kernels +// added in https://github.com/BallisticLA/RandBLAS/pull/163 : +// +// 1. Sparse: How much does the new RandBLAS::spsymm (one-triangle storage, +// MKL fast path via SPARSE_MATRIX_TYPE_SYMMETRIC) speed up over the +// "both-triangles + right_spmm" workaround that downstream code had to +// use before this PR? +// +// 2. Dense: How much does the rewritten RandBLAS::sketch_symmetric (now +// backed by blas::symm) speed up over its pre-PR behaviour, which +// silently forwarded to sketch_general -> lskge3 -> blas::gemm with no +// symmetry exploitation? +// +// NOTATION: +// A_symm - symmetric matrix of order n_A (dense or sparse, one triangle) +// B - dense matrix (n_A x d for side=Left; d x n_A for side=Right) +// Y - dense result matrix (same shape as B) +// density - fraction of upper-triangle entries of A_symm that are nonzero +// (the implied lower-triangle entries follow by symmetry) +// +// USAGE: +// ./spsymm_performance # default sweep +// ./spsymm_performance n_A d density [num_trials] # single config +// +// Defaults: num_trials=10, density=0.05. +// +// EXAMPLES: +// env OMP_NUM_THREADS=8 ./spsymm_performance +// env OMP_NUM_THREADS=8 ./spsymm_performance 2000 200 0.01 +// +// ============================================================================ + +#include + +#include +#include +#include +#include +#include +#include +#include + +using std::chrono::steady_clock; +using std::chrono::duration_cast; +using std::chrono::microseconds; +using blas::Layout; +using blas::Uplo; +using blas::Side; +using blas::Op; + +using T = double; +using sint_t = int64_t; +using SpMat = RandBLAS::sparse_data::CSRMatrix; + + +// ============================================================================ +// Helpers: random dense symmetric matrix + sparse counterparts. +// ============================================================================ + +void make_dense_symmetric(int64_t n, std::vector& A, uint64_t seed) { + A.assign(static_cast(n) * n, T(0)); + std::mt19937_64 rng(seed); + std::uniform_real_distribution uni(-1.0, 1.0); + for (int64_t j = 0; j < n; ++j) { + for (int64_t i = 0; i <= j; ++i) { + T v = uni(rng); + A[i + j * n] = v; + if (i != j) A[j + i * n] = v; // mirror + } + } +} + +// Sparsify the upper triangle of A (including diagonal) at the given density, +// then mirror into the lower triangle. +void sparsify_upper_then_mirror(int64_t n, std::vector& A, double density, uint64_t seed) { + std::mt19937_64 rng(seed); + std::uniform_real_distribution uni01(0.0, 1.0); + for (int64_t j = 0; j < n; ++j) { + for (int64_t i = 0; i <= j; ++i) { + if (uni01(rng) >= density) { + A[i + j * n] = T(0); + if (i != j) A[j + i * n] = T(0); + } + } + } +} + +// Build a one-triangle (Upper) CSR view by walking the upper triangle of A_dense. +SpMat build_csr_upper_only(int64_t n, const std::vector& A_dense) { + std::vector rowptr(n + 1, 0); + for (int64_t i = 0; i < n; ++i) { + for (int64_t j = i; j < n; ++j) + if (A_dense[i + j * n] != T(0)) ++rowptr[i + 1]; + } + for (int64_t i = 0; i < n; ++i) rowptr[i + 1] += rowptr[i]; + int64_t nnz = rowptr[n]; + + SpMat A_sparse(n, n); + A_sparse.reserve(nnz); + std::vector tmp_rp(rowptr); // running cursor per row + for (int64_t i = 0; i < n; ++i) { + for (int64_t j = i; j < n; ++j) { + T v = A_dense[i + j * n]; + if (v != T(0)) { + int64_t pos = tmp_rp[i]++; + A_sparse.colidxs[pos] = static_cast(j); + A_sparse.vals[pos] = v; + } + } + } + for (int64_t i = 0; i <= n; ++i) A_sparse.rowptr[i] = rowptr[i]; + return A_sparse; +} + +// Build a full (both-triangles-stored) CSR for the pre-PR workaround comparison. +SpMat build_csr_both_triangles(int64_t n, const std::vector& A_dense) { + std::vector rowptr(n + 1, 0); + for (int64_t i = 0; i < n; ++i) + for (int64_t j = 0; j < n; ++j) + if (A_dense[i + j * n] != T(0)) ++rowptr[i + 1]; + for (int64_t i = 0; i < n; ++i) rowptr[i + 1] += rowptr[i]; + int64_t nnz = rowptr[n]; + + SpMat A_sparse(n, n); + A_sparse.reserve(nnz); + std::vector tmp_rp(rowptr); + for (int64_t i = 0; i < n; ++i) { + for (int64_t j = 0; j < n; ++j) { + T v = A_dense[i + j * n]; + if (v != T(0)) { + int64_t pos = tmp_rp[i]++; + A_sparse.colidxs[pos] = static_cast(j); + A_sparse.vals[pos] = v; + } + } + } + for (int64_t i = 0; i <= n; ++i) A_sparse.rowptr[i] = rowptr[i]; + return A_sparse; +} + + +// ============================================================================ +// Timing helper: median + min over num_trials. +// ============================================================================ +template +std::pair run_trials(Func&& func, int num_trials) { + std::vector times; + times.reserve(num_trials); + for (int t = 0; t < num_trials; ++t) { + auto start = steady_clock::now(); + func(); + auto end = steady_clock::now(); + times.push_back(duration_cast(end - start).count()); + } + std::sort(times.begin(), times.end()); + return {times[0], times[num_trials / 2]}; +} + +void print_row(const std::string& name, long min_us, long med_us, long baseline) { + double ratio = (baseline > 0) ? (double)min_us / baseline : 1.0; + std::cout << " " << std::setw(38) << std::left << name + << std::setw(10) << std::right << min_us + << std::setw(10) << med_us + << std::setw(10) << std::fixed << std::setprecision(2) << ratio << "x\n"; +} + + +// ============================================================================ +// run_config: one (n_A, d, density) point. Compares the symm-aware path +// against the workaround / pre-PR baselines. +// ============================================================================ +void run_config(int64_t n_A, int64_t d, double density, int num_trials) { + uint64_t seed = 12345; + Layout layout = Layout::ColMajor; + T alpha = T(1.0), beta = T(0.0); + + std::cout << "--- A is " << n_A << "x" << n_A + << " symmetric, B is " << n_A << "x" << d + << ", density=" << std::setprecision(4) << density + << " (median + min over " << num_trials << " trials) ---\n"; + + std::vector A_full(static_cast(n_A) * n_A); + make_dense_symmetric(n_A, A_full, seed); + sparsify_upper_then_mirror(n_A, A_full, density, seed + 1); + + SpMat A_csr_upper = build_csr_upper_only(n_A, A_full); + SpMat A_csr_full = build_csr_both_triangles(n_A, A_full); + + int64_t actual_nnz_upper = A_csr_upper.rowptr[n_A]; + int64_t actual_nnz_full = A_csr_full.rowptr[n_A]; + std::cout << " upper-triangle nnz = " << actual_nnz_upper + << " (full = " << actual_nnz_full << ")\n"; + + std::vector B(static_cast(n_A) * d); + { + std::mt19937_64 rng(seed + 2); + std::uniform_real_distribution uni(-1.0, 1.0); + for (auto& x : B) x = uni(rng); + } + std::vector Y(static_cast(n_A) * d, T(0)); + + std::cout << "\n SPARSE (side=Left, Y = A*B):\n"; + std::cout << " " << std::setw(38) << std::left << "kernel" + << std::setw(10) << std::right << "min(us)" + << std::setw(10) << "med(us)" + << std::setw(10) << "ratio\n"; + + // 1. RandBLAS::spsymm on one-triangle CSR (PR #163's new path). + auto [t1_min, t1_med] = run_trials([&] { + RandBLAS::spsymm(layout, Uplo::Upper, n_A, d, + alpha, A_csr_upper, B.data(), n_A, + beta, Y.data(), n_A); + }, num_trials); + + // 2. RandBLAS::spmm with full storage (the pre-PR workaround: + // both triangles materialized into a general-sparse CSR). + auto [t2_min, t2_med] = run_trials([&] { + RandBLAS::spmm(layout, Op::NoTrans, Op::NoTrans, + n_A, d, n_A, + alpha, A_csr_full, B.data(), n_A, + beta, Y.data(), n_A); + }, num_trials); + + // 3. Dense blas::symm on the fully populated dense A (the "ideal" + // BLAS-symm baseline -- doesn't exploit sparsity, but is the + // fastest possible dense SYMM). + auto [t3_min, t3_med] = run_trials([&] { + blas::symm(layout, Side::Left, Uplo::Upper, n_A, d, + alpha, A_full.data(), n_A, B.data(), n_A, + beta, Y.data(), n_A); + }, num_trials); + + print_row("RandBLAS::spsymm (one tri + MKL)", t1_min, t1_med, t1_min); + print_row("RandBLAS::spmm (full + MKL)", t2_min, t2_med, t1_min); + print_row("blas::symm (dense ref)", t3_min, t3_med, t1_min); + + // Dense-symmetric sketching comparison: new SYMM-backed vs the + // pre-PR GEMM-forwarding behaviour. + std::cout << "\n DENSE (sketch_symmetric vs sketch_general on dense-symm A):\n"; + std::cout << " " << std::setw(38) << std::left << "kernel" + << std::setw(10) << std::right << "min(us)" + << std::setw(10) << "med(us)" + << std::setw(10) << "ratio\n"; + + RandBLAS::DenseDist DS(n_A, d, RandBLAS::ScalarDist::Uniform); + RandBLAS::DenseSkOp S(DS, static_cast(seed + 3)); + RandBLAS::fill_dense(S); + + // 1. New: RandBLAS::sketch_symmetric (SYMM-backed via blas::symm). + auto [s1_min, s1_med] = run_trials([&] { + RandBLAS::sketch_symmetric(layout, Uplo::Upper, n_A, d, + alpha, A_full.data(), n_A, S, 0, 0, + beta, Y.data(), n_A); + }, num_trials); + + // 2. Old: equivalent call via sketch_general (the pre-PR behaviour + // when sketch_symmetric silently forwarded to GEMM). + auto [s2_min, s2_med] = run_trials([&] { + RandBLAS::sketch_general(layout, Op::NoTrans, Op::NoTrans, + n_A, d, n_A, + alpha, A_full.data(), n_A, S, 0, 0, + beta, Y.data(), n_A); + }, num_trials); + + print_row("sketch_symmetric (new SYMM path)", s1_min, s1_med, s1_min); + print_row("sketch_general (pre-PR GEMM path)", s2_min, s2_med, s1_min); + + std::cout << "\n"; +} + + +// ============================================================================ +// main: default 4-point sweep, or single config from argv. +// ============================================================================ +int main(int argc, char** argv) { + if (argc >= 4) { + int64_t n_A = std::stoll(argv[1]); + int64_t d = std::stoll(argv[2]); + double density = std::stod(argv[3]); + int num_trials = (argc >= 5) ? std::stoi(argv[4]) : 10; + run_config(n_A, d, density, num_trials); + } else { + int num_trials = 10; + std::cout << "Default sweep (n_A in {500, 1000, 2000}, d=200, density=0.05)\n\n"; + for (int64_t n_A : {500, 1000, 2000}) { + run_config(n_A, /*d=*/200, /*density=*/0.05, num_trials); + } + } + return 0; +} From 27cb43922d1490721dd3c4ce35882480a97da371 Mon Sep 17 00:00:00 2001 From: mmelnich Date: Thu, 14 May 2026 16:45:03 -0400 Subject: [PATCH 16/19] sksy: hand-roll Case B (dense-symm A x sparse SkOp); promote stub to 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 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 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. --- RandBLAS/sksy.hh | 193 ++++++++++++++++++--- RandBLAS/sparse_data/DevNotes.md | 51 ++++-- rtd/source/FAQ.rst | 11 +- rtd/source/api_reference/sketch_dense.rst | 8 +- rtd/source/api_reference/sketch_sparse.rst | 18 +- test/linops/test_sketch_symmetric.cc | 114 ++++++++++++ 6 files changed, 342 insertions(+), 53 deletions(-) diff --git a/RandBLAS/sksy.hh b/RandBLAS/sksy.hh index c852e185..59846919 100644 --- a/RandBLAS/sksy.hh +++ b/RandBLAS/sksy.hh @@ -38,7 +38,8 @@ // Symmetric sketching helpers (SYMM-backed). See project-plans/randblas-symm-plan.md // for the four-case API design and the implementation status of each case. // - lsksy3, rsksy3: Case A (dense-symm A x dense Omega), via blas::symm. -// - SparseSkOp branches of sketch_symmetric: Case B stubs (not implemented). +// - lsksys, rsksys: Case B (dense-symm A x sparse SkOp), per-stored-entry +// two-axpy scatter that reads only the uplo triangle of A. // ============================================================================= namespace RandBLAS::dense { @@ -187,33 +188,171 @@ void rsksy3( } // end namespace RandBLAS::dense -namespace RandBLAS { +namespace RandBLAS::sparse { -using namespace RandBLAS::dense; -using namespace RandBLAS::sparse; +// ============================================================================= +/// LSKSYS: dense symmetric A on the right of a SparseSkOp. +/// +/// Computes B = alpha * submat(S) * mat(A) + beta * B, where: +/// - mat(A) is n-by-n dense symmetric, with only the `uplo` triangle stored. +/// - submat(S) is the d-by-n view of S at (ro_s, co_s); S is a SparseSkOp. +/// - mat(B) is d-by-n dense. +/// +/// Inner loop: for each stored nonzero (row_S, col_S, v) of S inside the +/// submatrix window, contribute alpha*v * (row col_S of A) to row (row_S - ro_s) +/// of B. The row of symmetric A is assembled from the stored triangle as two +/// blas::axpy calls (one along the stored column, one along the stored row past +/// the diagonal). +template +void lsksys( + blas::Layout layout, + blas::Uplo uplo, + int64_t d, + int64_t n, + T alpha, + const SparseSkOp &S, + int64_t ro_s, + int64_t co_s, + const T *A, + int64_t lda, + T beta, + T *B, + int64_t ldb +) { + if (S.nnz < 0) { + SparseSkOp shallowcopy(S.dist, S.seed_state); + fill_sparse(shallowcopy); + lsksys(layout, uplo, d, n, alpha, shallowcopy, ro_s, co_s, A, lda, beta, B, ldb); + return; + } + util::lascl(layout, d, n, beta, B, ldb); + if (alpha == T(0)) return; + + auto Scoo = coo_view_of_skop(S); + const bool col_major = (layout == blas::Layout::ColMajor); + + for (int64_t p = 0; p < Scoo.nnz; ++p) { + sint_t row_S = Scoo.rows[p]; + sint_t col_S = Scoo.cols[p]; + if (row_S < ro_s || row_S >= ro_s + d) continue; + if (col_S < co_s || col_S >= co_s + n) continue; + int64_t i = static_cast(row_S) - ro_s; // B row + int64_t j = static_cast(col_S) - co_s; // A row index (sym A read as row j) + T av = alpha * Scoo.vals[p]; + + // B[i, :] += av * row_j(A_sym), split by uplo into two contiguous + // ranges of A so each becomes a single axpy. + if (uplo == blas::Uplo::Upper) { + // row j of A: + // c in [0, j]: read A[c, j] (column j of stored Upper, rows 0..j) + // c in (j, n-1]: read A[j, c] (row j of stored Upper, cols j+1..n-1) + if (col_major) { + blas::axpy(j + 1, av, &A[j * lda], 1, &B[i], ldb); + blas::axpy(n - j - 1, av, &A[j + (j + 1) * lda], lda, &B[i + (j + 1) * ldb], ldb); + } else { + blas::axpy(j + 1, av, &A[j], lda, &B[i * ldb], 1); + blas::axpy(n - j - 1, av, &A[j * lda + j + 1], 1, &B[i * ldb + j + 1], 1); + } + } else { + // Lower-stored: + // c in [0, j]: read A[j, c] (row j of stored Lower, cols 0..j) + // c in (j, n-1]: read A[c, j] (column j of stored Lower, rows j+1..n-1) + if (col_major) { + blas::axpy(j + 1, av, &A[j], lda, &B[i], ldb); + blas::axpy(n - j - 1, av, &A[(j + 1) + j * lda], 1, &B[i + (j + 1) * ldb], ldb); + } else { + blas::axpy(j + 1, av, &A[j * lda], 1, &B[i * ldb], 1); + blas::axpy(n - j - 1, av, &A[(j + 1) * lda + j], lda, &B[i * ldb + j + 1], 1); + } + } + } +} -namespace detail { // ============================================================================= -/// Shared Case-B stub-throw helper for the four SparseSkOp-taking -/// sketch_symmetric specializations. The variadic parameter pack consumes -/// the caller's arguments so the compiler doesn't warn about unused -/// parameters in the (genuinely-unused, by design) stub bodies. -template -[[noreturn]] inline void throw_sketch_symmetric_case_b(Args&&...) { - randblas_require( - false && - "RandBLAS::sketch_symmetric with a sparse sketching operator is the " - "Case-B kernel from project-plans/randblas-symm-plan.md. Not implemented " - "in this PR --- the API signature is reserved so future PRs can fill in " - "the body without breaking source compatibility. Composition fallback: " - "densify the SkOp then call sketch_symmetric on the dense version." - ); - __builtin_unreachable(); +/// RSKSYS: dense symmetric A on the left of a SparseSkOp. +/// +/// Computes B = alpha * mat(A) * submat(S) + beta * B, where: +/// - mat(A) is n-by-n dense symmetric, with only the `uplo` triangle stored. +/// - submat(S) is the n-by-d view of S at (ro_s, co_s); S is a SparseSkOp. +/// - mat(B) is n-by-d dense. +/// +/// Same two-axpy scatter pattern as lsksys, but on columns of A (the column +/// of A indexed by row_S of S contributes into the column of B indexed by +/// col_S - co_s). +template +void rsksys( + blas::Layout layout, + blas::Uplo uplo, + int64_t n, + int64_t d, + T alpha, + const T *A, + int64_t lda, + const SparseSkOp &S, + int64_t ro_s, + int64_t co_s, + T beta, + T *B, + int64_t ldb +) { + if (S.nnz < 0) { + SparseSkOp shallowcopy(S.dist, S.seed_state); + fill_sparse(shallowcopy); + rsksys(layout, uplo, n, d, alpha, A, lda, shallowcopy, ro_s, co_s, beta, B, ldb); + return; + } + + util::lascl(layout, n, d, beta, B, ldb); + if (alpha == T(0)) return; + + auto Scoo = coo_view_of_skop(S); + const bool col_major = (layout == blas::Layout::ColMajor); + + for (int64_t p = 0; p < Scoo.nnz; ++p) { + sint_t row_S = Scoo.rows[p]; + sint_t col_S = Scoo.cols[p]; + if (row_S < ro_s || row_S >= ro_s + n) continue; + if (col_S < co_s || col_S >= co_s + d) continue; + int64_t i = static_cast(row_S) - ro_s; // A column index + int64_t j = static_cast(col_S) - co_s; // B column + T av = alpha * Scoo.vals[p]; + + // B[:, j] += av * col_i(A_sym), split by uplo: + if (uplo == blas::Uplo::Upper) { + // col i of A: + // r in [0, i]: read A[r, i] (column i of stored Upper, rows 0..i) + // r in (i, n-1]: read A[i, r] (row i of stored Upper, cols i+1..n-1) + if (col_major) { + blas::axpy(i + 1, av, &A[i * lda], 1, &B[j * ldb], 1); + blas::axpy(n - i - 1, av, &A[i + (i + 1) * lda], lda, &B[(i + 1) + j * ldb], 1); + } else { + blas::axpy(i + 1, av, &A[i], lda, &B[j], ldb); + blas::axpy(n - i - 1, av, &A[i * lda + i + 1], 1, &B[(i + 1) * ldb + j], ldb); + } + } else { + // Lower-stored col i: + // r in [0, i-1]: read A[i, r] (row i of stored Lower, cols 0..i-1) + // r in [i, n-1]: read A[r, i] (column i of stored Lower, rows i..n-1) + if (col_major) { + blas::axpy(i, av, &A[i], lda, &B[j * ldb], 1); + blas::axpy(n - i, av, &A[i + i * lda], 1, &B[i + j * ldb], 1); + } else { + blas::axpy(i, av, &A[i * lda], 1, &B[j], ldb); + blas::axpy(n - i, av, &A[i * lda + i], lda, &B[i * ldb + j], ldb); + } + } + } } -} // namespace detail +} // end namespace RandBLAS::sparse + + +namespace RandBLAS { + +using namespace RandBLAS::dense; +using namespace RandBLAS::sparse; // MARK: SUBMAT(S) @@ -343,7 +482,7 @@ inline void sketch_symmetric( T beta, T* B, int64_t ldb ) { - detail::throw_sketch_symmetric_case_b(layout, uplo, d, n, alpha, S, ro_s, co_s, A, lda, beta, B, ldb); + RandBLAS::sparse::lsksys(layout, uplo, d, n, alpha, S, ro_s, co_s, A, lda, beta, B, ldb); } @@ -402,7 +541,7 @@ inline void sketch_symmetric( T beta, T* B, int64_t ldb ) { - detail::throw_sketch_symmetric_case_b(layout, uplo, n, d, alpha, A, lda, S, ro_s, co_s, beta, B, ldb); + RandBLAS::sparse::rsksys(layout, uplo, n, d, alpha, A, lda, S, ro_s, co_s, beta, B, ldb); } @@ -460,7 +599,9 @@ inline void sketch_symmetric( T beta, T* B, int64_t ldb ) { - detail::throw_sketch_symmetric_case_b(layout, uplo, alpha, S, A, lda, beta, B, ldb); + int64_t d = S.dist.n_rows; + int64_t n = S.dist.n_cols; + RandBLAS::sparse::lsksys(layout, uplo, d, n, alpha, S, 0, 0, A, lda, beta, B, ldb); } @@ -515,7 +656,9 @@ inline void sketch_symmetric( T beta, T* B, int64_t ldb ) { - detail::throw_sketch_symmetric_case_b(layout, uplo, alpha, S, A, lda, beta, B, ldb); + int64_t n = S.dist.n_rows; + int64_t d = S.dist.n_cols; + RandBLAS::sparse::rsksys(layout, uplo, n, d, alpha, A, lda, S, 0, 0, beta, B, ldb); } } // end namespace RandBLAS diff --git a/RandBLAS/sparse_data/DevNotes.md b/RandBLAS/sparse_data/DevNotes.md index 2c3b017f..7369921d 100644 --- a/RandBLAS/sparse_data/DevNotes.md +++ b/RandBLAS/sparse_data/DevNotes.md @@ -62,8 +62,8 @@ two operands (``A`` symmetric vs. the second factor ``B``): | Tag | Operation | A storage | B storage | Status this PR | |-----|---------------------------------|---------------------|-----------|-------------------------------------------------| -| A | dense-symm × dense | dense, one triangle | dense | Implemented via ``blas::symm`` in ``sksy.hh`` | -| B | dense-symm × sparse | dense, one triangle | sparse | Stub in ``sksy.hh`` (``sketch_symmetric`` on SparseSkOp). Throws ``RandBLAS::Error``. | +| 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 | Stub in ``spsymm_dispatch.hh`` (overload taking two ``SparseMatrix`` args). Throws ``RandBLAS::Error``. | @@ -72,7 +72,7 @@ two operands (``A`` symmetric vs. the second factor ``B``): | Tag | MKL native? | Notes | |-----|----------------|-----------------------------------------------------------------------------------------------------------------| | A | No (BLAS++) | ``blas::symm`` directly. | -| B | No | The transpose trick puts the sparse op on the left of ``mkl_sparse_d_mm``, but the dense A has no ``matrix_descr``, so MKL can't be told A is symmetric. | +| B | No | The transpose trick puts the sparse op on the left of ``mkl_sparse_d_mm``, but the dense A has no ``matrix_descr``, so MKL can't be told A is symmetric. We hand-roll instead --- see ``lsksys`` / ``rsksys`` in ``sksy.hh``. | | C | Yes (mostly) | ``mkl_sparse_d_mm`` with ``descr.type = SPARSE_MATRIX_TYPE_SYMMETRIC``. RandBLAS falls back to a hand kernel for side=Right (MKL has no Side parameter) and for CSC (``mkl_sparse_d_mm`` returns NOT_SUPPORTED on CSC). | | D | Partial | ``mkl_sparse_sp2m`` accepts a symmetric descriptor on A but writes sparse output; densifying afterward is a workable composition fallback but not a single-call kernel. | @@ -107,19 +107,38 @@ The public-facing wrappers in the top-level ``RandBLAS::`` namespace are: -- routes via the ``Symmetric`` carrier so the uplo annotation travels with the matrix. -### Cases B and D: why stub-only - -No portable kernel for "dense-symm × sparse" or "sparse-symm × sparse → dense" -exists in MKL, Ginkgo, the SparseBLAS reference implementation, the C++ -Sparse BLAS proposal (arXiv:2411.13259), or MAGMA-sparse (surveyed 2026-05). -Hand-rolling is feasible but non-trivial: case B's access pattern is "for -each nonzero of B, gather a column of A from the stored triangle", which is -awkward to vectorize; case D layers a symmetric-triangle filter on top of -an spgemm-into-dense pattern. - -The stubs exist so the API surface is locked: future PRs can fill in the -bodies without breaking source compatibility for downstream callers. The -stub message points back to ``project-plans/randblas-symm-plan.md``. +### Case B: hand-rolled in ``lsksys`` / ``rsksys`` + +Lives in ``sksy.hh`` next to the dense-SkOp helpers ``lsksy3`` / ``rsksy3``. +For each stored nonzero ``(i_S, j_S, v)`` of the SparseSkOp's COO view (with +submatrix offsets ``(ro_s, co_s)`` filtered inline), the kernel applies +``alpha * v`` to one row or column of the symmetric dense A and accumulates +into the corresponding row or column of the output. Reading "row j of A" (or +"col i of A") splits into two ranges based on the diagonal: + + - For ``Uplo::Upper``: the part above the diagonal walks A's stored row / + column directly; the part below comes from the symmetric reflection + (reading the transposed-position entry from the stored triangle). + - For ``Uplo::Lower``: the roles swap. + +Each range becomes a single ``blas::axpy`` with the appropriate stride, so +the inner-loop body is exactly two AXPY calls per stored nonzero of S +(plus a uniform layout / uplo branch). No special handling for the diagonal +beyond consistent inclusion in one of the two ranges. SparseSkOp is COO +internally, so submatrix filtering is a direct ``if (row < ro_s ...) continue`` +on the COO triples. + +### Case D: stub-only + +No portable kernel for "sparse-symm × sparse → dense" exists in MKL, +Ginkgo, the SparseBLAS reference implementation, the C++ Sparse BLAS +proposal (arXiv:2411.13259), or MAGMA-sparse (surveyed 2026-05). +Hand-rolling is feasible but the design space is messier than Case B +(mixed A/B sparse formats balloon the dispatch grid; the symmetric +filter sits on top of an spgemm-into-dense pattern). Deferred to a +future PR. Composition fallbacks documented in the stub's throw +message: densify B then call Case C, or use ``mkl_sparse_sp2m`` with +a symmetric descriptor on A and densify the resulting sparse C. ## Sketching sparse data with dense operators diff --git a/rtd/source/FAQ.rst b/rtd/source/FAQ.rst index 3172781a..169d25dc 100644 --- a/rtd/source/FAQ.rst +++ b/rtd/source/FAQ.rst @@ -110,10 +110,13 @@ No support for negative values of "incx" or "incy" in sketch_vector. to extend our SPMV kernels to support that, then we'd happily accept such a contribution. (It shouldn't be hard! We just haven't gotten around to this.) -SparseSkOp is not yet supported in sketch_symmetric. - ``sketch_symmetric`` now takes a ``blas::Uplo`` and dispatches to ``blas::symm`` for ``DenseSkOp``, exploiting symmetry of :math:`A`. The ``SparseSkOp`` branch throws ``RandBLAS::Error`` --- this is Case B in the SYMM-kernels plan - (``RandBLAS/sparse_data/DevNotes.md``). No portable kernel for "dense-symm matrix times sparse operator" exists in MKL, Ginkgo, the SparseBLAS reference, or MAGMA-sparse, so we deferred the hand-roll to a future PR. - As a composition fallback, ``DenseSkOp(densify(sparse_op))`` can be passed to ``sketch_symmetric`` instead --- it loses the operator-side sparsity benefit but exploits A's symmetry. +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 hand-rolled kernel (``lsksys`` / + ``rsksys``) that emits two ``blas::axpy`` calls per stored nonzero of the SkOp, + reading 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`` falls back to GEMM. When the ``DenseSkOp``'s storage layout differs from the caller's ``layout`` parameter, ``sketch_symmetric`` falls back to ``blas::gemm`` with the transpose flag --- ``blas::symm`` has no on-the-fly transpose flag for the dense operand. The layout-matched case still gets the SYMM speedup. diff --git a/rtd/source/api_reference/sketch_dense.rst b/rtd/source/api_reference/sketch_dense.rst index 88ae0611..1b784b23 100644 --- a/rtd/source/api_reference/sketch_dense.rst +++ b/rtd/source/api_reference/sketch_dense.rst @@ -69,9 +69,11 @@ Analogs to SYMM 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. Currently only ``DenseSkOp`` is supported; -calling these with a ``SparseSkOp`` throws ``RandBLAS::Error`` (Case B of the -SYMM-kernels plan; see ``RandBLAS/sparse_data/DevNotes.md``). +by symmetry and is **not** read. Both ``DenseSkOp`` and ``SparseSkOp`` are +supported: DenseSkOp dispatches to ``blas::symm`` via ``lsksy3`` / ``rsksy3``, +SparseSkOp dispatches to a hand-rolled two-axpy-per-nonzero kernel via +``lsksys`` / ``rsksys``. See ``RandBLAS/sparse_data/DevNotes.md`` for the +SparseSkOp access-pattern detail. .. dropdown:: :math:`\mtxB = \alpha \cdot \mtxS \cdot \mtxA + \beta \cdot \mtxB` :animate: fade-in-slide-down diff --git a/rtd/source/api_reference/sketch_sparse.rst b/rtd/source/api_reference/sketch_sparse.rst index 1b56f55a..42a3b760 100644 --- a/rtd/source/api_reference/sketch_sparse.rst +++ b/rtd/source/api_reference/sketch_sparse.rst @@ -137,11 +137,19 @@ Deterministic operations prevents a symmetric sparse matrix from being accidentally passed to the general ``spmm`` / ``spgemm`` routines. - Companion stubs exist for the cases where the second factor is also - sparse (Case D) or where the symmetric factor is dense and the second - is sparse (Case B, via ``sketch_symmetric``). They throw - ``RandBLAS::Error`` with a pointer to the design plan; no portable - kernel for those shapes exists in current sparse-BLAS libraries. + The "dense-symm A times sparse SkOp" case (Case B) is also implemented: + it lives in ``sketch_symmetric`` (the SparseSkOp branch) and uses a + hand-rolled two-axpy-per-nonzero 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 + a throw-stub: the new two-``SparseMatrix``-arg ``spsymm`` overload + throws ``RandBLAS::Error``. No portable reference kernel exists for + this shape; composition fallbacks are listed in the stub's throw + message (densify the sparse RHS and call Case C, or call + ``mkl_sparse_sp2m`` with a symmetric descriptor on A and densify + the resulting sparse output). .. dropdown:: :math:`\mtxB = \alpha \cdot \op(\mtxA)^{-1} \cdot \mtxB,` with sparse triangular :math:`\mtxA` :animate: fade-in-slide-down diff --git a/test/linops/test_sketch_symmetric.cc b/test/linops/test_sketch_symmetric.cc index 94383cb0..e7697384 100644 --- a/test/linops/test_sketch_symmetric.cc +++ b/test/linops/test_sketch_symmetric.cc @@ -31,6 +31,7 @@ #include "RandBLAS/base.hh" #include "RandBLAS/random_gen.hh" #include "RandBLAS/dense_skops.hh" +#include "RandBLAS/sparse_skops.hh" #include "RandBLAS/util.hh" #include "RandBLAS/sksy.hh" @@ -39,6 +40,8 @@ using blas::Uplo; using RandBLAS::ScalarDist; using RandBLAS::DenseDist; using RandBLAS::DenseSkOp; +using RandBLAS::SparseDist; +using RandBLAS::SparseSkOp; using RandBLAS::RNGState; using RandBLAS::Axis; @@ -169,6 +172,70 @@ class TestSketchSymmetric : public ::testing::Test { // only the triangle named by `uplo` is read --- the opposite triangle // contributing wrong values is undefined behavior at the caller, not a // RandBLAS check. The test was removed accordingly. + + // ============================================================================= + // Case B exerciser: sparse SkOp x dense symmetric A. Reference is the + // densified-sparse-skop fed to blas::symm. layout is explicit (SparseSkOp + // has no S.layout the way DenseSkOp does --- the COO storage is layout- + // agnostic, and the test layout determines how both B and the dense + // reference of S are laid out). + template + static void test_sparse_skop( + Layout layout, + uint32_t seed_a, uint32_t seed_skop, Axis major_axis, + T alpha, int64_t d, int64_t n, int64_t lda, T beta, + blas::Side side_skop, + int64_t vec_nnz = 2, + Uplo uplo = Uplo::Upper + ) { + auto [rows_out, cols_out] = dims_of_sketch_symmetric_output(d, n, side_skop); + std::vector A(lda * lda, T(0)); + random_symmetric_mat(n, A.data(), lda, RNGState(seed_a)); + + SparseDist DS(rows_out, cols_out, vec_nnz, major_axis); + SparseSkOp S(DS, seed_skop); + RandBLAS::fill_sparse(S); + + // Densify S into a buffer matching the requested layout, for the SYMM + // reference. lds is the major-axis leading dim. + int64_t lds = (layout == Layout::ColMajor) ? rows_out : cols_out; + int64_t ldb = lds; + std::vector S_dense(static_cast(rows_out) * cols_out, T(0)); + auto Scoo = RandBLAS::coo_view_of_skop(S); + for (int64_t p = 0; p < Scoo.nnz; ++p) { + int64_t r = Scoo.rows[p], c = Scoo.cols[p]; + if (layout == Layout::ColMajor) { + S_dense[r + c * lds] = Scoo.vals[p]; + } else { + S_dense[r * lds + c] = Scoo.vals[p]; + } + } + + uint32_t seed_b = seed_a + 42; + std::vector B_actual(static_cast(rows_out) * cols_out); + DenseDist DB(rows_out, cols_out, ScalarDist::Uniform); + RandBLAS::fill_dense_unpacked(layout, DB, rows_out, cols_out, 0, 0, B_actual.data(), RNGState(seed_b)); + std::vector B_expect = B_actual; + + auto side_a = sketch_symmetric_side( + side_skop, layout, uplo, rows_out, cols_out, + alpha, A.data(), lda, S, 0, 0, beta, B_actual.data(), ldb + ); + blas::symm(layout, side_a, uplo, rows_out, cols_out, + alpha, A.data(), lda, S_dense.data(), lds, beta, B_expect.data(), ldb); + + // Same FMA-order-divergence tolerance as the spsymm tests: the sparse + // path emits two AXPYs per stored nonzero (different accumulation + // order than dense SYMM on a fully-stored matrix). + T atol = T(100) * std::numeric_limits::epsilon(); + T rtol = T(10) * std::numeric_limits::epsilon(); + auto msg = RandBLAS::testing::matrices_approx_equal( + layout, blas::Op::NoTrans, rows_out, cols_out, + B_actual.data(), ldb, B_expect.data(), ldb, + __PRETTY_FUNCTION__, __FILE__, __LINE__, atol, rtol + ); + if (!msg.empty()) FAIL() << msg; + } }; @@ -370,3 +437,50 @@ TEST_F(TestSketchSymmetric, right_lift_opposing_layouts) { test_opposing_layouts(31, 33, Axis::Long, 0.5, 50, 10, 19, -1.0, blas::Side::Right); } + +// MARK: SPARSE SkOp (Case B). 4-axis sweep: +// side x layout x uplo x beta-mode (zero/nonzero). +// One seed pair per cell to keep the per-config trial count modest. + +TEST_F(TestSketchSymmetric, sparse_skop_left_colmajor_upper) { + test_sparse_skop(Layout::ColMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, 0.0, blas::Side::Left, 2, Uplo::Upper); + test_sparse_skop(Layout::ColMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, -1.0, blas::Side::Left, 2, Uplo::Upper); + test_sparse_skop(Layout::ColMajor, 31, 33, Axis::Long, 0.5, 3, 10, 19, 0.0, blas::Side::Left, 2, Uplo::Upper); + test_sparse_skop( Layout::ColMajor, 0, 1, Axis::Short, 0.5f, 3, 10, 10, 0.0f, blas::Side::Left, 2, Uplo::Upper); +} + +TEST_F(TestSketchSymmetric, sparse_skop_left_rowmajor_upper) { + test_sparse_skop(Layout::RowMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, 0.0, blas::Side::Left, 2, Uplo::Upper); + test_sparse_skop(Layout::RowMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, -1.0, blas::Side::Left, 2, Uplo::Upper); + test_sparse_skop(Layout::RowMajor, 31, 33, Axis::Long, 0.5, 3, 10, 19, 0.0, blas::Side::Left, 2, Uplo::Upper); +} + +TEST_F(TestSketchSymmetric, sparse_skop_right_colmajor_upper) { + test_sparse_skop(Layout::ColMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, 0.0, blas::Side::Right, 2, Uplo::Upper); + test_sparse_skop(Layout::ColMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, -1.0, blas::Side::Right, 2, Uplo::Upper); + test_sparse_skop(Layout::ColMajor, 31, 33, Axis::Long, 0.5, 3, 10, 19, 0.0, blas::Side::Right, 2, Uplo::Upper); +} + +TEST_F(TestSketchSymmetric, sparse_skop_right_rowmajor_upper) { + test_sparse_skop(Layout::RowMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, 0.0, blas::Side::Right, 2, Uplo::Upper); + test_sparse_skop(Layout::RowMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, -1.0, blas::Side::Right, 2, Uplo::Upper); + test_sparse_skop(Layout::RowMajor, 31, 33, Axis::Long, 0.5, 3, 10, 19, 0.0, blas::Side::Right, 2, Uplo::Upper); +} + +TEST_F(TestSketchSymmetric, sparse_skop_lower_triangle) { + // Uplo::Lower coverage --- one cell per (side, layout) combo, since the + // Upper-vs-Lower difference exercises the same kernel branches as + // ColMajor-vs-RowMajor (different inner-loop axpy strides). + test_sparse_skop(Layout::ColMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, 0.0, blas::Side::Left, 2, Uplo::Lower); + test_sparse_skop(Layout::RowMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, 0.0, blas::Side::Left, 2, Uplo::Lower); + test_sparse_skop(Layout::ColMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, 0.0, blas::Side::Right, 2, Uplo::Lower); + test_sparse_skop(Layout::RowMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, 0.0, blas::Side::Right, 2, Uplo::Lower); +} + +TEST_F(TestSketchSymmetric, sparse_skop_lift) { + // Embedding goes the wrong way (d > n): tests that the kernel handles + // non-square sketches in both directions. + test_sparse_skop(Layout::ColMajor, 0, 1, Axis::Short, 0.5, 13, 10, 10, 0.0, blas::Side::Left, 2, Uplo::Upper); + test_sparse_skop(Layout::ColMajor, 0, 1, Axis::Short, 0.5, 13, 10, 10, 0.0, blas::Side::Right, 2, Uplo::Upper); +} + From d4c42b720cdd83d010ccd54fc9067ff38a74581b Mon Sep 17 00:00:00 2001 From: mmelnich Date: Thu, 14 May 2026 18:50:16 -0400 Subject: [PATCH 17/19] spsymm: promote Case D from stub to densify-B + Case-C composition The two-SparseMatrix-arg spsymm overload (sparse-symm A x sparse B -> dense Y) now allocates an m-by-n std::vector 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. --- RandBLAS/sparse_data/DevNotes.md | 44 ++++--- RandBLAS/sparse_data/spsymm_dispatch.hh | 83 ++++++++------ rtd/source/api_reference/sketch_sparse.rst | 16 ++- test/linops/test_spsymm.cc | 127 ++++++++++++++++++--- 4 files changed, 200 insertions(+), 70 deletions(-) diff --git a/RandBLAS/sparse_data/DevNotes.md b/RandBLAS/sparse_data/DevNotes.md index 7369921d..b971038a 100644 --- a/RandBLAS/sparse_data/DevNotes.md +++ b/RandBLAS/sparse_data/DevNotes.md @@ -65,7 +65,7 @@ two operands (``A`` symmetric vs. the second factor ``B``): | 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 | Stub in ``spsymm_dispatch.hh`` (overload taking two ``SparseMatrix`` args). Throws ``RandBLAS::Error``. | +| D | sparse-symm × sparse → dense | sparse, one triangle | sparse | Implemented via densify-B + Case-C composition in ``spsymm_dispatch.hh``. | ### MKL availability @@ -74,7 +74,7 @@ two operands (``A`` symmetric vs. the second factor ``B``): | A | No (BLAS++) | ``blas::symm`` directly. | | B | No | The transpose trick puts the sparse op on the left of ``mkl_sparse_d_mm``, but the dense A has no ``matrix_descr``, so MKL can't be told A is symmetric. We hand-roll instead --- see ``lsksys`` / ``rsksys`` in ``sksy.hh``. | | C | Yes (mostly) | ``mkl_sparse_d_mm`` with ``descr.type = SPARSE_MATRIX_TYPE_SYMMETRIC``. RandBLAS falls back to a hand kernel for side=Right (MKL has no Side parameter) and for CSC (``mkl_sparse_d_mm`` returns NOT_SUPPORTED on CSC). | -| D | Partial | ``mkl_sparse_sp2m`` accepts a symmetric descriptor on A but writes sparse output; densifying afterward is a workable composition fallback but not a single-call kernel. | +| D | No | ``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. Symmetric expansion has to happen on the RandBLAS side, so we don't gain anything by routing through MKL. | ### Case C dispatch (``spsymm_dispatch.hh``) @@ -128,17 +128,35 @@ beyond consistent inclusion in one of the two ranges. SparseSkOp is COO internally, so submatrix filtering is a direct ``if (row < ro_s ...) continue`` on the COO triples. -### Case D: stub-only - -No portable kernel for "sparse-symm × sparse → dense" exists in MKL, -Ginkgo, the SparseBLAS reference implementation, the C++ Sparse BLAS -proposal (arXiv:2411.13259), or MAGMA-sparse (surveyed 2026-05). -Hand-rolling is feasible but the design space is messier than Case B -(mixed A/B sparse formats balloon the dispatch grid; the symmetric -filter sits on top of an spgemm-into-dense pattern). Deferred to a -future PR. Composition fallbacks documented in the stub's throw -message: densify B then call Case C, or use ``mkl_sparse_sp2m`` with -a symmetric descriptor on A and densify the resulting sparse C. +### Case D: densify-B + Case-C composition + +Lives in ``spsymm_dispatch.hh`` next to the Case-C dispatcher. The +overload taking two ``SparseMatrix`` operands allocates an ``m`` by +``n`` ``std::vector`` (tight leading dim 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``, +then calls the existing Case-C ``spsymm`` overload on the densified +buffer. Works in any build (the MKL fast path or the per-format hand +kernel both apply, depending on the build), and across all 3 × 3 = 9 +sparse-format pairings for ``(A, B)`` since the densification picks +the right format-specific helper. + +Why this composition rather than a single MKL ``sp2m`` call: + + - ``mkl_sparse_sp2m`` returns ``SPARSE_STATUS_NOT_SUPPORTED`` when + the ``matrix_descr`` on either operand is + ``SPARSE_MATRIX_TYPE_SYMMETRIC``; only ``GENERAL`` is accepted + there. + - ``mkl_sparse_d_spmmd`` (which writes directly to dense ``C``) + accepts 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 a +temporary dense buffer for ``B`` (cost ``O(m*n)``), and for the +typical RandNLA workload where ``B`` is a sketching operator with +``nnz(B) << m*n`` the buffer cost is small relative to the work that +would have to happen anyway. ``Y`` itself is never touched until the +Case-C call. ## Sketching sparse data with dense operators diff --git a/RandBLAS/sparse_data/spsymm_dispatch.hh b/RandBLAS/sparse_data/spsymm_dispatch.hh index 4c5012f2..20354f85 100644 --- a/RandBLAS/sparse_data/spsymm_dispatch.hh +++ b/RandBLAS/sparse_data/spsymm_dispatch.hh @@ -125,31 +125,30 @@ void spsymm( } // ============================================================================= -/// Case D stub: sparse-symmetric A times sparse B, dense output. +/// Case D: sparse-symmetric A times sparse B, dense output. /// /// @verbatim embed:rst:leading-slashes -/// Not implemented in this PR. The function signature is reserved here so that -/// future PRs can fill in the body without breaking source compatibility for -/// downstream users. Calling this overload triggers ``RandBLAS::Error`` via -/// ``randblas_require(false, ...)``. +/// Implemented as a composition: densify ``B`` into a temporary dense buffer +/// matching the caller's ``layout``, then call the Case-C ``spsymm`` overload +/// (sparse-symm × dense) on the result. This works in any build (the MKL +/// fast path or the per-format hand kernel both apply, depending on the +/// build), and across all 3 × 3 = 9 sparse-format pairings for ``(A, B)``, +/// since the densification picks the right format-specific helper. /// -/// Why no implementation? No portable kernel for "sparse-symmetric A times -/// sparse B into a dense result" exists in MKL, Ginkgo, the SparseBLAS -/// reference (``SparseBLAS/spblas-reference``), or MAGMA-sparse (all surveyed -/// 2026-05). The C++ Sparse BLAS standardization proposal (arXiv:2411.13259) -/// does not specify a SYMM-shaped sparse op either. Hand-rolling the kernel is -/// feasible but ~1-2 weeks of work, deferred outside the current PR. +/// Why this composition rather than a single MKL ``sp2m`` call: MKL's +/// ``mkl_sparse_sp2m`` returns ``SPARSE_STATUS_NOT_SUPPORTED`` when the +/// ``matrix_descr`` on either operand is ``SPARSE_MATRIX_TYPE_SYMMETRIC``; +/// only ``GENERAL`` is supported there. ``mkl_sparse_d_spmmd`` (which writes +/// directly to dense ``C``) accepts no descriptor at all. So the symmetric +/// expansion has to happen on the RandBLAS side. Composing through Case C +/// gets it for free at the cost of a temporary dense buffer for ``B`` +/// (cost: ``O(m * n)``), and for the typical RandNLA workload where ``B`` +/// is a sketching operator with ``nnz(B) << m*n`` the cost is small. /// -/// Composition fallbacks callers can use today: -/// -/// - Densify B and call the Case-C spsymm (above). Exploits A's symmetry -/// but loses B's sparsity benefit. -/// - With Intel MKL: ``mkl_sparse_sp2m`` accepts a symmetric ``matrix_descr`` -/// on A and produces a sparse C, which can be densified afterward. -/// Exploits both structures but pays an intermediate sparse-C plus a -/// dense-fill step. -/// -/// See project-plans/randblas-symm-plan.md for the full design context. +/// The dimensions of the densified ``B`` are the same as the user's ``Y``: +/// ``m``-by-``n``, with the user's ``layout`` and a tight leading dim. The +/// densified buffer is freed when the temporary ``std::vector`` goes out +/// of scope. ``Y`` itself is never touched until the Case-C call. /// @endverbatim template @@ -164,20 +163,34 @@ void spsymm( T beta, T* Y, int64_t ldy ) { - (void) layout; (void) side; (void) uplo; - (void) m; (void) n; (void) alpha; - (void) A; (void) B; (void) beta; - (void) Y; (void) ldy; - randblas_require( - false && - "RandBLAS::sparse_data::spsymm(..., const SpMatA& A, const SpMatB& B, ...) " - "(Case D: sparse-symmetric A times sparse B into dense Y) is not " - "implemented. No portable reference kernel exists in MKL, Ginkgo, " - "spblas-reference, or MAGMA-sparse (as of 2026-05). Composition " - "fallbacks: (a) densify B and call the Case-C spsymm, or (b) call " - "mkl_sparse_sp2m with a symmetric descriptor on A and densify the " - "resulting sparse C. See project-plans/randblas-symm-plan.md." - ); + static_assert(std::is_same_v, + "Case D: A and B must share scalar_t."); + + randblas_require(A.n_rows == A.n_cols); + int64_t k = (side == blas::Side::Left) ? m : n; + randblas_require(A.n_rows == k); + randblas_require(B.n_rows == m); + randblas_require(B.n_cols == n); + + // Densify B into a tight buffer in the caller's layout. + int64_t ldb_dense = (layout == blas::Layout::ColMajor) ? m : n; + std::vector B_dense(static_cast(m) * static_cast(n), T(0)); + + using sint_B = typename SpMatB::index_t; + if constexpr (std::is_same_v>) { + coo::coo_to_dense(B, layout, B_dense.data()); + } else if constexpr (std::is_same_v>) { + csr::csr_to_dense(B, layout, B_dense.data()); + } else if constexpr (std::is_same_v>) { + csc::csc_to_dense(B, layout, B_dense.data()); + } else { + static_assert(sizeof(SpMatB) == 0, + "RandBLAS::sparse_data::spsymm: SpMatB must be COO, CSR, or CSC."); + } + + // Compose into Case C with the densified B. + spsymm(layout, side, uplo, m, n, + alpha, A, B_dense.data(), ldb_dense, beta, Y, ldy); } } // end namespace RandBLAS::sparse_data diff --git a/rtd/source/api_reference/sketch_sparse.rst b/rtd/source/api_reference/sketch_sparse.rst index 42a3b760..8505c694 100644 --- a/rtd/source/api_reference/sketch_sparse.rst +++ b/rtd/source/api_reference/sketch_sparse.rst @@ -144,12 +144,16 @@ Deterministic operations access pattern. The "sparse-symm A times sparse RHS, dense output" case (Case D) is - a throw-stub: the new two-``SparseMatrix``-arg ``spsymm`` overload - throws ``RandBLAS::Error``. No portable reference kernel exists for - this shape; composition fallbacks are listed in the stub's throw - message (densify the sparse RHS and call Case C, or call - ``mkl_sparse_sp2m`` with a symmetric descriptor on A and densify - the resulting sparse output). + implemented via a two-``SparseMatrix``-arg ``spsymm`` overload that + densifies the sparse RHS into a temporary ``m``-by-``n`` buffer + (in the caller's layout) and then calls the Case-C dispatcher on + the densified buffer. This covers all 3 x 3 = 9 sparse-format + pairings for ``(A, B)``. ``mkl_sparse_sp2m`` rejects symmetric + descriptors and ``mkl_sparse_d_spmmd`` takes no descriptor, so + routing through MKL would not avoid the symmetric expansion --- + composing through Case C costs an ``O(m*n)`` temporary, which is + small for the typical RandNLA workload where ``B`` is a sketching + operator with ``nnz(B) << m*n``. .. dropdown:: :math:`\mtxB = \alpha \cdot \op(\mtxA)^{-1} \cdot \mtxB,` with sparse triangular :math:`\mtxA` :animate: fade-in-slide-down diff --git a/test/linops/test_spsymm.cc b/test/linops/test_spsymm.cc index 17f33b17..5792df3e 100644 --- a/test/linops/test_spsymm.cc +++ b/test/linops/test_spsymm.cc @@ -36,8 +36,10 @@ #include "RandBLAS/testing/comparison.hh" #include -#include +#include +#include #include +#include using namespace RandBLAS::sparse_data; using namespace RandBLAS::sparse_data::coo; @@ -175,6 +177,84 @@ class TestSpsymm : public ::testing::Test { } } } + + // Case D: sparse-symmetric A times sparse B -> dense Y. Reference is + // dense blas::symm on a fully-populated A and a densified B. + template + static void run_case_d( + Layout layout, Side side, Uplo uplo, + int64_t n_A, int64_t d, + T alpha, T beta, + uint32_t seed_A, uint32_t seed_B, + double density_B = 0.3 + ) { + int64_t m_BY, n_BY; + if (side == Side::Left) { m_BY = n_A; n_BY = d; } + else { m_BY = d; n_BY = n_A; } + + // Build dense symm A (full-storage reference). + int64_t lda = n_A; + std::vector A_full(lda * n_A, T(0)); + fill_sym_dense(n_A, A_full.data(), lda, seed_A); + // Build sparse A from a one-triangle copy. + std::vector A_tri(A_full); + zero_other_triangle(n_A, A_tri.data(), lda, uplo); + SpMatA A_sparse(n_A, n_A); + dense_to_sparse_format(Layout::ColMajor, A_tri.data(), T(0), A_sparse); + + // Random sparse B as a ColMajor dense buffer first, then convert to SpMatB. + std::vector B_dense(m_BY * n_BY, T(0)); + { + std::mt19937_64 rng(static_cast(seed_B)); + std::uniform_real_distribution uni01(0.0, 1.0); + std::uniform_real_distribution univ(-1.0, 1.0); + for (int64_t j = 0; j < n_BY; ++j) { + for (int64_t i = 0; i < m_BY; ++i) { + if (uni01(rng) < density_B) { + B_dense[i + j * m_BY] = static_cast(univ(rng)); + } + } + } + } + SpMatB B_sparse(m_BY, n_BY); + dense_to_sparse_format(Layout::ColMajor, B_dense.data(), T(0), B_sparse); + + int64_t ldb = (layout == Layout::ColMajor) ? m_BY : n_BY; + int64_t ldy = ldb; + // The dense reference call to blas::symm uses the requested layout. + std::vector B_dense_layout(m_BY * n_BY); + if (layout == Layout::ColMajor) { + std::copy(B_dense.begin(), B_dense.end(), B_dense_layout.begin()); + } else { + for (int64_t i = 0; i < m_BY; ++i) + for (int64_t j = 0; j < n_BY; ++j) + B_dense_layout[i * ldb + j] = B_dense[i + j * m_BY]; + } + + std::vector Y_actual(m_BY * n_BY); + DenseDist DY(m_BY, n_BY, ScalarDist::Uniform); + RandBLAS::fill_dense_unpacked(layout, DY, m_BY, n_BY, 0, 0, Y_actual.data(), RNGState(seed_B + 13)); + std::vector Y_expect = Y_actual; + + // Reference: dense blas::symm on full-storage A and dense B. + blas::symm(layout, side, uplo, m_BY, n_BY, + alpha, A_full.data(), lda, B_dense_layout.data(), ldb, + beta, Y_expect.data(), ldy); + + // Under test: sparse-symm A times sparse B via Case D. + RandBLAS::sparse_data::spsymm(layout, side, uplo, m_BY, n_BY, + alpha, A_sparse, B_sparse, + beta, Y_actual.data(), ldy); + + T atol = T(100) * std::numeric_limits::epsilon(); + T rtol = T(10) * std::numeric_limits::epsilon(); + auto msg = RandBLAS::testing::matrices_approx_equal( + layout, blas::Op::NoTrans, m_BY, n_BY, + Y_actual.data(), ldy, Y_expect.data(), ldy, + __PRETTY_FUNCTION__, __FILE__, __LINE__, atol, rtol + ); + if (!msg.empty()) FAIL() << msg; + } }; @@ -203,21 +283,36 @@ TEST_F(TestSpsymm, CSR_Left_Alpha0) { run_case>(Layout::ColMajor, Side::Left, Uplo::Upper, 10, 4, 0.0, 0.5, 0, 3); } -// Symmetric wrapper routing: covers the public RandBLAS::spsymm(layout, Symmetric, ...) overload -// Case D stub: sparse-symmetric A times sparse B must throw RandBLAS::Error. -// The API surface is reserved; the body is "not implemented" pending a future PR. -TEST_F(TestSpsymm, CaseD_SparseSparseThrows) { - CSRMatrix A_sparse(4, 4); - CSRMatrix B_sparse(4, 3); - std::vector Y(4 * 3, 0.0); - EXPECT_THROW({ - RandBLAS::sparse_data::spsymm( - Layout::ColMajor, Side::Left, Uplo::Upper, - 4, 3, 1.0, - A_sparse, B_sparse, - 0.0, Y.data(), 4 - ); - }, RandBLAS::Error); + +// Format-pair sweep: 3 A-formats x 3 B-formats x both sides + an uplo and +// edge-case sample. MKL handles all 9 pairs via make_mkl_handle. +TEST_F(TestSpsymm, CaseD_CSR_CSR_Left) { + run_case_d, CSRMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 1.5, -0.5, 0, 1); + run_case_d, CSRMatrix>(Layout::RowMajor, Side::Left, Uplo::Lower, 8, 3, 1.5, -0.5, 0, 1); +} +TEST_F(TestSpsymm, CaseD_CSC_CSC_Left) { + run_case_d, CSCMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 1.5, -0.5, 0, 1); + run_case_d, CSCMatrix>(Layout::RowMajor, Side::Left, Uplo::Lower, 8, 3, 1.5, -0.5, 0, 1); +} +TEST_F(TestSpsymm, CaseD_COO_COO_Left) { + run_case_d, COOMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 1.5, -0.5, 0, 1); + run_case_d, COOMatrix>(Layout::RowMajor, Side::Left, Uplo::Lower, 8, 3, 1.5, -0.5, 0, 1); +} +TEST_F(TestSpsymm, CaseD_Mixed_Format_Left) { + run_case_d, CSCMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 1.0, 0.0, 0, 1); + run_case_d, CSRMatrix>(Layout::ColMajor, Side::Left, Uplo::Lower, 8, 3, 1.0, 0.0, 0, 1); + run_case_d, CSRMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 1.0, 0.0, 0, 1); +} +TEST_F(TestSpsymm, CaseD_CSR_CSR_Right) { + run_case_d, CSRMatrix>(Layout::ColMajor, Side::Right, Uplo::Upper, 8, 3, 1.5, -0.5, 0, 1); + run_case_d, CSRMatrix>(Layout::RowMajor, Side::Right, Uplo::Lower, 8, 3, 1.5, -0.5, 0, 1); +} +TEST_F(TestSpsymm, CaseD_Float) { + run_case_d, CSRMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 1.5f, -0.5f, 0, 1); +} +TEST_F(TestSpsymm, CaseD_AlphaZero_BetaScale) { + // alpha=0 path: just beta-scales Y, doesn't even touch A or B. + run_case_d, CSRMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 0.0, 0.5, 0, 1); } // Routes through the public RandBLAS::spsymm(Symmetric) wrapper From 0af25b0a44c791b6365ffbe0cc17ea6abf45ae92 Mon Sep 17 00:00:00 2001 From: mmelnich Date: Fri, 15 May 2026 16:06:04 -0400 Subject: [PATCH 18/19] Address Riley review pass: factor COO body, drop perf claim, fix stale 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. --- RandBLAS/sksy.hh | 99 +++---------- RandBLAS/sparse_data/DevNotes.md | 16 ++- RandBLAS/sparse_data/coo_sksys_impl.hh | 184 +++++++++++++++++++++++++ rtd/source/FAQ.rst | 4 +- 4 files changed, 214 insertions(+), 89 deletions(-) create mode 100644 RandBLAS/sparse_data/coo_sksys_impl.hh diff --git a/RandBLAS/sksy.hh b/RandBLAS/sksy.hh index 59846919..68ae37a1 100644 --- a/RandBLAS/sksy.hh +++ b/RandBLAS/sksy.hh @@ -32,14 +32,17 @@ #include "RandBLAS/util.hh" #include "RandBLAS/base.hh" #include "RandBLAS/skge.hh" +#include "RandBLAS/sparse_data/coo_sksys_impl.hh" // ============================================================================= // Symmetric sketching helpers (SYMM-backed). See project-plans/randblas-symm-plan.md // for the four-case API design and the implementation status of each case. // - lsksy3, rsksy3: Case A (dense-symm A x dense Omega), via blas::symm. -// - lsksys, rsksys: Case B (dense-symm A x sparse SkOp), per-stored-entry -// two-axpy scatter that reads only the uplo triangle of A. +// - lsksys, rsksys: Case B (dense-symm A x sparse SkOp). Thin wrappers +// handling the SparseSkOp materialization-and-recurse; the actual COO +// walk + two-axpy scatter lives in sparse_data/coo_sksys_impl.hh +// as coo_lsksys / coo_rsksys. // ============================================================================= namespace RandBLAS::dense { @@ -105,7 +108,7 @@ void lsksy3( } else { // Layout mismatch: transpose-copy S into a tight buffer in the // caller's layout, then SYMM. Costs O(d*n) for the copy, but keeps - // the SYMM speedup (~1.3-1.8x over GEMM) on the d*n*n_A matvec. + // the SYMM speedup over GEMM on the d*n*n_A matvec. int64_t lds_new = (layout == blas::Layout::ColMajor) ? d : n; std::vector S_copy(static_cast(d) * static_cast(n)); auto [irs_in, ics_in] = layout_to_strides(S.layout, lds); @@ -227,46 +230,10 @@ void lsksys( } util::lascl(layout, d, n, beta, B, ldb); - if (alpha == T(0)) return; - auto Scoo = coo_view_of_skop(S); - const bool col_major = (layout == blas::Layout::ColMajor); - - for (int64_t p = 0; p < Scoo.nnz; ++p) { - sint_t row_S = Scoo.rows[p]; - sint_t col_S = Scoo.cols[p]; - if (row_S < ro_s || row_S >= ro_s + d) continue; - if (col_S < co_s || col_S >= co_s + n) continue; - int64_t i = static_cast(row_S) - ro_s; // B row - int64_t j = static_cast(col_S) - co_s; // A row index (sym A read as row j) - T av = alpha * Scoo.vals[p]; - - // B[i, :] += av * row_j(A_sym), split by uplo into two contiguous - // ranges of A so each becomes a single axpy. - if (uplo == blas::Uplo::Upper) { - // row j of A: - // c in [0, j]: read A[c, j] (column j of stored Upper, rows 0..j) - // c in (j, n-1]: read A[j, c] (row j of stored Upper, cols j+1..n-1) - if (col_major) { - blas::axpy(j + 1, av, &A[j * lda], 1, &B[i], ldb); - blas::axpy(n - j - 1, av, &A[j + (j + 1) * lda], lda, &B[i + (j + 1) * ldb], ldb); - } else { - blas::axpy(j + 1, av, &A[j], lda, &B[i * ldb], 1); - blas::axpy(n - j - 1, av, &A[j * lda + j + 1], 1, &B[i * ldb + j + 1], 1); - } - } else { - // Lower-stored: - // c in [0, j]: read A[j, c] (row j of stored Lower, cols 0..j) - // c in (j, n-1]: read A[c, j] (column j of stored Lower, rows j+1..n-1) - if (col_major) { - blas::axpy(j + 1, av, &A[j], lda, &B[i], ldb); - blas::axpy(n - j - 1, av, &A[(j + 1) + j * lda], 1, &B[i + (j + 1) * ldb], ldb); - } else { - blas::axpy(j + 1, av, &A[j * lda], 1, &B[i * ldb], 1); - blas::axpy(n - j - 1, av, &A[(j + 1) * lda + j], lda, &B[i * ldb + j + 1], 1); - } - } - } + RandBLAS::sparse_data::coo_lsksys( + layout, uplo, d, n, alpha, Scoo, ro_s, co_s, A, lda, B, ldb + ); } @@ -305,45 +272,10 @@ void rsksys( } util::lascl(layout, n, d, beta, B, ldb); - if (alpha == T(0)) return; - auto Scoo = coo_view_of_skop(S); - const bool col_major = (layout == blas::Layout::ColMajor); - - for (int64_t p = 0; p < Scoo.nnz; ++p) { - sint_t row_S = Scoo.rows[p]; - sint_t col_S = Scoo.cols[p]; - if (row_S < ro_s || row_S >= ro_s + n) continue; - if (col_S < co_s || col_S >= co_s + d) continue; - int64_t i = static_cast(row_S) - ro_s; // A column index - int64_t j = static_cast(col_S) - co_s; // B column - T av = alpha * Scoo.vals[p]; - - // B[:, j] += av * col_i(A_sym), split by uplo: - if (uplo == blas::Uplo::Upper) { - // col i of A: - // r in [0, i]: read A[r, i] (column i of stored Upper, rows 0..i) - // r in (i, n-1]: read A[i, r] (row i of stored Upper, cols i+1..n-1) - if (col_major) { - blas::axpy(i + 1, av, &A[i * lda], 1, &B[j * ldb], 1); - blas::axpy(n - i - 1, av, &A[i + (i + 1) * lda], lda, &B[(i + 1) + j * ldb], 1); - } else { - blas::axpy(i + 1, av, &A[i], lda, &B[j], ldb); - blas::axpy(n - i - 1, av, &A[i * lda + i + 1], 1, &B[(i + 1) * ldb + j], ldb); - } - } else { - // Lower-stored col i: - // r in [0, i-1]: read A[i, r] (row i of stored Lower, cols 0..i-1) - // r in [i, n-1]: read A[r, i] (column i of stored Lower, rows i..n-1) - if (col_major) { - blas::axpy(i, av, &A[i], lda, &B[j * ldb], 1); - blas::axpy(n - i, av, &A[i + i * lda], 1, &B[i + j * ldb], 1); - } else { - blas::axpy(i, av, &A[i * lda], 1, &B[j], ldb); - blas::axpy(n - i, av, &A[i * lda + i], lda, &B[i * ldb + j], ldb); - } - } - } + RandBLAS::sparse_data::coo_rsksys( + layout, uplo, n, d, alpha, A, lda, Scoo, ro_s, co_s, B, ldb + ); } } // end namespace RandBLAS::sparse @@ -414,9 +346,10 @@ using namespace RandBLAS::sparse; /// /// S - [in] /// * A SketchingOperator object (DenseSkOp or SparseSkOp). -/// * 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`. +/// * Defines :math:`\submat(\mtxS).` DenseSkOp dispatches to a SYMM-backed +/// kernel (Case A); SparseSkOp dispatches to a hand-rolled +/// per-stored-nonzero scatter that reads only the named triangle of +/// :math:`A` (Case B). See `project-plans/randblas-symm-plan.md`. /// /// ro_s - [in] /// * A nonnegative integer. diff --git a/RandBLAS/sparse_data/DevNotes.md b/RandBLAS/sparse_data/DevNotes.md index b971038a..7b15c958 100644 --- a/RandBLAS/sparse_data/DevNotes.md +++ b/RandBLAS/sparse_data/DevNotes.md @@ -63,7 +63,7 @@ two operands (``A`` symmetric vs. the second factor ``B``): | Tag | Operation | A storage | B storage | Status this PR | |-----|---------------------------------|---------------------|-----------|-------------------------------------------------| | 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). | +| B | dense-symm × sparse | dense, one triangle | sparse | Implemented via hand-rolled ``lsksys`` / ``rsksys`` wrappers in ``sksy.hh`` that handle the SparseSkOp materialization-and-recurse, and the ``coo_lsksys`` / ``coo_rsksys`` COO kernels in ``sparse_data/coo_sksys_impl.hh`` that do the 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``. | @@ -72,7 +72,7 @@ two operands (``A`` symmetric vs. the second factor ``B``): | Tag | MKL native? | Notes | |-----|----------------|-----------------------------------------------------------------------------------------------------------------| | A | No (BLAS++) | ``blas::symm`` directly. | -| B | No | The transpose trick puts the sparse op on the left of ``mkl_sparse_d_mm``, but the dense A has no ``matrix_descr``, so MKL can't be told A is symmetric. We hand-roll instead --- see ``lsksys`` / ``rsksys`` in ``sksy.hh``. | +| B | No | The transpose trick puts the sparse op on the left of ``mkl_sparse_d_mm``, but the dense A has no ``matrix_descr``, so MKL can't be told A is symmetric. We hand-roll instead --- see ``coo_lsksys`` / ``coo_rsksys`` in ``sparse_data/coo_sksys_impl.hh`` (with thin SparseSkOp wrappers in ``sksy.hh``). | | C | Yes (mostly) | ``mkl_sparse_d_mm`` with ``descr.type = SPARSE_MATRIX_TYPE_SYMMETRIC``. RandBLAS falls back to a hand kernel for side=Right (MKL has no Side parameter) and for CSC (``mkl_sparse_d_mm`` returns NOT_SUPPORTED on CSC). | | D | No | ``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. Symmetric expansion has to happen on the RandBLAS side, so we don't gain anything by routing through MKL. | @@ -107,9 +107,17 @@ The public-facing wrappers in the top-level ``RandBLAS::`` namespace are: -- routes via the ``Symmetric`` carrier so the uplo annotation travels with the matrix. -### Case B: hand-rolled in ``lsksys`` / ``rsksys`` +### Case B: hand-rolled COO kernels in ``coo_sksys_impl.hh`` + +The COO walk and scatter live in ``sparse_data/coo_sksys_impl.hh`` as +``coo_lsksys`` / ``coo_rsksys``, taking a ``COOMatrix`` plus +submatrix offsets and writing into a dense buffer. Thin wrappers +``lsksys`` / ``rsksys`` in ``sksy.hh`` handle the SparseSkOp +materialization-and-recurse pattern, beta-scale ``B`` via ``util::lascl``, +unpack the SparseSkOp into its COO view, and forward into the kernel. +The split keeps ``sksy.hh`` focused on SkOp dispatch and puts the +format-specific work next to the other COO kernels. -Lives in ``sksy.hh`` next to the dense-SkOp helpers ``lsksy3`` / ``rsksy3``. For each stored nonzero ``(i_S, j_S, v)`` of the SparseSkOp's COO view (with submatrix offsets ``(ro_s, co_s)`` filtered inline), the kernel applies ``alpha * v`` to one row or column of the symmetric dense A and accumulates diff --git a/RandBLAS/sparse_data/coo_sksys_impl.hh b/RandBLAS/sparse_data/coo_sksys_impl.hh new file mode 100644 index 00000000..cd7f5ebc --- /dev/null +++ b/RandBLAS/sparse_data/coo_sksys_impl.hh @@ -0,0 +1,184 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include "RandBLAS/base.hh" +#include "RandBLAS/util.hh" +#include "RandBLAS/sparse_data/base.hh" +#include "RandBLAS/sparse_data/coo_matrix.hh" + +#include + + +namespace RandBLAS::sparse_data { + +// ============================================================================= +// COO sparse-from-left times dense symmetric A into dense B: +// B = alpha * submat(Scoo) * sym(A, uplo) + beta * B +// +// - sym(A, uplo) is n-by-n dense symmetric, only the `uplo` triangle stored. +// - submat(Scoo) is the d-by-n window of Scoo at (ro_s, co_s). +// - B is d-by-n dense, layout-matched. +// +// For each stored nonzero (row_S, col_S, v) of Scoo inside the window, the +// kernel contributes alpha*v * (row col_S of A) to row (row_S - ro_s) of B. +// The row of symmetric A is assembled from the stored triangle as two +// blas::axpy calls -- one along the stored column up to (and including) the +// diagonal, one along the stored row past the diagonal -- chosen by `uplo` +// and the matrix layout. +// +// Beta-scaling of B and the alpha==0 short-circuit are the caller's +// responsibility (so that this kernel can be composed with other accumulators +// without redundant scaling). +// ============================================================================= +template +void coo_lsksys( + blas::Layout layout, + blas::Uplo uplo, + int64_t d, + int64_t n, + T alpha, + const COOMatrix &Scoo, + int64_t ro_s, + int64_t co_s, + const T *A, + int64_t lda, + T *B, + int64_t ldb +) { + if (alpha == T(0)) return; + + const bool col_major = (layout == blas::Layout::ColMajor); + + for (int64_t p = 0; p < Scoo.nnz; ++p) { + sint_t row_S = Scoo.rows[p]; + sint_t col_S = Scoo.cols[p]; + if (row_S < ro_s || row_S >= ro_s + d) continue; + if (col_S < co_s || col_S >= co_s + n) continue; + int64_t i = static_cast(row_S) - ro_s; // B row + int64_t j = static_cast(col_S) - co_s; // A row index (sym A read as row j) + T av = alpha * Scoo.vals[p]; + + // B[i, :] += av * row_j(A_sym), split by uplo into two contiguous + // ranges of A so each becomes a single axpy. + if (uplo == blas::Uplo::Upper) { + // row j of A: + // c in [0, j]: read A[c, j] (column j of stored Upper, rows 0..j) + // c in (j, n-1]: read A[j, c] (row j of stored Upper, cols j+1..n-1) + if (col_major) { + blas::axpy(j + 1, av, &A[j * lda], 1, &B[i], ldb); + blas::axpy(n - j - 1, av, &A[j + (j + 1) * lda], lda, &B[i + (j + 1) * ldb], ldb); + } else { + blas::axpy(j + 1, av, &A[j], lda, &B[i * ldb], 1); + blas::axpy(n - j - 1, av, &A[j * lda + j + 1], 1, &B[i * ldb + j + 1], 1); + } + } else { + // Lower-stored: + // c in [0, j]: read A[j, c] (row j of stored Lower, cols 0..j) + // c in (j, n-1]: read A[c, j] (column j of stored Lower, rows j+1..n-1) + if (col_major) { + blas::axpy(j + 1, av, &A[j], lda, &B[i], ldb); + blas::axpy(n - j - 1, av, &A[(j + 1) + j * lda], 1, &B[i + (j + 1) * ldb], ldb); + } else { + blas::axpy(j + 1, av, &A[j * lda], 1, &B[i * ldb], 1); + blas::axpy(n - j - 1, av, &A[(j + 1) * lda + j], lda, &B[i * ldb + j + 1], 1); + } + } + } +} + + +// ============================================================================= +// COO sparse-from-right times dense symmetric A into dense B: +// B = alpha * sym(A, uplo) * submat(Scoo) + beta * B +// +// Same two-axpy scatter pattern as coo_lsksys, but on columns of A: the +// column of A indexed by row_S of S contributes into the column of B indexed +// by col_S - co_s. +// +// Beta-scaling of B and the alpha==0 short-circuit are the caller's +// responsibility. +// ============================================================================= +template +void coo_rsksys( + blas::Layout layout, + blas::Uplo uplo, + int64_t n, + int64_t d, + T alpha, + const T *A, + int64_t lda, + const COOMatrix &Scoo, + int64_t ro_s, + int64_t co_s, + T *B, + int64_t ldb +) { + if (alpha == T(0)) return; + + const bool col_major = (layout == blas::Layout::ColMajor); + + for (int64_t p = 0; p < Scoo.nnz; ++p) { + sint_t row_S = Scoo.rows[p]; + sint_t col_S = Scoo.cols[p]; + if (row_S < ro_s || row_S >= ro_s + n) continue; + if (col_S < co_s || col_S >= co_s + d) continue; + int64_t i = static_cast(row_S) - ro_s; // A column index + int64_t j = static_cast(col_S) - co_s; // B column + T av = alpha * Scoo.vals[p]; + + // B[:, j] += av * col_i(A_sym), split by uplo: + if (uplo == blas::Uplo::Upper) { + // col i of A: + // r in [0, i]: read A[r, i] (column i of stored Upper, rows 0..i) + // r in (i, n-1]: read A[i, r] (row i of stored Upper, cols i+1..n-1) + if (col_major) { + blas::axpy(i + 1, av, &A[i * lda], 1, &B[j * ldb], 1); + blas::axpy(n - i - 1, av, &A[i + (i + 1) * lda], lda, &B[(i + 1) + j * ldb], 1); + } else { + blas::axpy(i + 1, av, &A[i], lda, &B[j], ldb); + blas::axpy(n - i - 1, av, &A[i * lda + i + 1], 1, &B[(i + 1) * ldb + j], ldb); + } + } else { + // Lower-stored col i: + // r in [0, i-1]: read A[i, r] (row i of stored Lower, cols 0..i-1) + // r in [i, n-1]: read A[r, i] (column i of stored Lower, rows i..n-1) + if (col_major) { + blas::axpy(i, av, &A[i], lda, &B[j * ldb], 1); + blas::axpy(n - i, av, &A[i + i * lda], 1, &B[i + j * ldb], 1); + } else { + blas::axpy(i, av, &A[i * lda], 1, &B[j], ldb); + blas::axpy(n - i, av, &A[i * lda + i], lda, &B[i * ldb + j], ldb); + } + } + } +} + +} // end namespace RandBLAS::sparse_data diff --git a/rtd/source/FAQ.rst b/rtd/source/FAQ.rst index 169d25dc..f4bed3bb 100644 --- a/rtd/source/FAQ.rst +++ b/rtd/source/FAQ.rst @@ -118,8 +118,8 @@ sketch_symmetric supports both DenseSkOp and SparseSkOp. reading 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`` falls back to GEMM. - When the ``DenseSkOp``'s storage layout differs from the caller's ``layout`` parameter, ``sketch_symmetric`` falls back to ``blas::gemm`` with the transpose flag --- ``blas::symm`` has no on-the-fly transpose flag for the dense operand. The layout-matched case still gets the SYMM speedup. +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. Language interoperability From e7986dcd3694d2ce2a7677c92a3af84aa06419f7 Mon Sep 17 00:00:00 2001 From: mmelnich Date: Fri, 15 May 2026 16:23:09 -0400 Subject: [PATCH 19/19] docs: strip em-dashes from comments and notes in symm-kernel files 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 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. --- RandBLAS/sparse_data/DevNotes.md | 8 ++++---- RandBLAS/sparse_data/coo_sksys_impl.hh | 6 +++--- rtd/source/api_reference/sketch_sparse.rst | 4 ++-- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/RandBLAS/sparse_data/DevNotes.md b/RandBLAS/sparse_data/DevNotes.md index 7b15c958..3ba86cb7 100644 --- a/RandBLAS/sparse_data/DevNotes.md +++ b/RandBLAS/sparse_data/DevNotes.md @@ -72,7 +72,7 @@ two operands (``A`` symmetric vs. the second factor ``B``): | Tag | MKL native? | Notes | |-----|----------------|-----------------------------------------------------------------------------------------------------------------| | A | No (BLAS++) | ``blas::symm`` directly. | -| B | No | The transpose trick puts the sparse op on the left of ``mkl_sparse_d_mm``, but the dense A has no ``matrix_descr``, so MKL can't be told A is symmetric. We hand-roll instead --- see ``coo_lsksys`` / ``coo_rsksys`` in ``sparse_data/coo_sksys_impl.hh`` (with thin SparseSkOp wrappers in ``sksy.hh``). | +| B | No | The transpose trick puts the sparse op on the left of ``mkl_sparse_d_mm``, but the dense A has no ``matrix_descr``, so MKL can't be told A is symmetric. We hand-roll instead. See ``coo_lsksys`` / ``coo_rsksys`` in ``sparse_data/coo_sksys_impl.hh`` (with thin SparseSkOp wrappers in ``sksy.hh``). | | C | Yes (mostly) | ``mkl_sparse_d_mm`` with ``descr.type = SPARSE_MATRIX_TYPE_SYMMETRIC``. RandBLAS falls back to a hand kernel for side=Right (MKL has no Side parameter) and for CSC (``mkl_sparse_d_mm`` returns NOT_SUPPORTED on CSC). | | D | No | ``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. Symmetric expansion has to happen on the RandBLAS side, so we don't gain anything by routing through MKL. | @@ -101,10 +101,10 @@ the ``Y <- beta * Y`` pass on entry. The public-facing wrappers in the top-level ``RandBLAS::`` namespace are: - - ``spsymm(layout, uplo, m, n, alpha, A, B, ldb, beta, Y, ldy)`` -- + - ``spsymm(layout, uplo, m, n, alpha, A, B, ldb, beta, Y, ldy)``, convenience for side=Left. - - ``spsymm(layout, m, n, alpha, Symmetric A_sym, B, ldb, beta, Y, ldy)`` - -- routes via the ``Symmetric`` carrier so the uplo annotation + - ``spsymm(layout, m, n, alpha, Symmetric A_sym, B, ldb, beta, Y, ldy)``, + routes via the ``Symmetric`` carrier so the uplo annotation travels with the matrix. ### Case B: hand-rolled COO kernels in ``coo_sksys_impl.hh`` diff --git a/RandBLAS/sparse_data/coo_sksys_impl.hh b/RandBLAS/sparse_data/coo_sksys_impl.hh index cd7f5ebc..7aa82dd2 100644 --- a/RandBLAS/sparse_data/coo_sksys_impl.hh +++ b/RandBLAS/sparse_data/coo_sksys_impl.hh @@ -50,9 +50,9 @@ namespace RandBLAS::sparse_data { // For each stored nonzero (row_S, col_S, v) of Scoo inside the window, the // kernel contributes alpha*v * (row col_S of A) to row (row_S - ro_s) of B. // The row of symmetric A is assembled from the stored triangle as two -// blas::axpy calls -- one along the stored column up to (and including) the -// diagonal, one along the stored row past the diagonal -- chosen by `uplo` -// and the matrix layout. +// blas::axpy calls: one along the stored column up to (and including) the +// diagonal, one along the stored row past the diagonal. The split is chosen +// by `uplo` and the matrix layout. // // Beta-scaling of B and the alpha==0 short-circuit are the caller's // responsibility (so that this kernel can be composed with other accumulators diff --git a/rtd/source/api_reference/sketch_sparse.rst b/rtd/source/api_reference/sketch_sparse.rst index 8505c694..51890b21 100644 --- a/rtd/source/api_reference/sketch_sparse.rst +++ b/rtd/source/api_reference/sketch_sparse.rst @@ -150,8 +150,8 @@ Deterministic operations the densified buffer. This covers all 3 x 3 = 9 sparse-format pairings for ``(A, B)``. ``mkl_sparse_sp2m`` rejects symmetric descriptors and ``mkl_sparse_d_spmmd`` takes no descriptor, so - routing through MKL would not avoid the symmetric expansion --- - composing through Case C costs an ``O(m*n)`` temporary, which is + routing through MKL would not avoid the symmetric expansion. The + composition through Case C costs an ``O(m*n)`` temporary, which is small for the typical RandNLA workload where ``B`` is a sketching operator with ``nnz(B) << m*n``.