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/RandBLAS/sksy.hh b/RandBLAS/sksy.hh index e49208b6..1f0ac639 100644 --- a/RandBLAS/sksy.hh +++ b/RandBLAS/sksy.hh @@ -32,162 +32,310 @@ #include "RandBLAS/util.hh" #include "RandBLAS/base.hh" #include "RandBLAS/skge.hh" +#include "RandBLAS/sparse_data/coo_sksys_impl.hh" -namespace RandBLAS { +#include +#include -using namespace RandBLAS::dense; -using namespace RandBLAS::sparse; +// ============================================================================= +// Symmetric sketching helpers (SYMM-backed). See sparse_data/DevNotes.md for +// the four-case design (dense/sparse symmetric operand x dense/sparse factor). +// - lsksy3, rsksy3: Case A (dense-symm A x dense Omega), via blas::symm. +// - lsksys, rsksys: Case B (dense-symm A x sparse SkOp). Thin wrappers +// handling validation, beta, and SparseSkOp materialization; the actual +// column-driven accumulation kernel lives in +// sparse_data/coo_sksys_impl.hh as coo_lsksys / coo_rsksys. +// ============================================================================= -// MARK: SUBMAT(S) +namespace RandBLAS::dense { // ============================================================================= -/// \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).` -/// -/// ldb - [in] -/// * A nonnegative integer. -/// * Leading dimension of :math:`\mat(B)` when reading from :math:`B.` -/// -/// @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. +// Copy the n_rows-by-n_cols matrix at src (read in src_layout with leading +// dimension ld_src) into a tight buffer laid out in target_layout, whose +// leading dimension is written to ld_out. Used by lsksy3 / rsksy3 when the +// SkOp's storage layout mismatches the caller's: blas::symm cannot transpose +// an operand on the fly, so the copy keeps the SYMM speedup at an O(size) +// one-time cost. +// ============================================================================= +template +inline std::vector transpose_copy_to_layout( + blas::Layout target_layout, int64_t n_rows, int64_t n_cols, + const T* src, blas::Layout src_layout, int64_t ld_src, + int64_t &ld_out +) { + ld_out = (target_layout == blas::Layout::ColMajor) ? n_rows : n_cols; + std::vector out(static_cast(n_rows) * static_cast(n_cols)); + auto [irs_in, ics_in] = layout_to_strides(src_layout, ld_src); + auto [irs_out, ics_out] = layout_to_strides(target_layout, ld_out); + util::omatcopy(n_rows, n_cols, src, irs_in, ics_in, out.data(), irs_out, ics_out); + return out; +} + +// ============================================================================= +// 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, 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, - int64_t n, // number of rows in B - int64_t d, // number of columns in B + 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 T* A, + const DenseSkOp &S, + int64_t ro_s, + int64_t co_s, + const T *A, int64_t lda, - const SKOP &S, + 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 ); + validate_submat_dims(S.n_rows, S.n_cols, d, n, ro_s, 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 { + int64_t lds_new; + auto S_copy = transpose_copy_to_layout(layout, d, n, S_ptr, S.layout, lds, lds_new); + blas::symm(layout, blas::Side::Right, uplo, d, n, + alpha, A, lda, S_copy.data(), lds_new, beta, B, ldb); + } + return; +} + + +// ============================================================================= +// RSKSY3: SYMM-backed right-sketch with a symmetric matrix A. +// +// 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. +// +// Same materialization and layout-mismatch fallback semantics as `lsksy3`. +// Final call (matching layout): `blas::symm` with `side = Left`. +template +void rsksy3( + blas::Layout layout, + 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, + int64_t lda, + const DenseSkOp &S, int64_t ro_s, int64_t co_s, T beta, - T* B, - int64_t ldb, - T sym_check_tol = 0 + 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, n, d, n, alpha, A, lda, S, ro_s, co_s, beta, B, 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 ); + validate_submat_dims(S.n_rows, S.n_cols, n, d, ro_s, 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 { + int64_t lds_new; + auto S_copy = transpose_copy_to_layout(layout, n, d, S_ptr, S.layout, lds, lds_new); + blas::symm(layout, blas::Side::Left, uplo, n, d, + alpha, A, lda, S_copy.data(), lds_new, beta, B, ldb); + } + return; } +} // end namespace RandBLAS::dense + + +namespace RandBLAS::sparse { // ============================================================================= -/// \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 -/// ) +// 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. +// +// Validation mirrors the dense-path lsksy3. When S is unmaterialized, only +// the requested d-by-n window is sampled (submatrix_as_coo, the same pattern +// lskges uses); a materialized S is consumed through a lightweight COO view +// with the window filtered inside the kernel. The kernel itself +// (sparse_data::coo_lsksys) is a column-driven pure accumulator; beta is +// applied here, exactly once. +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 +) { + validate_submat_dims(S.n_rows, S.n_cols, d, n, ro_s, co_s); + if (layout == blas::Layout::ColMajor) { + randblas_require(lda >= n); + randblas_require(ldb >= d); + } else { + randblas_require(lda >= n); + randblas_require(ldb >= n); + } + + util::lascl(layout, d, n, beta, B, ldb); + if (alpha == T(0)) return; + + if (S.nnz < 0) { + // Sample only the requested window rather than materializing all of S. + auto Ssub = submatrix_as_coo(S, d, n, ro_s, co_s); + RandBLAS::sparse_data::coo_lsksys( + layout, uplo, d, n, alpha, Ssub, 0, 0, A, lda, B, ldb + ); + return; + } + auto Scoo = coo_view_of_skop(S); + RandBLAS::sparse_data::coo_lsksys( + layout, uplo, d, n, alpha, Scoo, ro_s, co_s, A, lda, B, ldb + ); +} + + +// ============================================================================= +// 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 validation, window-sampling, and beta conventions as lsksys; the +// kernel (sparse_data::coo_rsksys) reduces to coo_lsksys via the transpose +// identity. +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 +) { + validate_submat_dims(S.n_rows, S.n_cols, n, d, ro_s, co_s); + if (layout == blas::Layout::ColMajor) { + randblas_require(lda >= n); + randblas_require(ldb >= n); + } else { + randblas_require(lda >= n); + randblas_require(ldb >= d); + } + + util::lascl(layout, n, d, beta, B, ldb); + if (alpha == T(0)) return; + + if (S.nnz < 0) { + // Sample only the requested window rather than materializing all of S. + auto Ssub = submatrix_as_coo(S, n, d, ro_s, co_s); + RandBLAS::sparse_data::coo_rsksys( + layout, uplo, n, d, alpha, A, lda, Ssub, 0, 0, B, ldb + ); + return; + } + auto Scoo = coo_view_of_skop(S); + RandBLAS::sparse_data::coo_rsksys( + layout, uplo, n, d, alpha, A, lda, Scoo, ro_s, co_s, B, ldb + ); +} + +} // end namespace RandBLAS::sparse + + +namespace RandBLAS { + +using namespace RandBLAS::dense; +using namespace RandBLAS::sparse; + + +// MARK: SUBMAT(S) + +// ============================================================================= +/// \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 +343,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. +/// 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. /// -/// **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 +/// 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 +372,27 @@ 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).` DenseSkOp dispatches to a SYMM-backed +/// kernel (Case A); SparseSkOp dispatches to a column-driven +/// accumulation kernel that reads only the named triangle of +/// :math:`A`. See ``RandBLAS/sparse_data/DevNotes.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,249 +406,239 @@ 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 ) { - 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); + if constexpr (requires { S.buff; S.layout; }) { + RandBLAS::dense::lsksy3(layout, uplo, d, n, alpha, S, ro_s, co_s, A, lda, beta, B, ldb); + } else if constexpr (requires { S.nnz; }) { + RandBLAS::sparse::lsksys(layout, uplo, d, n, alpha, S, ro_s, co_s, A, lda, beta, B, ldb); + } else { + static_assert(sizeof(SKOP) == 0, + "sketch_symmetric supports DenseSkOp and SparseSkOp. For other " + "SketchingOperator types, apply the operator to a fully-stored A " + "with sketch_general."); + // see GitHub PR #155 for why we don't use static_assert(false, ...). + } } -// 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 -/// -/// .. 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). +/// Sketch from the right in a SYMM-like operation /// -/// 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) = \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$} /// -/// .. math:: -/// \mat(B)_{ij} = B[i \cdot \ldb + j]. +/// where :math:`\alpha` and :math:`\beta` are real scalars and :math:`\mtxS` is a sketching operator. /// -/// In this case, :math:`\ldb` must be :math:`\geq d.` +/// 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).` /// -/// .. dropdown:: Full parameter descriptions -/// :animate: fade-in-slide-down +/// @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 +) { + if constexpr (requires { S.buff; S.layout; }) { + RandBLAS::dense::rsksy3(layout, uplo, n, d, alpha, A, lda, S, ro_s, co_s, beta, B, ldb); + } else if constexpr (requires { S.nnz; }) { + RandBLAS::sparse::rsksys(layout, uplo, n, d, alpha, A, lda, S, ro_s, co_s, beta, B, ldb); + } else { + static_assert(sizeof(SKOP) == 0, + "sketch_symmetric supports DenseSkOp and SparseSkOp. For other " + "SketchingOperator types, apply the operator to a fully-stored A " + "with sketch_general."); + // see GitHub PR #155 for why we don't use static_assert(false, ...). + } +} + + +// 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). /// -/// layout - [in] -/// * Either Layout::ColMajor or Layout::RowMajor -/// * Matrix storage for :math:`\mat(B).` +/// .. 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$} /// -/// alpha - [in] -/// * A real scalar. -/// * If zero, then :math:`A` is not accessed. +/// 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}`). /// -/// A - [in] -/// * Pointer to a 1D array of real scalars. -/// * Defines :math:`\mat(A).` +/// See the submatrix overload above for the meaning of :math:`\mat(A)` and +/// :math:`\texttt{uplo}`. /// -/// lda - [in] -/// * A nonnegative integer. -/// * Leading dimension of :math:`\mat(A)` when reading from :math:`A.` -/// -/// S - [in] -/// * A DenseSkOp or SparseSkOp object. +/// @endverbatim +template +inline void 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 +) { + int64_t d = S.dist.n_rows; + int64_t n = S.dist.n_cols; + sketch_symmetric(layout, uplo, d, n, alpha, S, 0, 0, A, lda, beta, B, ldb); +} + + +// ============================================================================= +/// \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 +/// Sketch from the right in a SYMM-like operation, with :math:`\mtxS` used in full. /// -/// beta - [in] -/// * A real scalar. -/// * If zero, then :math:`B` need not be set on input. +/// .. 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$} /// -/// 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:`n` and :math:`d` are taken from :math:`\mtxS` directly +/// (:math:`n = \mtxS.\ttt{dist.n\_rows}`, :math:`d = \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 T* A, int64_t lda, const SKOP &S, T beta, - T* B, - int64_t ldb, - T sym_check_tol = 0 + T* B, int64_t ldb ) { int64_t n = S.dist.n_rows; int64_t d = S.dist.n_cols; + sketch_symmetric(layout, uplo, n, d, alpha, A, lda, S, 0, 0, beta, B, ldb); +} + + +// MARK: LEGACY (sym_check_tol) + +// ============================================================================= +// The four overloads below reproduce the pre-Uplo sketch_symmetric API and +// semantics exactly: mat(A) must be stored in the format of a general matrix +// with BOTH triangles populated (both are read), a runtime symmetry check +// runs first (skip it, at your own peril, by passing sym_check_tol < 0), and +// the operation forwards to sketch_general. They are retained so that +// existing call sites keep compiling and behaving identically; new code +// should prefer the blas::Uplo overloads above, which read only the named +// triangle and skip the O(n^2) runtime check. +// ============================================================================= + +// ============================================================================= +/// \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 +/// ) +/// @verbatim embed:rst:leading-slashes +/// Legacy overload, retained for API compatibility. Check that :math:`\mat(A)` +/// is symmetric up to tolerance :math:`\texttt{sym_check_tol}` (pass a negative +/// tolerance to skip the check), then sketch from the left: +/// :math:`\mat(B) = \alpha \cdot \submat(\mtxS) \cdot \mat(A) + \beta \cdot \mat(B)`. +/// Requires both triangles of :math:`\mat(A)` populated; both are read. +/// Prefer the ``blas::Uplo`` overloads for new code. +/// @endverbatim +template +inline void sketch_symmetric( + // B = alpha*S*A + beta*B + 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 +) { 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); + sketch_general(layout, blas::Op::NoTrans, blas::Op::NoTrans, d, n, n, alpha, S, ro_s, co_s, A, lda, beta, B, ldb); } +// ============================================================================= +/// \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 +/// Legacy overload, retained for API compatibility. Same contract as its +/// left-sketch sibling above, sketching from the right: +/// :math:`\mat(B) = \alpha \cdot \mat(A) \cdot \submat(\mtxS) + \beta \cdot \mat(B)`. +/// @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. + 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 +) { + 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); +} // ============================================================================= /// \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 -/// ) +/// ) /// @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. -/// -/// **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.` -/// -/// 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).` -/// -/// ldb - [in] -/// * A nonnegative integer. -/// * Leading dimension of :math:`\mat(B)` when reading from :math:`B.` -/// +/// Legacy overload, retained for API compatibility; :math:`\mtxS` used in +/// full. Same contract as the submatrix legacy overloads. /// @endverbatim -template +template inline void sketch_symmetric( // B = alpha*S*A + beta*B blas::Layout layout, T alpha, const SKOP &S, - const T* A, - int64_t lda, + const T* A, int64_t lda, T beta, - T* B, - int64_t ldb, + T* B, int64_t ldb, T sym_check_tol = 0 ) { int64_t d = S.dist.n_rows; @@ -529,4 +647,30 @@ inline void sketch_symmetric( sketch_general(layout, blas::Op::NoTrans, blas::Op::NoTrans, d, n, n, alpha, S, 0, 0, A, lda, beta, B, ldb); } +// ============================================================================= +/// \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 +/// ) +/// @verbatim embed:rst:leading-slashes +/// Legacy overload, retained for API compatibility; :math:`\mtxS` used in +/// full, sketching from the right. Same contract as the submatrix legacy +/// overloads. +/// @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. + 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 +) { + 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); +} + } // end namespace RandBLAS diff --git a/RandBLAS/sparse_data/DevNotes.md b/RandBLAS/sparse_data/DevNotes.md index c575abf2..88a09913 100644 --- a/RandBLAS/sparse_data/DevNotes.md +++ b/RandBLAS/sparse_data/DevNotes.md @@ -60,6 +60,134 @@ we're inside one of those functions. This heuristic can differ from that used if we had called `[left/right]_spmm` directly on a COO matrix. +## 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 | +|-----|---------------------------------|---------------------|-----------|-------------------------------------------------| +| 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 ``lsksys`` / ``rsksys`` wrappers in ``sksy.hh`` (validation, beta, window sampling) over the column-driven ``coo_lsksys`` kernel in ``sparse_data/coo_sksys_impl.hh`` (``coo_rsksys`` reduces to it via the transpose identity). Only the named triangle of A is read. | +| C | sparse-symm × dense (→ dense) | sparse, one triangle | dense | Implemented in ``spsymm_dispatch.hh`` (side=Right normalized at entry; MKL fast path covers all three formats; column-driven per-format fallbacks). | +| D | sparse-symm × sparse → dense | sparse, one triangle | sparse | Implemented in ``spsymm_dispatch.hh``: expand A's triangle to a general sparse matrix and reuse the sparse-times-sparse path (MKL builds); densify-B + Case-C composition as the non-MKL / index-width-mismatch fallback. | + +### 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. 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 | ``mkl_sparse_?_mm`` with ``descr.type = SPARSE_MATRIX_TYPE_SYMMETRIC``. side=Right is normalized to side=Left by the dispatcher before MKL is reached (layout-flip identity, valid since ``A == A^T``); CSC is consumed inside ``mkl_spsymm`` as a CSR-of-transpose view with ``uplo`` flipped. The hand kernels run only on non-MKL builds, on index-width mismatch with ``MKL_INT``, or on a runtime ``NOT_SUPPORTED``. | +| D | Yes (via expansion) | ``mkl_sparse_sp2m`` returns ``SPARSE_STATUS_NOT_SUPPORTED`` when ``descrA.type == SPARSE_MATRIX_TYPE_SYMMETRIC`` (only ``GENERAL`` is accepted there), and ``mkl_sparse_?_spmmd`` takes no descriptor at all, so the symmetric expansion happens on the RandBLAS side: ``expand_symmetric_to_general`` (O(nnz)) followed by ``mkl_spgemm_to_dense`` with a GENERAL descriptor. | + +### Case C dispatch (``spsymm_dispatch.hh``) + +``RandBLAS::sparse_data::spsymm(layout, side, uplo, m, n, alpha, A, B, ldb, beta, C, ldc)`` +dispatches as follows: + +1. Normalize side=Right to side=Left: since ``A == A^T``, ``C = B*A`` equals + ``C^T = A*B^T``, and reinterpreting the B and C buffers in the opposite + layout presents them as ``B^T`` and ``C^T`` with the same leading + dimensions. ``uplo`` is unchanged. Everything below assumes side=Left. +2. Validate: zero-based indices (``A.index_base == IndexBase::Zero``), A + square of order ``m``, and leading-dimension lower bounds for B and C + (mirroring ``left_spmm``). +3. Handle empty products (a zero dimension, alpha == 0, or a structurally + empty operand) by leaving beta * C, before MKL can reject a valid empty + sparse matrix at handle creation (the left_spmm contract). +4. If RandBLAS was built with MKL and the index width matches ``MKL_INT``, + try ``mkl::mkl_spsymm`` with the caller's beta (MKL applies alpha and + beta itself, as ``left_spmm`` hands beta to ``mkl_left_spmm``). It applies + the ``SPARSE_MATRIX_TYPE_SYMMETRIC`` descriptor and calls + ``mkl_sparse_?_mm``; CSC goes through the CSR-of-transpose view with + ``uplo`` flipped. It returns false only on a runtime ``NOT_SUPPORTED`` (a + parameter-validation result, so C is untouched when it happens); control + then falls to step 5. +5. Format-specific fallback: apply beta via ``util::lascl``, then run + ``csr_spsymm`` / ``coo_spsymm`` (and + ``csc_spsymm``, a three-line delegation to ``csr_spsymm`` on the + transpose view with ``uplo`` flipped). The kernels are column-driven: + an OpenMP-parallel outer loop over the n right-hand-side columns of B/C + (each column is owned by one thread, so there are no races), with a scan + of the stored triangle inside. Each stored off-diagonal entry + ``A(i,j) = v`` contributes twice per column (the entry and its implied + symmetric counterpart); diagonal entries once. Entries outside the named + triangle are silently skipped, so a caller that mistakenly stored both + triangles still gets the correct answer. + +The public-facing wrappers in the top-level ``RandBLAS::`` namespace are: + + - ``spsymm(layout, uplo, n, alpha, A, B, ldb, beta, C, ldc)``, + convenience for side=Left. The order of A comes from the matrix itself + (A must be square), so the single dimension argument is ``n``, the + column count of B and C -- the same convention as MKL's + ``mkl_sparse_?_mm``. + - ``spmm(layout, n, alpha, Symmetric A_sym, B, ldb, beta, C, ldc)``, + an overload of the general ``spmm``: the ``Symmetric`` carrier + requests symmetric semantics through the type system, so the uplo + annotation travels with the matrix. + +### Case B: column-driven COO kernel in ``coo_sksys_impl.hh`` + +The accumulation kernel lives in ``sparse_data/coo_sksys_impl.hh`` as +``coo_lsksys``, taking a ``COOMatrix`` plus submatrix offsets and +writing into a dense buffer. ``coo_rsksys`` is a delegation: ``B = A*S`` +implies ``B^T = S^T * A`` (A symmetric), so it calls ``coo_lsksys`` with the +layout and ``uplo`` flipped, the ``Scoo.transpose()`` view, and the window +offsets swapped. Wrappers ``lsksys`` / ``rsksys`` in ``sksy.hh`` validate +the arguments (mirroring the dense-path ``lsksy3`` / ``rsksy3``), beta-scale +``B`` via ``util::lascl`` exactly once, 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. + +The kernel is column-driven: an OpenMP-parallel outer loop over the n +columns of B (each column owned by one thread, race-free), with a scan of +the window nonzeros of S inside. The contribution of nonzero +``(row_S, col_S, v)`` to column c is +``B(row_S - ro_s, c) += alpha * v * sym(A, uplo)(col_S - co_s, c)``, where +the symmetric element read resolves to the stored triangle by swapping the +index pair when it falls outside it. There is one address computation for A +(a two-way in-triangle test) instead of a per-``uplo``-per-layout grid of +strided AXPY range splits. SparseSkOp is COO internally, so submatrix +filtering is a direct ``if (row < ro_s ...) continue`` on the COO triples. + +For an unmaterialized SparseSkOp, ``lsksys`` / ``rsksys`` sample only the +requested window via ``submatrix_as_coo`` (the same pattern ``lskges`` +uses), so the RNG and memory cost is proportional to the window, not to the +full operator; a materialized SparseSkOp is consumed through +``coo_view_of_skop`` with the window filtered inside the kernel. + +### Case D: expand-A + spgemm reuse (densify-B as the fallback) + +Lives in ``spsymm_dispatch.hh`` next to the Case-C dispatcher. side=Right is +normalized to side=Left at entry via the same layout-flip identity as +Case C, with ``B^T`` obtained as the lightweight ``B.transpose()`` view. +Validation and the single beta application follow the Case-C pattern. + +Why MKL cannot be handed the symmetric A directly: + + - ``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_?_spmmd`` (which writes directly to dense ``C``) + accepts no descriptor at all. + +So the symmetric expansion happens on the RandBLAS side. Primary path +(MKL builds with both index widths matching ``MKL_INT``): +``expand_symmetric_to_general(A, uplo)`` (in ``symmetric.hh``) builds an +owning general COO from the stored triangle in ``O(nnz)`` memory, and +``mkl_spgemm_to_dense`` runs on it with a GENERAL descriptor. ``B`` stays +sparse, so the cost is proportional to the actual nonzero structure. + +Fallback path (non-MKL builds, or index width mismatched with +``MKL_INT``): densify ``B`` into an ``m`` by ``n`` temporary in the +caller's ``layout`` via the format-specific ``coo_to_dense`` / +``csr_to_dense`` / ``csc_to_dense`` helper, then compose through the +Case-C overload. The ``O(m*n)`` temporary is a fallback-only cost. + ## 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/base.hh b/RandBLAS/sparse_data/base.hh index ecf98471..ed687d9c 100644 --- a/RandBLAS/sparse_data/base.hh +++ b/RandBLAS/sparse_data/base.hh @@ -32,6 +32,7 @@ #include "RandBLAS/base.hh" #include #include +#include // std::iota, used below #ifdef __cpp_concepts #include diff --git a/RandBLAS/sparse_data/coo_matrix.hh b/RandBLAS/sparse_data/coo_matrix.hh index 33b15f88..11b0e44a 100644 --- a/RandBLAS/sparse_data/coo_matrix.hh +++ b/RandBLAS/sparse_data/coo_matrix.hh @@ -380,8 +380,8 @@ void dense_to_coo(Layout layout, T* mat, T abs_tol, COOMatrix &spmat) { } } -template -void coo_to_dense(const COOMatrix &spmat, int64_t stride_row, int64_t stride_col, T *mat) { +template +void coo_to_dense(const COOMatrix &spmat, int64_t stride_row, int64_t stride_col, T *mat) { #define MAT(_i, _j) mat[(_i) * stride_row + (_j) * stride_col] for (int64_t i = 0; i < spmat.n_rows; ++i) { for (int64_t j = 0; j < spmat.n_cols; ++j) { @@ -404,8 +404,8 @@ void coo_to_dense(const COOMatrix &spmat, int64_t stride_row, int64_t stride_ return; } -template -void coo_to_dense(const COOMatrix &spmat, Layout layout, T *mat) { +template +void coo_to_dense(const COOMatrix &spmat, Layout layout, T *mat) { if (layout == Layout::ColMajor) { coo_to_dense(spmat, 1, spmat.n_rows, mat); } else { diff --git a/RandBLAS/sparse_data/coo_sksys_impl.hh b/RandBLAS/sparse_data/coo_sksys_impl.hh new file mode 100644 index 00000000..32d0fd58 --- /dev/null +++ b/RandBLAS/sparse_data/coo_sksys_impl.hh @@ -0,0 +1,147 @@ +// 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) + 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. +// +// The loop is column-driven: each column c of B accumulates, over the window +// nonzeros (row_S, col_S, v) of Scoo, the contribution +// B(row_S - ro_s, c) += alpha * v * sym(A, uplo)(col_S - co_s, c), +// where the symmetric element read resolves to the stored triangle by +// swapping the index pair when it falls outside it. Columns are independent, +// so the outer loop parallelizes without races, and the accesses to A and B +// walk down single columns in ColMajor. +// +// 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; + + // Plain locals rather than structured bindings: clang does not (yet) + // support referencing structured bindings inside OpenMP regions. + stride_64t sa = layout_to_strides(layout, lda); + stride_64t sb = layout_to_strides(layout, ldb); + int64_t irs_a = sa.inter_row_stride, ics_a = sa.inter_col_stride; + int64_t irs_b = sb.inter_row_stride, ics_b = sb.inter_col_stride; + bool upper = (uplo == blas::Uplo::Upper); + + #pragma omp parallel for schedule(static) + for (int64_t c = 0; c < n; ++c) { + T* B_c = &B[c * ics_b]; + for (int64_t p = 0; p < Scoo.nnz; ++p) { + int64_t row_S = (int64_t) Scoo.rows[p]; + int64_t col_S = (int64_t) 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 = row_S - ro_s; // B row + int64_t r = col_S - co_s; // row of sym(A) contributing to column c + // sym(A, uplo)(r, c): read (r, c) if it lies in the stored + // triangle, else the mirrored (c, r). + bool mirrored = upper ? (r > c) : (r < c); + T a_rc = mirrored ? A[c * irs_a + r * ics_a] + : A[r * irs_a + c * ics_a]; + B_c[i * irs_b] += alpha * Scoo.vals[p] * a_rc; + } + } +} + + +// ============================================================================= +// COO sparse-from-right times dense symmetric A into dense B: +// B = alpha * sym(A, uplo) * submat(Scoo) + B +// +// - sym(A, uplo) is n-by-n; submat(Scoo) is the n-by-d window at +// (ro_s, co_s); B is n-by-d. +// +// Reduction to coo_lsksys via the transpose identity: B = A * S implies +// B^T = S^T * A^T = S^T * A (A symmetric). Reading the A and B buffers in +// the flipped layout presents them as A^T (= A, with the stored triangle +// name flipped) and B^T; S^T is the lightweight Scoo.transpose() view, with +// the window offsets swapped. +// +// 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 +) { + auto flipped_uplo = (uplo == blas::Uplo::Upper) + ? blas::Uplo::Lower + : blas::Uplo::Upper; + auto St = Scoo.transpose(); + coo_lsksys(flipped_layout(layout), flipped_uplo, d, n, alpha, + St, co_s, ro_s, A, lda, B, ldb); +} + +} // end namespace RandBLAS::sparse_data diff --git a/RandBLAS/sparse_data/coo_spsymm_impl.hh b/RandBLAS/sparse_data/coo_spsymm_impl.hh new file mode 100644 index 00000000..1a1e1ea2 --- /dev/null +++ b/RandBLAS/sparse_data/coo_spsymm_impl.hh @@ -0,0 +1,88 @@ +// 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 + +namespace RandBLAS::sparse_data { + +// ============================================================================= +// COO fallback for symmetric sparse-times-dense (side=Left). +// +// Accumulates C += alpha * A * B over the stored (row, col, val) triples; +// no assumption on their order. For uplo=Upper, structurally populated +// entries satisfy row <= col; for uplo=Lower, row >= col. Entries outside +// the named triangle are silently skipped. The caller (the spsymm +// dispatcher) has already validated the arguments, applied beta scaling to +// C, and normalized side=Right away, so this kernel is a pure accumulator. +// +// Column-driven for the same reasons as csr_spsymm: the outer loop over +// dense right-hand-side columns is race-free under OpenMP and the inner +// updates are unit-stride in ColMajor. +template +void coo_spsymm( + blas::Layout layout, + blas::Uplo uplo, + int64_t m, int64_t n, + T alpha, + const COOMatrix& A, + const T* B, int64_t ldb, + T* C, int64_t ldc +) { + (void) m; + // Plain locals rather than structured bindings: clang does not (yet) + // support referencing structured bindings inside OpenMP regions. + stride_64t sb = layout_to_strides(layout, ldb); + stride_64t sc = layout_to_strides(layout, ldc); + int64_t irs_b = sb.inter_row_stride, ics_b = sb.inter_col_stride; + int64_t irs_c = sc.inter_row_stride, ics_c = sc.inter_col_stride; + bool upper = (uplo == blas::Uplo::Upper); + + #pragma omp parallel for schedule(static) + for (int64_t c = 0; c < n; ++c) { + const T* B_c = &B[c * ics_b]; + T* C_c = &C[c * ics_c]; + for (int64_t p = 0; p < A.nnz; ++p) { + int64_t i = (int64_t) A.rows[p]; + int64_t j = (int64_t) A.cols[p]; + if (upper ? (j < i) : (j > i)) continue; + T av = alpha * A.vals[p]; + C_c[i * irs_c] += av * B_c[j * irs_b]; + if (i != j) + C_c[j * irs_c] += av * B_c[i * irs_b]; + } + } +} + +} // namespace RandBLAS::sparse_data diff --git a/RandBLAS/sparse_data/csc_matrix.hh b/RandBLAS/sparse_data/csc_matrix.hh index 7789b52c..913da69a 100644 --- a/RandBLAS/sparse_data/csc_matrix.hh +++ b/RandBLAS/sparse_data/csc_matrix.hh @@ -262,8 +262,8 @@ namespace RandBLAS::sparse_data::csc { using namespace RandBLAS::sparse_data; using blas::Layout; -template -void csc_to_dense(const CSCMatrix &spmat, int64_t stride_row, int64_t stride_col, T *mat) { +template +void csc_to_dense(const CSCMatrix &spmat, int64_t stride_row, int64_t stride_col, T *mat) { randblas_require(spmat.index_base == IndexBase::Zero); #define MAT(_i, _j) mat[(_i) * stride_row + (_j) * stride_col] for (int64_t i = 0; i < spmat.n_rows; ++i) { @@ -282,8 +282,8 @@ void csc_to_dense(const CSCMatrix &spmat, int64_t stride_row, int64_t stride_ return; } -template -void csc_to_dense(const CSCMatrix &spmat, Layout layout, T *mat) { +template +void csc_to_dense(const CSCMatrix &spmat, Layout layout, T *mat) { if (layout == Layout::ColMajor) { csc_to_dense(spmat, 1, spmat.n_rows, mat); } else { diff --git a/RandBLAS/sparse_data/csc_spsymm_impl.hh b/RandBLAS/sparse_data/csc_spsymm_impl.hh new file mode 100644 index 00000000..5eb142a8 --- /dev/null +++ b/RandBLAS/sparse_data/csc_spsymm_impl.hh @@ -0,0 +1,63 @@ +// 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/csc_matrix.hh" +#include "RandBLAS/sparse_data/csr_spsymm_impl.hh" +#include + +namespace RandBLAS::sparse_data { + +// ============================================================================= +// CSC fallback for symmetric sparse-times-dense (side=Left). +// +// A symmetric CSC matrix's arrays, read as CSR, describe A^T = A over the +// same buffers, with the stored triangle name flipped: a CSC Upper entry at +// (i, j) with i <= j appears in the CSR view at (j, i), which is Lower. +// So this kernel is a delegation to csr_spsymm on the lightweight +// A.transpose() view with uplo flipped. Same identity as the MKL path. +template +void csc_spsymm( + blas::Layout layout, + blas::Uplo uplo, + int64_t m, int64_t n, + T alpha, + const CSCMatrix& A, + const T* B, int64_t ldb, + T* C, int64_t ldc +) { + auto At = A.transpose(); + blas::Uplo uplo_flipped = (uplo == blas::Uplo::Upper) + ? blas::Uplo::Lower + : blas::Uplo::Upper; + csr_spsymm(layout, uplo_flipped, m, n, alpha, At, B, ldb, C, ldc); +} + +} // namespace RandBLAS::sparse_data diff --git a/RandBLAS/sparse_data/csr_matrix.hh b/RandBLAS/sparse_data/csr_matrix.hh index 026e1b66..b2224f3a 100644 --- a/RandBLAS/sparse_data/csr_matrix.hh +++ b/RandBLAS/sparse_data/csr_matrix.hh @@ -263,8 +263,8 @@ namespace RandBLAS::sparse_data::csr { using namespace RandBLAS::sparse_data; using blas::Layout; -template -void csr_to_dense(const CSRMatrix &spmat, int64_t stride_row, int64_t stride_col, T *mat) { +template +void csr_to_dense(const CSRMatrix &spmat, int64_t stride_row, int64_t stride_col, T *mat) { randblas_require(spmat.index_base == IndexBase::Zero); auto rowptr = spmat.rowptr; auto colidxs = spmat.colidxs; @@ -286,8 +286,8 @@ void csr_to_dense(const CSRMatrix &spmat, int64_t stride_row, int64_t stride_ return; } -template -void csr_to_dense(const CSRMatrix &spmat, Layout layout, T *mat) { +template +void csr_to_dense(const CSRMatrix &spmat, Layout layout, T *mat) { if (layout == Layout::ColMajor) { csr_to_dense(spmat, 1, spmat.n_rows, mat); } else { diff --git a/RandBLAS/sparse_data/csr_spsymm_impl.hh b/RandBLAS/sparse_data/csr_spsymm_impl.hh new file mode 100644 index 00000000..46b95b4f --- /dev/null +++ b/RandBLAS/sparse_data/csr_spsymm_impl.hh @@ -0,0 +1,93 @@ +// 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/csr_matrix.hh" +#include + +namespace RandBLAS::sparse_data { + + +// ============================================================================= +// CSR fallback for the symmetric sparse-times-dense kernel (side=Left). +// +// Accumulates C += alpha * A * B, where A is m-by-m symmetric with only the +// triangle named by uplo structurally stored, and B, C are dense m-by-n. +// The caller (the spsymm dispatcher) has already validated the arguments, +// applied beta scaling to C, and normalized side=Right away, so this kernel +// is a pure accumulator. Entries 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). +// +// The loop is column-driven: each dense right-hand-side column c of B/C is +// processed independently, so the outer loop parallelizes without races and +// the inner updates are unit-stride in ColMajor. Each stored off-diagonal +// entry A(i, j) = v contributes twice per column (the entry itself and the +// implied symmetric A(j, i) = v); diagonal entries contribute once. +template +void csr_spsymm( + blas::Layout layout, + blas::Uplo uplo, + int64_t m, int64_t n, + T alpha, + const CSRMatrix& A, + const T* B, int64_t ldb, + T* C, int64_t ldc +) { + (void) m; + // Plain locals rather than structured bindings: clang does not (yet) + // support referencing structured bindings inside OpenMP regions. + stride_64t sb = layout_to_strides(layout, ldb); + stride_64t sc = layout_to_strides(layout, ldc); + int64_t irs_b = sb.inter_row_stride, ics_b = sb.inter_col_stride; + int64_t irs_c = sc.inter_row_stride, ics_c = sc.inter_col_stride; + bool upper = (uplo == blas::Uplo::Upper); + + #pragma omp parallel for schedule(static) + for (int64_t c = 0; c < n; ++c) { + const T* B_c = &B[c * ics_b]; + T* C_c = &C[c * ics_c]; + for (int64_t i = 0; i < A.n_rows; ++i) { + for (int64_t p = A.rowptr[i]; p < A.rowptr[i+1]; ++p) { + int64_t j = (int64_t) A.colidxs[p]; + if (upper ? (j < i) : (j > i)) continue; + T av = alpha * A.vals[p]; + C_c[i * irs_c] += av * B_c[j * irs_b]; + if (i != j) + C_c[j * irs_c] += av * B_c[i * irs_b]; + } + } + } +} + +} // namespace RandBLAS::sparse_data diff --git a/RandBLAS/sparse_data/mkl_spmm_impl.hh b/RandBLAS/sparse_data/mkl_spmm_impl.hh index 12c0f5a3..3a77e818 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" @@ -257,6 +258,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. @@ -334,25 +367,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( - mkl_op, 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( - mkl_op, alpha, h.handle, descr, - to_mkl_layout(layout), - B, (MKL_INT)n, (MKL_INT)ldb, - beta, C, (MKL_INT)ldc - ); - } else { - // unsupported floating point type. - return false; - } + sparse_status_t status = mkl_sparse_mm_call( + mkl_op, 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 } @@ -394,20 +413,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; } @@ -437,17 +444,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/sparse_data/mkl_spsymm_impl.hh b/RandBLAS/sparse_data/mkl_spsymm_impl.hh new file mode 100644 index 00000000..9d90b84a --- /dev/null +++ b/RandBLAS/sparse_data/mkl_spsymm_impl.hh @@ -0,0 +1,128 @@ +// 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: C = alpha * A * B + beta * C, side=Left by +// contract. The spsymm dispatcher normalizes side=Right to side=Left (via the +// layout-flip identity) before this function is reached, and passes the +// caller's beta through: MKL applies alpha and beta itself, matching how +// left_spmm hands beta to mkl_left_spmm. +// A is symmetric sparse (one triangle stored, named by uplo); B and C dense. +// +// Returns true if MKL handled the call, false to signal fallback to the +// hand-rolled per-format kernel. Fallback triggers: +// - a runtime SPARSE_STATUS_NOT_SUPPORTED from mkl_sparse_?_mm (some MKL +// versions return it for combinations we could not predict). +// The dispatcher additionally skips this function entirely when the index +// type width does not match MKL_INT, or on non-MKL builds. +// +// 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, +// valid since A is symmetric so A == A^T) and recursing. The triangle the +// user named in the CSC is 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( + 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 *C, + int64_t ldc +) { + (void) m; + using sint_t = typename SpMat::index_t; + constexpr bool is_csc = std::is_same_v>; + + if constexpr (is_csc) { + auto At = A.transpose(); + blas::Uplo uplo_flipped = (uplo == blas::Uplo::Upper) + ? blas::Uplo::Lower + : blas::Uplo::Upper; + return mkl_spsymm(layout, uplo_flipped, m, n, + alpha, At, B, ldb, beta, C, ldc); + } + + 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 = mkl_sparse_mm_call( + SPARSE_OPERATION_NON_TRANSPOSE, alpha, h.handle, descr, + to_mkl_layout(layout), + B, n, ldb, beta, C, ldc + ); + + // Some MKL versions return NOT_SUPPORTED for combinations we could not + // predict. Do not throw; signal fallback. + if (status == SPARSE_STATUS_NOT_SUPPORTED) + return false; + + check_mkl_status(status, "mkl_sparse_mm (symmetric)"); + return true; +} + +} // namespace RandBLAS::sparse_data::mkl + +#endif // RandBLAS_HAS_MKL diff --git a/RandBLAS/sparse_data/sksp.hh b/RandBLAS/sparse_data/sksp.hh index 37984ce8..6a4fc5d5 100644 --- a/RandBLAS/sparse_data/sksp.hh +++ b/RandBLAS/sparse_data/sksp.hh @@ -33,6 +33,7 @@ #include "RandBLAS/dense_skops.hh" #include "RandBLAS/exceptions.hh" #include "RandBLAS/util.hh" +#include "RandBLAS/sparse_skops.hh" namespace RandBLAS::sparse_data { @@ -322,6 +323,163 @@ void rsksp3( return; } +// MARK: LSKSPS + +// ============================================================================= +/// \fn lsksps(blas::Layout layout, blas::Op opS, blas::Op opA, int64_t d, +/// int64_t n, int64_t m, T alpha, const SparseSkOp &S, int64_t ro_s, int64_t co_s, +/// const SpMat &A, int64_t ro_a, int64_t co_a, T beta, T *B, int64_t ldb +/// ) +/// @verbatim embed:rst:leading-slashes +/// Sketch from the left in an SpMM-like operation +/// +/// .. math:: +/// \mat(B) = \alpha \cdot \underbrace{\op(\submat(\mtxS))}_{d \times m} \cdot \underbrace{\op(\mtxA)}_{m \times n} + \beta \cdot \underbrace{\mat(B)}_{d \times n}, \tag{$\star$} +/// +/// where :math:`\alpha` and :math:`\beta` are real scalars, :math:`\op(\mtxX)` either returns a matrix :math:`\mtxX` +/// or its transpose, and both :math:`\mtxA` and :math:`\mtxS` are sparse. +/// @endverbatim +// +// The sparse-times-sparse backend (spgemm, MKL-backed) has no submatrix +// parameters for either operand and no op flag for its second operand, so +// this function differs from lsksp3 in two ways: A must be passed in full +// (ro_a == co_a == 0; spgemm has nothing to extract a submatrix into without +// first copying A, which this PR does not do), and op(A) is folded into a +// transpose *view* of A (A.transpose(), the same lightweight technique +// spsymm_dispatch.hh's Case D uses) rather than passed as a flag. S keeps +// full submatrix support: submatrix_as_coo / coo_view_of_skop already build +// exactly the owning or non-owning COOMatrix window spgemm's first operand +// needs, mirroring how LSKGES handles a SparseSkOp against a dense operand. +template +void lsksps( + blas::Layout layout, + blas::Op opS, + blas::Op opA, + int64_t d, // B is d-by-n + int64_t n, // op(A) is m-by-n + int64_t m, // op(submat(S)) is d-by-m + T alpha, + const SparseSkOp &S, + int64_t ro_s, + int64_t co_s, + const SpMat &A, + int64_t ro_a, + int64_t co_a, + T beta, + T *B, + int64_t ldb +) { + randblas_require(ro_a == 0 && co_a == 0); + randblas_require(A.index_base == IndexBase::Zero); + auto [rows_submat_S, cols_submat_S] = dims_before_op(d, m, opS); + validate_submat_dims(S.n_rows, S.n_cols, rows_submat_S, cols_submat_S, ro_s, co_s); + randblas_require(A.n_rows == m && A.n_cols == n); + if (layout == blas::Layout::ColMajor) { + randblas_require(ldb >= d); + } else { + randblas_require(ldb >= n); + } + + if (d == 0 || n == 0) + return; + bool full_operator = (S.n_rows == rows_submat_S && S.n_cols == cols_submat_S && ro_s == 0 && co_s == 0); + if (alpha == T(0) || A.nnz == 0) { + RandBLAS::util::lascl(layout, d, n, beta, B, ldb); + return; + } + + // A submatrix of S must be materialized to know its nnz (there is no + // cheap structural check for a windowed operator), so this second empty + // check runs after building Scoo rather than before. + COOMatrix Scoo = full_operator + ? RandBLAS::sparse::coo_view_of_skop(S) + : RandBLAS::sparse::submatrix_as_coo(S, rows_submat_S, cols_submat_S, ro_s, co_s); + if (Scoo.nnz == 0) { + RandBLAS::util::lascl(layout, d, n, beta, B, ldb); + return; + } + + if (opA == blas::Op::NoTrans) { + spgemm(layout, opS, alpha, Scoo, A, beta, B, ldb); + } else { + auto At = A.transpose(); + spgemm(layout, opS, alpha, Scoo, At, beta, B, ldb); + } + return; +} + +// MARK: RSKSPS + +// ============================================================================= +/// \fn rsksps(blas::Layout layout, blas::Op opA, blas::Op opS, int64_t m, +/// int64_t d, int64_t n, T alpha, const SpMat &A, int64_t ro_a, int64_t co_a, +/// const SparseSkOp &S, int64_t ro_s, int64_t co_s, T beta, T *B, int64_t ldb +/// ) +/// @verbatim embed:rst:leading-slashes +/// Sketch from the right in an SpMM-like operation +/// +/// .. math:: +/// \mat(B) = \alpha \cdot \underbrace{\op(\mtxA)}_{m \times n} \cdot \underbrace{\op(\submat(\mtxS))}_{n \times d} + \beta \cdot \underbrace{\mat(B)}_{m \times d}, \tag{$\star$} +/// +/// where :math:`\alpha` and :math:`\beta` are real scalars, :math:`\op(\mtxX)` either returns a matrix :math:`\mtxX` +/// or its transpose, and both :math:`\mtxA` and :math:`\mtxS` are sparse. Same restrictions as +/// lsksps: A must be passed in full, and op(submat(S)) is folded into a transpose view when needed. +/// @endverbatim +template +void rsksps( + blas::Layout layout, + blas::Op opA, + blas::Op opS, + int64_t m, // B is m-by-d + int64_t d, // op(submat(S)) is n-by-d + int64_t n, // op(A) is m-by-n + T alpha, + const SpMat &A, + int64_t ro_a, + int64_t co_a, + const SparseSkOp &S, + int64_t ro_s, + int64_t co_s, + T beta, + T *B, + int64_t ldb +) { + randblas_require(ro_a == 0 && co_a == 0); + randblas_require(A.index_base == IndexBase::Zero); + auto [rows_submat_S, cols_submat_S] = dims_before_op(n, d, opS); + validate_submat_dims(S.n_rows, S.n_cols, rows_submat_S, cols_submat_S, ro_s, co_s); + randblas_require(A.n_rows == m && A.n_cols == n); + if (layout == blas::Layout::ColMajor) { + randblas_require(ldb >= m); + } else { + randblas_require(ldb >= d); + } + + if (m == 0 || d == 0) + return; + if (alpha == T(0) || A.nnz == 0) { + RandBLAS::util::lascl(layout, m, d, beta, B, ldb); + return; + } + + bool full_operator = (S.n_rows == rows_submat_S && S.n_cols == cols_submat_S && ro_s == 0 && co_s == 0); + COOMatrix Scoo = full_operator + ? RandBLAS::sparse::coo_view_of_skop(S) + : RandBLAS::sparse::submatrix_as_coo(S, rows_submat_S, cols_submat_S, ro_s, co_s); + if (Scoo.nnz == 0) { + RandBLAS::util::lascl(layout, m, d, beta, B, ldb); + return; + } + + if (opS == blas::Op::NoTrans) { + spgemm(layout, opA, alpha, A, Scoo, beta, B, ldb); + } else { + auto St = Scoo.transpose(); + spgemm(layout, opA, alpha, A, St, beta, B, ldb); + } + return; +} + } // end namespace RandBLAS::sparse_data @@ -535,4 +693,70 @@ inline void sketch_sparse( return; } + +// MARK: SKSP overloads, sparse sketching operator + +// ============================================================================= +/// \fn sketch_sparse(blas::Layout layout, blas::Op opS, blas::Op opA, int64_t d, int64_t n, int64_t m, +/// T alpha, const SparseSkOp &S, int64_t ro_s, int64_t co_s, const SpMat &A, T beta, T *B, int64_t ldb +/// ) +/// @verbatim embed:rst:leading-slashes +/// Sketch a sparse matrix from the left with a sparse sketching operator. +/// Same equation as the DenseSkOp overload above. Requires Intel MKL (the +/// sparse-times-sparse backend has no non-MKL fallback) and :math:`\mtxA` in +/// full: a submatrix of :math:`\mtxA` is not supported through this +/// overload. :math:`\submat(\mtxS)` is fully supported. +/// @endverbatim +template +inline void sketch_sparse( + blas::Layout layout, + blas::Op opS, + blas::Op opA, + int64_t d, // B is d-by-n + int64_t n, // A is m-by-n + int64_t m, // op(submat(\mtxS)) is d-by-m + T alpha, + const SparseSkOp &S, + int64_t ro_s, + int64_t co_s, + const SpMat &A, + T beta, + T *B, + int64_t ldb +) { + sparse_data::lsksps(layout, opS, opA, d, n, m, alpha, S, ro_s, co_s, A, 0, 0, beta, B, ldb); + return; +} + +// ============================================================================= +/// \fn sketch_sparse(blas::Layout layout, blas::Op opA, blas::Op opS, int64_t m, int64_t d, int64_t n, +/// T alpha, const SpMat &A, const SparseSkOp &S, int64_t ro_s, int64_t co_s, T beta, T *B, int64_t ldb +/// ) +/// @verbatim embed:rst:leading-slashes +/// Sketch a sparse matrix from the right with a sparse sketching operator. +/// Same equation as the DenseSkOp overload above, and the same restrictions +/// as the left-sketching overload just above: requires Intel MKL, and +/// :math:`\mtxA` must be passed in full. +/// @endverbatim +template +inline void sketch_sparse( + blas::Layout layout, + blas::Op opA, + blas::Op opS, + int64_t m, // B is m-by-d + int64_t d, // op(submat(\mtxA)) is m-by-n + int64_t n, // op(submat(\mtxS)) is n-by-d + T alpha, + const SpMat &A, + const SparseSkOp &S, + int64_t ro_s, + int64_t co_s, + T beta, + T *B, + int64_t ldb +) { + sparse_data::rsksps(layout, opA, opS, m, d, n, alpha, A, 0, 0, S, ro_s, co_s, beta, B, ldb); + return; +} + } // end namespace RandBLAS diff --git a/RandBLAS/sparse_data/spmm_dispatch.hh b/RandBLAS/sparse_data/spmm_dispatch.hh index 8df989a5..50dbd40b 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" diff --git a/RandBLAS/sparse_data/spsymm_dispatch.hh b/RandBLAS/sparse_data/spsymm_dispatch.hh new file mode 100644 index 00000000..ea0fe957 --- /dev/null +++ b/RandBLAS/sparse_data/spsymm_dispatch.hh @@ -0,0 +1,325 @@ +// 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/util.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" +// Pulls in the two general RandBLAS::spmm overloads so any TU that sees the +// Symmetric overload below sees the full RandBLAS::spmm overload set, +// not just the one declared in this file. +#include "RandBLAS/sparse_data/spmm_dispatch.hh" +#include "RandBLAS/config.h" + +#if defined(RandBLAS_HAS_MKL) +#include "RandBLAS/sparse_data/mkl_spsymm_impl.hh" +#endif + +#include +#include + + +namespace RandBLAS::sparse_data { + +// ============================================================================= +/// Dispatched symmetric sparse-times-dense multiplication. +/// +/// @verbatim embed:rst:leading-slashes +/// Computes +/// +/// .. math:: +/// \mat(C) = \alpha \cdot \op_{\ttt{side}}(\mat(A), \mat(B)) + \beta \cdot \mat(C), +/// +/// where: +/// +/// - :math:`\mat(A)` is a sparse symmetric matrix stored in COO, CSR, or +/// CSC format, with zero-based indices. 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(C)` 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(C) = \alpha A B + \beta C`. +/// - For :math:`\ttt{side} = \ttt{Right}`, :math:`\mat(A)` is n-by-n and +/// the operation is :math:`\mat(C) = \alpha B A + \beta C`. +/// +/// Arguments are validated at entry. Empty products (a zero dimension, a +/// structurally empty :math:`\mat(A)`, or :math:`\alpha = 0`) leave +/// :math:`\beta \cdot \mat(C)`. +/// @endverbatim +// +// Dispatch mechanics (see also sparse_data/DevNotes.md): side=Right is +// normalized to side=Left at entry (A == A^T, so C = B*A is C^T = A*B^T with +// the B/C buffers reinterpreted in the flipped layout). On MKL builds with +// the index width matching MKL_INT, mkl_spsymm runs with the caller's beta +// (SPARSE_MATRIX_TYPE_SYMMETRIC descriptor; CSC consumed as a +// CSR-of-transpose view with uplo flipped). The per-format hand kernels are +// pure accumulators reached on non-MKL builds, index-width mismatch, or a +// runtime NOT_SUPPORTED from MKL; beta is applied by lascl just before them, +// mirroring left_spmm. +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* C, int64_t ldc +) { + using blas::Layout; + if (side == blas::Side::Right) { + // A == A^T, so C = alpha B A + beta C is C^T = alpha A B^T + beta C^T + // with B^T and C^T read from the same buffers in the flipped layout. + // The dimensions of C^T are n-by-m; uplo is unchanged (the equation + // transpose does not change which physical entries of A are stored). + spsymm(flipped_layout(layout), blas::Side::Left, uplo, n, m, alpha, A, B, ldb, beta, C, ldc); + return; + } + + 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."); + + // side is Left from here on: A is m-by-m, B and C are m-by-n. + randblas_require(A.index_base == IndexBase::Zero); + randblas_require(A.n_rows == A.n_cols); + randblas_require(A.n_rows == m); + if (layout == Layout::ColMajor) { + randblas_require(ldb >= m); + randblas_require(ldc >= m); + } else { + randblas_require(ldb >= n); + randblas_require(ldc >= n); + } + + // Empty products: no output elements, or nothing to accumulate beyond + // beta * C. Handled here because MKL rejects some valid empty sparse + // matrices at handle creation (same contract as left_spmm). + if (m == 0 || n == 0) + return; + if (alpha == T(0) || A.nnz == 0) { + RandBLAS::util::lascl(layout, m, n, beta, C, ldc); + return; + } + +#if defined(RandBLAS_HAS_MKL) + if constexpr (sizeof(sint_t) == sizeof(MKL_INT)) { + // MKL applies beta itself. + bool handled = mkl::mkl_spsymm( + layout, uplo, m, n, alpha, A, B, ldb, beta, C, ldc + ); + if (handled) return; + } +#endif + + // Fallback path: apply beta here, then run a pure-accumulator hand + // kernel. Safe after an MKL NOT_SUPPORTED, which is a parameter + // validation result: C is untouched when it fires. + RandBLAS::util::lascl(layout, m, n, beta, C, ldc); + if constexpr (is_csr) { + csr_spsymm(layout, uplo, m, n, alpha, A, B, ldb, C, ldc); + } else if constexpr (is_csc) { + csc_spsymm(layout, uplo, m, n, alpha, A, B, ldb, C, ldc); + } else if constexpr (is_coo) { + coo_spsymm(layout, uplo, m, n, alpha, A, B, ldb, C, ldc); + } +} + +// ============================================================================= +/// Case D: sparse-symmetric A times sparse B, dense output. +/// +/// @verbatim embed:rst:leading-slashes +/// Computes the same operation as the Case-C overload, with :math:`\mat(B)` +/// sparse (COO, CSR, or CSC; any index type). side=Right is normalized to +/// side=Left at entry via the same layout-flip identity, with +/// :math:`B^T` obtained as the lightweight ``B.transpose()`` view. +/// +/// Primary path (MKL builds with index widths matching MKL_INT): expand +/// :math:`A`'s stored triangle into a general sparse matrix in +/// :math:`O(\ttt{nnz})` memory (``expand_symmetric_to_general``), then call +/// the existing sparse-times-sparse dense-output routine +/// (``mkl_spgemm_to_dense``) with a GENERAL descriptor. The expansion has to +/// happen on the RandBLAS side either way: MKL's ``mkl_sparse_sp2m`` returns +/// ``SPARSE_STATUS_NOT_SUPPORTED`` when either operand's ``matrix_descr`` is +/// ``SPARSE_MATRIX_TYPE_SYMMETRIC``, and ``mkl_sparse_?_spmmd`` accepts no +/// descriptor at all. Keeping :math:`B` sparse keeps the cost proportional +/// to the actual nonzero structure. +/// +/// Fallback path (non-MKL builds, or index width mismatched with MKL_INT): +/// densify :math:`B` into a temporary dense buffer in the caller's layout +/// (cost: an :math:`O(m \cdot n)` temporary, a fallback-only property) and +/// compose through the Case-C overload above. +/// @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* C, int64_t ldc +) { + using blas::Layout; + static_assert(std::is_same_v, + "Case D: A and B must share scalar_t."); + + if (side == blas::Side::Right) { + // Same identity as Case C; B^T is a lightweight transpose view. + auto Bt = B.transpose(); + spsymm(flipped_layout(layout), blas::Side::Left, uplo, n, m, alpha, A, Bt, beta, C, ldc); + return; + } + + // side is Left from here on: A is m-by-m, B and C are m-by-n. + randblas_require(A.index_base == IndexBase::Zero); + randblas_require(B.index_base == IndexBase::Zero); + randblas_require(A.n_rows == A.n_cols); + randblas_require(A.n_rows == m); + randblas_require(B.n_rows == m); + randblas_require(B.n_cols == n); + if (layout == Layout::ColMajor) { + randblas_require(ldc >= m); + } else { + randblas_require(ldc >= n); + } + + // Empty products: no output elements, or nothing to accumulate beyond + // beta * C. Handled here because MKL rejects some valid empty sparse + // matrices at handle creation (same contract as left_spmm). + if (m == 0 || n == 0) + return; + if (alpha == T(0) || A.nnz == 0 || B.nnz == 0) { + RandBLAS::util::lascl(layout, m, n, beta, C, ldc); + return; + } + + using sint_A = typename SpMatA::index_t; + using sint_B = typename SpMatB::index_t; + +#if defined(RandBLAS_HAS_MKL) + if constexpr (sizeof(sint_A) == sizeof(MKL_INT) && sizeof(sint_B) == sizeof(MKL_INT)) { + // mkl_spgemm_to_dense applies beta itself. + auto A_general = expand_symmetric_to_general(A, uplo); + mkl::mkl_spgemm_to_dense( + layout, blas::Op::NoTrans, alpha, A_general, B, beta, C, ldc + ); + return; + } +#endif + + // Fallback: densify B into a tight buffer in the caller's layout and + // compose through Case C, which handles beta itself. + int64_t ldb_dense = (layout == Layout::ColMajor) ? m : n; + std::vector B_dense(static_cast(m) * static_cast(n), T(0)); + + 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."); + } + + spsymm(layout, blas::Side::Left, uplo, m, n, + alpha, A, B_dense.data(), ldb_dense, beta, C, ldc); +} + +} // end namespace RandBLAS::sparse_data + + +namespace RandBLAS { + +// ============================================================================= +/// Convenience wrapper for symmetric sparse-times-dense matmul. +/// +/// Computes \math{C = \alpha A B + \beta C}, where \math{A} is the symmetric +/// sparse matrix and only the triangle named by uplo is read. \math{A} must +/// be square; its order is taken from the matrix itself, so the only +/// dimension argument is \math{n}, the number of columns in \math{B} and +/// \math{C}. +template +inline void spsymm( + blas::Layout layout, + blas::Uplo uplo, + int64_t n, + T alpha, + const SpMat& A, + const T* B, int64_t ldb, + T beta, + T* C, int64_t ldc +) { + RandBLAS::sparse_data::spsymm( + layout, blas::Side::Left, uplo, A.n_rows, n, + alpha, A, B, ldb, beta, C, ldc + ); +} + +// ============================================================================= +/// Overload of spmm for a sparse matrix tagged as symmetric. +/// +/// Computes \math{C = \alpha A B + \beta C}, where the Symmetric wrapper +/// supplies the matrix and the triangle to read; the opposite triangle is +/// implied by symmetry. Same dimension convention as the spsymm overload +/// above: the only dimension argument is \math{n}, the number of columns in +/// \math{B} and \math{C}. +template +inline void spmm( + blas::Layout layout, + int64_t n, + T alpha, + const Symmetric& A_sym, + const T* B, int64_t ldb, + T beta, + T* C, int64_t ldc +) { + RandBLAS::spsymm(layout, A_sym.uplo, n, alpha, A_sym.A, B, ldb, beta, C, ldc); +} + +} // end namespace RandBLAS diff --git a/RandBLAS/sparse_data/symmetric.hh b/RandBLAS/sparse_data/symmetric.hh new file mode 100644 index 00000000..be2a67cd --- /dev/null +++ b/RandBLAS/sparse_data/symmetric.hh @@ -0,0 +1,207 @@ +// 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/sparse_data/coo_matrix.hh" +#include "RandBLAS/sparse_data/csr_matrix.hh" +#include "RandBLAS/sparse_data/csc_matrix.hh" +#include "RandBLAS/exceptions.hh" + +#include + + +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}`. +/// +/// Symmetric semantics on a sparse matrix are requested either by name +/// (``spsymm``, which takes the bare :math:`\ttt{SparseMatrix}` plus an +/// explicit ``uplo``) or, for the general product, by wrapping the matrix in +/// ``Symmetric`` and passing it to :math:`\ttt{spmm}` — an overload +/// that a bare :math:`\ttt{SparseMatrix}` cannot bind to, since ``spmm`` / +/// ``spgemm`` always treat a bare :math:`\ttt{SparseMatrix}` as general. The +/// wrapper type is intentionally separate from the :math:`\ttt{SparseMatrix}` +/// concept so passing a symmetric matrix to ``spmm`` without the wrapper +/// fails to compile rather than silently dropping the opposite triangle. +/// +/// The wrapper is non-owning: it holds a reference to :math:`A`, not a copy. +/// Keep the named object alive for the lifetime of the wrapper. Wrapping a +/// temporary (for example, the by-value view returned by ``.transpose()``) +/// is rejected at compile time via deleted rvalue overloads, since the +/// temporary would be destroyed at the end of the statement and leave the +/// wrapper dangling. +/// @endverbatim +template +struct Symmetric { + const SpMat& A; + const blas::Uplo uplo; + // Alias for A.n_rows (== A.n_cols, enforced below), matching the + // n_rows/n_cols convention DenseSkOp and SparseSkOp use for their own + // wrapped distributions. + const int64_t n_rows; + const int64_t n_cols; + + 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), n_rows(A_in.n_rows), n_cols(A_in.n_cols) { + randblas_require(A_in.n_rows == A_in.n_cols); + } + + // Binding a temporary would dangle at the end of the statement. The + // const&& form catches const and non-const rvalues alike (transpose() + // views are const prvalues), while lvalues still bind to the const& + // constructor above. + Symmetric(const SpMat&&, blas::Uplo) = delete; +}; + + +// ============================================================================= +/// 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); +} + +// Wrapping a temporary (e.g. as_symmetric(A.transpose(), uplo)) would leave +// the wrapper referencing a destroyed object; reject it at compile time. +// The const&& form catches const and non-const rvalues alike without the +// forwarding-reference trap that a plain && overload would create. +template +Symmetric as_symmetric(const SpMat&& A, blas::Uplo uplo) = delete; + + +// ============================================================================= +/// Expand the stored triangle of a symmetric sparse matrix into an owning +/// general (both triangles populated) COOMatrix. +/// +/// @verbatim embed:rst:leading-slashes +/// Reads only the triangle of :math:`A` named by :math:`\ttt{uplo}`; entries +/// outside it are skipped, matching the semantics of the spsymm kernels. +/// Each stored off-diagonal entry :math:`(i, j, v)` is emitted twice (as +/// :math:`(i, j, v)` and :math:`(j, i, v)`); diagonal entries once. Memory +/// cost is :math:`O(\ttt{nnz})`. +/// +/// This is the bridge from one-triangle symmetric storage to consumers that +/// only accept general sparse matrices (notably MKL's sparse-times-sparse +/// routines, which reject ``SPARSE_MATRIX_TYPE_SYMMETRIC`` descriptors). +/// @endverbatim +template +COOMatrix expand_symmetric_to_general(const SpMat& A, blas::Uplo uplo) { + randblas_require(A.n_rows == A.n_cols); + randblas_require(A.index_base == IndexBase::Zero); + + 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, + "expand_symmetric_to_general requires COO, CSR, or CSC."); + + bool upper = (uplo == blas::Uplo::Upper); + // in_triangle(i, j): the entry belongs to the named triangle. + auto in_triangle = [upper](int64_t i, int64_t j) { + return upper ? (j >= i) : (j <= i); + }; + + // Pass 1: count the general-matrix entries. + int64_t nnz_general = 0; + auto count_entry = [&](int64_t i, int64_t j) { + if (!in_triangle(i, j)) return; + nnz_general += (i == j) ? 1 : 2; + }; + if constexpr (is_coo) { + for (int64_t p = 0; p < A.nnz; ++p) + count_entry((int64_t) A.rows[p], (int64_t) A.cols[p]); + } else if constexpr (is_csr) { + for (int64_t i = 0; i < A.n_rows; ++i) + for (int64_t p = A.rowptr[i]; p < A.rowptr[i+1]; ++p) + count_entry(i, (int64_t) A.colidxs[p]); + } else { + for (int64_t j = 0; j < A.n_cols; ++j) + for (int64_t p = A.colptr[j]; p < A.colptr[j+1]; ++p) + count_entry((int64_t) A.rowidxs[p], j); + } + + // Pass 2: fill. + COOMatrix G(A.n_rows, A.n_cols); + if (nnz_general == 0) return G; + reserve_coo(nnz_general, G); + int64_t q = 0; + auto emit_entry = [&](int64_t i, int64_t j, T v) { + if (!in_triangle(i, j)) return; + G.rows[q] = (sint_t) i; G.cols[q] = (sint_t) j; G.vals[q] = v; ++q; + if (i != j) { + G.rows[q] = (sint_t) j; G.cols[q] = (sint_t) i; G.vals[q] = v; ++q; + } + }; + if constexpr (is_coo) { + for (int64_t p = 0; p < A.nnz; ++p) + emit_entry((int64_t) A.rows[p], (int64_t) A.cols[p], A.vals[p]); + } else if constexpr (is_csr) { + for (int64_t i = 0; i < A.n_rows; ++i) + for (int64_t p = A.rowptr[i]; p < A.rowptr[i+1]; ++p) + emit_entry(i, (int64_t) A.colidxs[p], A.vals[p]); + } else { + for (int64_t j = 0; j < A.n_cols; ++j) + for (int64_t p = A.colptr[j]; p < A.colptr[j+1]; ++p) + emit_entry((int64_t) A.rowidxs[p], j, A.vals[p]); + } + return G; +} + +} // end namespace RandBLAS::sparse_data + + +namespace RandBLAS { + using RandBLAS::sparse_data::Symmetric; + using RandBLAS::sparse_data::as_symmetric; +} 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/RandBLAS/util.hh b/RandBLAS/util.hh index 62733f9e..823f13dd 100644 --- a/RandBLAS/util.hh +++ b/RandBLAS/util.hh @@ -165,7 +165,10 @@ 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. +/// (The legacy sketch_symmetric overloads, which take a trailing +/// \math{\ttt{sym_check_tol}}, call this function before sketching. The +/// blas::Uplo overloads read only the named triangle, so they never invoke +/// it; it also 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/examples/CMakeLists.txt b/examples/CMakeLists.txt index f4d32fd2..965b1880 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -140,6 +140,16 @@ target_link_libraries( saso_sampling_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 +) + foreach(example_target IN ITEMS tls_dense_skop tls_sparse_skop @@ -149,6 +159,7 @@ foreach(example_target IN ITEMS spmm_performance sketch_general_performance saso_sampling_performance + spsymm_performance ) randblas_stage_runtime_dlls(${example_target}) endforeach() diff --git a/examples/simple-kernel-benchmarks/spmm_performance.cc b/examples/simple-kernel-benchmarks/spmm_performance.cc index ff3ab37e..1f73d228 100644 --- a/examples/simple-kernel-benchmarks/spmm_performance.cc +++ b/examples/simple-kernel-benchmarks/spmm_performance.cc @@ -103,13 +103,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/simple-kernel-benchmarks/spsymm_performance.cc b/examples/simple-kernel-benchmarks/spsymm_performance.cc new file mode 100644 index 00000000..6063879c --- /dev/null +++ b/examples/simple-kernel-benchmarks/spsymm_performance.cc @@ -0,0 +1,367 @@ +// 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. +// + +// ============================================================================ +// SPSYMM / SKETCH_SYMMETRIC PERFORMANCE BENCHMARK +// ============================================================================ +// +// This benchmark answers two performance questions about symmetric-aware +// multiplication kernels: +// +// 1. Sparse: how much does the one-triangle RandBLAS::spsymm path (MKL +// SPARSE_MATRIX_TYPE_SYMMETRIC fast path where available) gain over the +// both-triangles workaround (materialize both triangles into a general +// sparse matrix and call RandBLAS::spmm)? +// +// 2. Dense: how much does the SYMM-backed RandBLAS::sketch_symmetric gain +// over the equivalent sketch_general call (a GEMM that reads the full +// matrix and exploits no symmetry)? +// +// Every timed method is verified against a dense blas::symm reference before +// its timing rows are printed; a FAIL note marks any row whose output +// diverged. +// +// NOTATION: +// A_symm - symmetric matrix of order n_A (dense or sparse, one triangle) +// B - dense matrix (n_A x d for side=Left) +// C - 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 [--help] [--threads T] # default sweep +// ./spsymm_performance [--threads T] n_A d density [trials] # single config +// +// Defaults: trials=10, density=0.05, ambient OpenMP thread count. +// +// ============================================================================ + +#include +#include "RandBLAS/testing/benchmarking.hh" + +#include +#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; + +namespace bench = RandBLAS::testing; + + +// ============================================================================ +// 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 rather than a RandBLAS sampler: A is arbitrary fixed + // input data here, not a sketching operator. + 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) + A[i + j * n] = uni(rng); + RandBLAS::symmetrize(blas::Layout::ColMajor, blas::Uplo::Upper, n, A.data(), n); +} + +// 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 by zeroing the strict lower triangle of a +// copy of A_dense and converting with the library's dense_to_csr. +SpMat build_csr_upper_only(int64_t n, const std::vector& A_dense) { + std::vector upper(A_dense); + RandBLAS::overwrite_triangle(blas::Layout::ColMajor, blas::Uplo::Lower, n, 1, upper.data(), n); + SpMat A_sparse(n, n); + RandBLAS::sparse_data::csr::dense_to_csr(blas::Layout::ColMajor, upper.data(), T(0), A_sparse); + return A_sparse; +} + +// Build a full (both-triangles-stored) CSR for the general-spmm workaround. +SpMat build_csr_both_triangles(int64_t n, const std::vector& A_dense) { + std::vector full(A_dense); + SpMat A_sparse(n, n); + RandBLAS::sparse_data::csr::dense_to_csr(blas::Layout::ColMajor, full.data(), T(0), A_sparse); + return A_sparse; +} + + +// ============================================================================ +// Timing helper: (min, median) microseconds 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(static_cast(duration_cast(end - start).count())); + } + std::sort(times.begin(), times.end()); + return {times[0], times[num_trials / 2]}; +} + +// Max relative elementwise deviation between two equally-shaped buffers. +double max_rel_error(const std::vector& actual, const std::vector& expect) { + double err = 0.0; + for (size_t i = 0; i < expect.size(); ++i) { + double scale = std::max(1.0, std::abs(static_cast(expect[i]))); + err = std::max(err, std::abs(static_cast(actual[i] - expect[i])) / scale); + } + return err; +} + +void print_row(const std::string& name, int64_t min_us, int64_t med_us, int64_t baseline, bool pass) { + double ratio = (baseline > 0) ? static_cast(min_us) / static_cast(baseline) : 1.0; + std::cout << " " << std::setw(40) << std::left << name + << std::setw(10) << std::right << min_us + << std::setw(10) << med_us + << std::setw(9) << std::fixed << std::setprecision(2) << ratio << "x" + << " " << (pass ? "PASS" : "FAIL") << "\n"; +} + + +// ============================================================================ +// run_config: one (n_A, d, density) point. Every method writes its own copy +// of C, is checked against the dense blas::symm reference, and is then timed. +// ============================================================================ +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); + // Verification tolerance: the sparse and dense paths accumulate in + // different orders; a loose 1e-10 relative bound catches routing bugs + // without tripping on roundoff. + const double check_tol = 1e-10; + + std::cout << "--- A is " << n_A << "x" << n_A + << " symmetric, B is " << n_A << "x" << d + << ", density=" << std::setprecision(4) << density + << " (min + median 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); + + std::cout << " upper-triangle nnz = " << A_csr_upper.rowptr[n_A] + << " (full = " << A_csr_full.rowptr[n_A] << ")\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); + } + + // Dense blas::symm reference output (also the timed dense baseline). + std::vector C_ref(static_cast(n_A) * d, T(0)); + blas::symm(layout, Side::Left, Uplo::Upper, n_A, d, + alpha, A_full.data(), n_A, B.data(), n_A, beta, C_ref.data(), n_A); + + std::cout << "\n SPARSE (side=Left, C = A*B):\n"; + std::cout << " " << std::setw(40) << std::left << "kernel" + << std::setw(10) << std::right << "min(us)" + << std::setw(10) << "med(us)" + << std::setw(10) << "ratio" << " check\n"; + + std::vector C(static_cast(n_A) * d, T(0)); + + // 1. One-triangle storage through RandBLAS::spsymm. + RandBLAS::spsymm(layout, Uplo::Upper, d, + alpha, A_csr_upper, B.data(), n_A, beta, C.data(), n_A); + bool ok1 = max_rel_error(C, C_ref) <= check_tol; + auto [t1_min, t1_med] = run_trials([&] { + RandBLAS::spsymm(layout, Uplo::Upper, d, + alpha, A_csr_upper, B.data(), n_A, + beta, C.data(), n_A); + }, num_trials); + + // 2. Both-triangles workaround through general RandBLAS::spmm. + RandBLAS::spmm(layout, Op::NoTrans, Op::NoTrans, n_A, d, n_A, + alpha, A_csr_full, B.data(), n_A, beta, C.data(), n_A); + bool ok2 = max_rel_error(C, C_ref) <= check_tol; + 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, C.data(), n_A); + }, num_trials); + + // 3. Dense blas::symm on the fully populated dense A: does not 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, C.data(), n_A); + }, num_trials); + + print_row("RandBLAS::spsymm (one triangle)", t1_min, t1_med, t1_min, ok1); + print_row("RandBLAS::spmm (both triangles)", t2_min, t2_med, t1_min, ok2); + print_row("blas::symm (dense reference)", t3_min, t3_med, t1_min, true); + + // Dense-symmetric sketching comparison: SYMM-backed sketch_symmetric vs + // the GEMM-forwarding sketch_general equivalent. + std::cout << "\n DENSE (sketch_symmetric vs sketch_general on dense-symm A):\n"; + std::cout << " " << std::setw(40) << std::left << "kernel" + << std::setw(10) << std::right << "min(us)" + << std::setw(10) << "med(us)" + << std::setw(10) << "ratio" << " check\n"; + + RandBLAS::DenseDist DS(n_A, d, RandBLAS::ScalarDist::Uniform); + RandBLAS::DenseSkOp S(DS, static_cast(seed + 3)); + RandBLAS::fill_dense(S); + + // Reference for the sketching comparison: the GEMM-forwarding + // sketch_general result (reads the full symmetric matrix, so it is the + // trusted baseline; the check below is SYMM-path vs GEMM-path agreement). + std::vector C_sk_ref(static_cast(n_A) * d, T(0)); + RandBLAS::sketch_general(layout, Op::NoTrans, Op::NoTrans, n_A, d, n_A, + alpha, A_full.data(), n_A, S, 0, 0, beta, C_sk_ref.data(), n_A); + + // 1. SYMM-backed sketch_symmetric. + RandBLAS::sketch_symmetric(layout, Uplo::Upper, n_A, d, + alpha, A_full.data(), n_A, S, 0, 0, beta, C.data(), n_A); + bool ok3 = max_rel_error(C, C_sk_ref) <= check_tol; + 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, C.data(), n_A); + }, num_trials); + + // 2. Equivalent call via sketch_general (GEMM, no symmetry exploited; + // this is the reference, so its check is definitionally PASS). + 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, C.data(), n_A); + }, num_trials); + + print_row("sketch_symmetric (SYMM path)", s1_min, s1_med, s1_min, ok3); + print_row("sketch_general (GEMM path)", s2_min, s2_med, s1_min, true); + + std::cout << "\n"; +} + + +// ============================================================================ +// main: default 3-point sweep, or single config from positional arguments. +// ============================================================================ +void print_usage(const char* prog) { + std::cout << "Usage:\n" + << " " << prog << " [--help] [--threads T] default sweep\n" + << " " << prog << " [--threads T] n_A d density [trials] single configuration\n" + << "n_A, d, trials are positive integers; density is in (0, 1].\n"; +} + +int main(int argc, char** argv) { + bench::OpenMPSettingsGuard omp_guard; + + std::vector args(argv + 1, argv + argc); + int requested_threads = 0; + for (size_t i = 0; i < args.size(); ) { + if (args[i] == "--help") { + print_usage(argv[0]); + return 0; + } else if (args[i] == "--threads" && i + 1 < args.size()) { + requested_threads = std::atoi(args[i + 1].c_str()); + if (requested_threads <= 0) { + std::cerr << "Invalid configuration. Expected a positive --threads value.\n"; + return 1; + } + args.erase(args.begin() + i, args.begin() + i + 2); + } else { + ++i; + } + } + if (requested_threads > 0) + bench::set_threads(requested_threads); + std::cout << "OpenMP threads: " << bench::current_threads(); + if (requested_threads > 0) + std::cout << " (requested " << requested_threads + << ", effective " << bench::effective_threads(requested_threads) << ")"; + std::cout << "\n\n"; + + if (args.size() >= 3) { + int64_t n_A = std::atoll(args[0].c_str()); + int64_t d = std::atoll(args[1].c_str()); + double density = std::atof(args[2].c_str()); + int num_trials = (args.size() >= 4) ? std::atoi(args[3].c_str()) : 10; + if (n_A <= 0 || d <= 0 || density <= 0.0 || density > 1.0 || num_trials <= 0) { + std::cerr << "Invalid configuration. Expected positive n_A, d, trials and density in (0, 1].\n"; + print_usage(argv[0]); + return 1; + } + run_config(n_A, d, density, num_trials); + } else if (args.empty()) { + 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); + } + } else { + std::cerr << "Invalid configuration. Expected zero or 3-4 positional arguments.\n"; + print_usage(argv[0]); + return 1; + } + return 0; +} diff --git a/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc b/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc index 4c82f68e..27955aec 100644 --- a/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc +++ b/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc @@ -339,9 +339,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]; diff --git a/rtd/source/FAQ.rst b/rtd/source/FAQ.rst index ea9a0e1c..b7a0d02a 100644 --- a/rtd/source/FAQ.rst +++ b/rtd/source/FAQ.rst @@ -7,10 +7,15 @@ How do I do this and that? -------------------------- How do I sketch a const symmetric matrix that's only stored in an upper or lower triangle? - You can only do this with dense sketching operators. - You'll have to prepare the plain buffer representation yourself with - :cpp:any:`RandBLAS::fill_dense_unpacked`, - and then you'll have to use that buffer in your own SYMM function. + If the matrix is dense, call sketch_symmetric and pass the ``blas::Uplo`` value naming the + stored triangle. This works with dense and sparse sketching operators alike, and the opposite + triangle is never read. + + There is no symmetry-exploiting sketch for a sparse matrix. Either sketch it as a general + sparse matrix with sketch_sparse (which reads exactly the entries you populate, so this only + gives a correct answer if you first materialize both triangles), or, if what you actually want + is a deterministic product against a sparse symmetric matrix rather than a sketch of it, wrap + it with ``as_symmetric`` and pass the wrapper to ``spmm``. How do I sketch a submatrix of a sparse matrix? You can only do this if the sparse matrix in COO format. @@ -106,14 +111,15 @@ No support for DenseSkOps with Rademachers: No support for negative values of "incx" or "incy" in sketch_vector. The BLAS function GEMV supports negative strides between input and output vector elements. It would be easy to extend sketch_vector to support this if we had a proper - SPMV implementation that supported negative increments. If someone wants to volunteer + SPMV implementation that supported negative increments. If someone wants to volunteer 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. +No zero-copy path for a layout-mismatched DenseSkOp in sketch_symmetric. + ``blas::symm`` has no on-the-fly transpose flag for the dense operand, so if a ``DenseSkOp``'s + storage layout differs from the caller's ``layout`` argument, ``sketch_symmetric`` pays an + ``O(d * n)`` allocation and transpose-copy of the operand to keep the SYMM speedup over a + ``blas::gemm`` fallback. Matching the operand's layout to the call's ``layout`` avoids the copy. Language interoperability @@ -168,12 +174,14 @@ We have a consistent naming convention for functions that involve sketching oper Functions that implement the overload-free conventions * [L/R]skge[X] for sketching a general matrix from the left (L) or right (R) with a matrix whose structure is indicated by [X]. C++ code should prefer overloaded sketch_general - * [L/R]sksp3 for sketching a sparse matrix from the left (L) (L) or right (R) with a DenseSkOp. + * [L/R]sksp[X] for sketching a sparse matrix from the left (L) or right (R), both for a DenseSkOp (3) and a SparseSkOp (s). C++ code should prefer overloaded sketch_sparse, unless operating on a submatrix of a COO-format sparse data matrix is needed. + The SparseSkOp form requires Intel MKL and, unlike the DenseSkOp form, does not support a submatrix of the sparse matrix being sketched. + * [L/R]sksy[X] for sketching a matrix with *explicit symmetry*, both for a DenseSkOp (3) and a SparseSkOp (s). + C++ code should prefer overloaded sketch_symmetric. Functions that are missing implementations of this convention * [L/R]skve[X] for sketching vectors. This functionality is availabile in C++ with sketch_vector - * [L/R]sksy[X] for sketching a matrix with *explicit symmetry*. This functionality is availabile in C++ with sketch_symmetric. Some discussion diff --git a/rtd/source/api_reference/sketch_dense.rst b/rtd/source/api_reference/sketch_dense.rst index 7dd0a571..83db2602 100644 --- a/rtd/source/api_reference/sketch_dense.rst +++ b/rtd/source/api_reference/sketch_dense.rst @@ -67,18 +67,22 @@ Analogs to GEMM Analogs to SYMM --------------- +These overloads accept a ``blas::Uplo`` parameter naming the triangle of +:math:`\mtxA` that is stored; the opposite triangle is implied by symmetry +and is not read. Both dense and sparse sketching operators are supported. + .. 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,6 +90,28 @@ Analogs to SYMM :animate: fade-in-slide-down :color: light + .. 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, 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 + + +.. dropdown:: Legacy variants (``sym_check_tol``) + :animate: fade-in-slide-down + :color: light + + These reproduce the pre-``Uplo`` API exactly: :math:`\mtxA` must be stored + with both triangles populated (both are read), and a runtime symmetry + check runs first (pass a negative tolerance to skip it). They forward to + ``sketch_general``. Prefer the ``blas::Uplo`` overloads for new code. + + .. 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) + :project: RandBLAS + + .. 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) + :project: RandBLAS + .. 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) :project: RandBLAS diff --git a/rtd/source/api_reference/sketch_sparse.rst b/rtd/source/api_reference/sketch_sparse.rst index 98a14df9..32ccb7b0 100644 --- a/rtd/source/api_reference/sketch_sparse.rst +++ b/rtd/source/api_reference/sketch_sparse.rst @@ -70,18 +70,37 @@ Operations with sparse matrices Sketching ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -.. dropdown:: :math:`\mtxB = \alpha \cdot \op(\submat(\mtxS))\cdot \op(\mtxA) + \beta \cdot \mtxB` +``sketch_sparse`` accepts either a dense or a sparse sketching operator. The dense-operator +overloads dispatch to ``left_spmm`` / ``right_spmm`` (a materialized dense operand); the +sparse-operator overloads dispatch to ``spgemm`` and so require Intel MKL, and (unlike the +dense-operator overloads) :math:`\mtxA` must be passed in full, not as a submatrix. + +.. dropdown:: :math:`\mtxB = \alpha \cdot \op(\submat(\mtxS))\cdot \op(\mtxA) + \beta \cdot \mtxB,` :math:`\mtxS` dense + :animate: fade-in-slide-down + :color: light + + .. doxygenfunction:: RandBLAS::sketch_sparse(blas::Layout layout, blas::Op opS, blas::Op opA, int64_t d, int64_t n, int64_t m, T alpha, const DenseSkOp &S, int64_t S_ro, int64_t S_co, const SpMat &A, T beta, T *B, int64_t ldb) + :project: RandBLAS + +.. dropdown:: :math:`\mtxB = \alpha \cdot \op(\mtxA)\cdot \op(\submat(\mtxS)) + \beta \cdot \mtxB,` :math:`\mtxS` dense :animate: fade-in-slide-down :color: light - .. doxygenfunction:: RandBLAS::sketch_sparse(blas::Layout layout, blas::Op opS, blas::Op opA, int64_t d, int64_t n, int64_t m, T alpha, const DenseSkOp &S, int64_t S_ro, int64_t S_co, const SpMat &A, T beta, T *B, int64_t ldb) + .. doxygenfunction:: RandBLAS::sketch_sparse(blas::Layout layout, blas::Op opA, blas::Op opS, int64_t m, int64_t d, int64_t n, T alpha, const SpMat &A, const DenseSkOp &S, int64_t S_ro, int64_t S_co, T beta, T *B, int64_t ldb) :project: RandBLAS -.. dropdown:: :math:`\mtxB = \alpha \cdot \op(\mtxA)\cdot \op(\submat(\mtxS)) + \beta \cdot \mtxB` +.. dropdown:: :math:`\mtxB = \alpha \cdot \op(\submat(\mtxS))\cdot \op(\mtxA) + \beta \cdot \mtxB,` :math:`\mtxS` sparse :animate: fade-in-slide-down :color: light - .. doxygenfunction:: RandBLAS::sketch_sparse(blas::Layout layout, blas::Op opA, blas::Op opS, int64_t m, int64_t d, int64_t n, T alpha, const SpMat &A, const DenseSkOp &S, int64_t S_ro, int64_t S_co, T beta, T *B, int64_t ldb) + .. doxygenfunction:: RandBLAS::sketch_sparse(blas::Layout layout, blas::Op opS, blas::Op opA, int64_t d, int64_t n, int64_t m, T alpha, const SparseSkOp &S, int64_t ro_s, int64_t co_s, const SpMat &A, T beta, T *B, int64_t ldb) + :project: RandBLAS + +.. dropdown:: :math:`\mtxB = \alpha \cdot \op(\mtxA)\cdot \op(\submat(\mtxS)) + \beta \cdot \mtxB,` :math:`\mtxS` sparse + :animate: fade-in-slide-down + :color: light + + .. doxygenfunction:: RandBLAS::sketch_sparse(blas::Layout layout, blas::Op opA, blas::Op opS, int64_t m, int64_t d, int64_t n, T alpha, const SpMat &A, const SparseSkOp &S, int64_t ro_s, int64_t co_s, T beta, T *B, int64_t ldb) :project: RandBLAS @@ -114,6 +133,38 @@ 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:`\mtxC = \alpha \cdot \mtxA \cdot \mtxB + \beta \cdot \mtxC,` with sparse symmetric :math:`\mtxA` + :animate: fade-in-slide-down + :color: light + + .. doxygenfunction:: RandBLAS::spsymm(blas::Layout layout, blas::Uplo uplo, int64_t n, T alpha, const SpMat &A, const T *B, int64_t ldb, T beta, T *C, int64_t ldc) + :project: RandBLAS + + .. doxygenfunction:: RandBLAS::spmm(blas::Layout layout, int64_t n, T alpha, const Symmetric &A_sym, const T *B, int64_t ldb, T beta, T *C, int64_t ldc) + :project: RandBLAS + + .. doxygenstruct:: RandBLAS::sparse_data::Symmetric + :project: RandBLAS + :members: + + .. doxygenfunction:: RandBLAS::sparse_data::as_symmetric(const SpMat &A, blas::Uplo uplo) + :project: RandBLAS + + Only the triangle of :math:`\mtxA` named by ``uplo`` is read; the opposite + triangle is implied by symmetry. :math:`\mtxA` must be square, and its + order is taken from the matrix itself. + +.. dropdown:: :math:`\mtxC = \alpha \cdot \mtxA \cdot \mtxB + \beta \cdot \mtxC,` with sparse symmetric :math:`\mtxA` and sparse :math:`\mtxB` + :animate: fade-in-slide-down + :color: light + + .. doxygenfunction:: RandBLAS::sparse_data::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 *C, int64_t ldc) + :project: RandBLAS + + Only the triangle of :math:`\mtxA` named by ``uplo`` is read; the opposite + triangle is implied by symmetry. Requires Intel MKL for the accelerated + path; without it, :math:`\mtxB` is densified into a temporary buffer. + .. dropdown:: :math:`\mtxB = \alpha \cdot \op(\mtxA)^{-1} \cdot \mtxB,` with sparse triangular :math:`\mtxA` :animate: fade-in-slide-down :color: light diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 1e9a4baf..77a0e68c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -39,12 +39,14 @@ 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 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/datastructures/test_symmetric_wrapper.cc b/test/datastructures/test_symmetric_wrapper.cc new file mode 100644 index 00000000..cb3311b2 --- /dev/null +++ b/test/datastructures/test_symmetric_wrapper.cc @@ -0,0 +1,96 @@ +// 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, constructs_from_coo) { + 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); + EXPECT_EQ(wrapped.n_rows, 4); + EXPECT_EQ(wrapped.n_cols, 4); + static_assert(std::is_same_v, + "Symmetric>::scalar_t must be double"); +} + +TEST_F(TestSymmetricWrapper, constructs_from_csr) { + 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, constructs_from_csc) { + 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, exports_at_randblas_scope) { + // 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, non_square_rejected_at_construction) { + COOMatrix A(3, 4); // not square + EXPECT_THROW({ + auto wrapped = as_symmetric(A, blas::Uplo::Upper); + (void) wrapped; + }, RandBLAS::Error); +} diff --git a/test/linops/test_sketch_sparse.cc b/test/linops/test_sketch_sparse.cc index f511c738..ea5702fc 100644 --- a/test/linops/test_sketch_sparse.cc +++ b/test/linops/test_sketch_sparse.cc @@ -639,3 +639,237 @@ TEST_F(TestRSKSP3, submatrix_s_single) Layout::ColMajor ); } + + +//////////////////////////////////////////////////////////////////////// +// +// +// LSKSPS / RSKSPS: sketching a sparse matrix with a SparseSkOp +// +// +//////////////////////////////////////////////////////////////////////// + +// These dispatch through spgemm, which is a compile-time error without +// Intel MKL, so the whole block is gated the same way test_spgemm.cc is. +#if defined(RandBLAS_HAS_MKL) + +using RandBLAS::SparseSkOp; +using RandBLAS::SparseDist; + +// Adapted from test_left_submat_sketch_of_eye above, with the DenseSkOp +// replaced by a SparseSkOp and lsksp3 replaced by lsksps. +template > +void test_left_submat_sketch_of_eye_sparse( + T alpha, SparseSkOp &S0, int64_t d1, int64_t m1, int64_t S_ro, int64_t S_co, Layout layout, T beta = 0.0 +) { + auto [d0, m0] = dimensions(S0); + randblas_require(d0 >= d1); + randblas_require(m0 >= m1); + bool is_colmajor = layout == Layout::ColMajor; + int64_t ldb = (is_colmajor) ? d1 : m1; + + auto I = eye(m1); + auto B = std::get<0>(random_matrix(d1, m1, RNGState(42))); + std::vector B_backup(B); + + lsksps( + layout, Op::NoTrans, Op::NoTrans, d1, m1, m1, + alpha, S0, S_ro, S_co, I, 0, 0, beta, B.data(), ldb + ); + + T *expect = new T[d0 * m0]; + to_explicit_buffer(S0, expect, layout); + int64_t ld_expect = (is_colmajor) ? d0 : m0; + auto [row_stride_s, col_stride_s] = layout_to_strides(layout, ld_expect); + auto [row_stride_b, col_stride_b] = layout_to_strides(layout, ldb); + int64_t offset = row_stride_s * S_ro + col_stride_s * S_co; + #define MAT_ES(_i, _j) expect[offset + (_i)*row_stride_s + (_j)*col_stride_s] + #define MAT_BS(_i, _j) B_backup[ (_i)*row_stride_b + (_j)*col_stride_b] + for (int i = 0; i < d1; ++i) { + for (int j = 0; j < m1; ++j) { + MAT_ES(i,j) = alpha * MAT_ES(i,j) + beta * MAT_BS(i, j); + } + } + + auto msg = RandBLAS::testing::matrices_approx_equal( + layout, Op::NoTrans, + d1, m1, + B.data(), ldb, + &expect[offset], ld_expect, + __RANDBLAS_PRETTY_FUNCTION__, __FILE__, __LINE__ + ); + if (msg.size() > 0) { + FAIL() << msg; + } + delete [] expect; +} + +// B = S^T * eye, mirroring test_left_transposed_sketch_of_eye above. +template > +void test_left_transposed_sketch_of_eye_sparse(SparseSkOp &S, Layout layout) { + auto [m, d] = dimensions(S); + auto I = eye(m); + std::vector B(d * m, 0.0); + bool is_colmajor = (Layout::ColMajor == layout); + int64_t ldb = (is_colmajor) ? d : m; + int64_t lds = (is_colmajor) ? m : d; + + lsksps( + layout, Op::Trans, Op::NoTrans, d, m, m, + (T) 1.0, S, 0, 0, I, 0, 0, (T) 0.0, B.data(), ldb + ); + + std::vector S_dense(m * d, 0.0); + to_explicit_buffer(S, S_dense.data(), layout); + auto msg = RandBLAS::testing::matrices_approx_equal( + layout, Op::Trans, d, m, + B.data(), ldb, S_dense.data(), lds, + __RANDBLAS_PRETTY_FUNCTION__, __FILE__, __LINE__ + ); + if (msg.size() > 0) { + FAIL() << msg; + } +} + +// B = eye * S, via the public sketch_sparse(SparseSkOp) overload rather than +// the low-level rsksps kernel, so the public entry point gets covered too. +template > +void test_right_sketch_of_eye_sparse(SparseSkOp &S, Layout layout) { + auto [n, d] = dimensions(S); + auto I = eye(n); + std::vector B(n * d, 0.0); + bool is_colmajor = (Layout::ColMajor == layout); + int64_t ldb = (is_colmajor) ? n : d; + int64_t lds = (is_colmajor) ? n : d; + + sketch_sparse( + layout, Op::NoTrans, Op::NoTrans, n, d, n, + (T) 1.0, I, S, 0, 0, (T) 0.0, B.data(), ldb + ); + + std::vector S_dense(n * d, 0.0); + to_explicit_buffer(S, S_dense.data(), layout); + auto msg = RandBLAS::testing::matrices_approx_equal( + layout, Op::NoTrans, n, d, + B.data(), ldb, S_dense.data(), lds, + __RANDBLAS_PRETTY_FUNCTION__, __FILE__, __LINE__ + ); + if (msg.size() > 0) { + FAIL() << msg; + } +} + + +class TestLSKSPS : public ::testing::Test +{ + protected: + + virtual void SetUp(){}; + + virtual void TearDown(){}; + + template + static void sketch_eye(uint32_t seed, int64_t m, int64_t d, Layout layout) { + SparseDist D(d, m, std::min(int64_t(3), std::min(d, m))); + SparseSkOp S0(D, seed); + fill_sparse(S0); + test_left_submat_sketch_of_eye_sparse(1.0, S0, d, m, 0, 0, layout, 0.0); + } + + template + static void transpose_S(uint32_t seed, int64_t m, int64_t d, Layout layout) { + SparseDist Dt(m, d, std::min(int64_t(3), std::min(m, d))); + SparseSkOp S0(Dt, seed); + fill_sparse(S0); + test_left_transposed_sketch_of_eye_sparse(S0, layout); + } + + template + static void submatrix_S( + uint32_t seed, int64_t d, int64_t m, int64_t d0, int64_t m0, + int64_t S_ro, int64_t S_co, Layout layout + ) { + randblas_require(d0 > d); + randblas_require(m0 > m); + SparseDist D(d0, m0, std::min(int64_t(3), std::min(d0, m0))); + SparseSkOp S0(D, seed); + fill_sparse(S0); + test_left_submat_sketch_of_eye_sparse(1.0, S0, d, m, S_ro, S_co, layout, 0.0); + } + + template + static void nontrivial_alpha_beta(uint32_t seed, int64_t m, int64_t d, Layout layout) { + SparseDist D(d, m, std::min(int64_t(3), std::min(d, m))); + SparseSkOp S0(D, seed); + fill_sparse(S0); + test_left_submat_sketch_of_eye_sparse(1.7, S0, d, m, 0, 0, layout, -0.4); + } +}; + +TEST_F(TestLSKSPS, sketch_eye_double_colmajor) { + sketch_eye(0, 12, 5, Layout::ColMajor); +} + +TEST_F(TestLSKSPS, sketch_eye_double_rowmajor) { + sketch_eye(0, 12, 5, Layout::RowMajor); +} + +TEST_F(TestLSKSPS, sketch_eye_single) { + sketch_eye(0, 12, 5, Layout::ColMajor); +} + +TEST_F(TestLSKSPS, transpose_double_colmajor) { + transpose_S(0, 12, 5, Layout::ColMajor); +} + +TEST_F(TestLSKSPS, transpose_double_rowmajor) { + transpose_S(0, 12, 5, Layout::RowMajor); +} + +TEST_F(TestLSKSPS, submatrix_s_double_colmajor) { + submatrix_S(0, 5, 8, 9, 12, 2, 1, Layout::ColMajor); +} + +TEST_F(TestLSKSPS, submatrix_s_double_rowmajor) { + submatrix_S(0, 5, 8, 9, 12, 2, 1, Layout::RowMajor); +} + +TEST_F(TestLSKSPS, nontrivial_alpha_beta_colmajor) { + nontrivial_alpha_beta(0, 12, 5, Layout::ColMajor); +} + +TEST_F(TestLSKSPS, nontrivial_alpha_beta_rowmajor) { + nontrivial_alpha_beta(0, 12, 5, Layout::RowMajor); +} + + +class TestRSKSPS : public ::testing::Test +{ + protected: + + virtual void SetUp(){}; + + virtual void TearDown(){}; + + template + static void sketch_eye(uint32_t seed, int64_t n, int64_t d, Layout layout) { + SparseDist D(n, d, std::min(int64_t(3), std::min(n, d))); + SparseSkOp S(D, seed); + fill_sparse(S); + test_right_sketch_of_eye_sparse(S, layout); + } +}; + +TEST_F(TestRSKSPS, sketch_eye_double_colmajor) { + sketch_eye(0, 12, 5, Layout::ColMajor); +} + +TEST_F(TestRSKSPS, sketch_eye_double_rowmajor) { + sketch_eye(0, 12, 5, Layout::RowMajor); +} + +TEST_F(TestRSKSPS, sketch_eye_single) { + sketch_eye(0, 12, 5, Layout::ColMajor); +} + +#endif // RandBLAS_HAS_MKL diff --git a/test/linops/test_sketch_symmetric.cc b/test/linops/test_sketch_symmetric.cc index 4156c3a5..356bd151 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; @@ -60,14 +63,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 +92,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 +109,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 +125,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 +150,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,35 +166,75 @@ class TestSketchSymmetric : public ::testing::Test { return; } + // Note on symmetry checking: the blas::Uplo overloads read only the + // named triangle, so they perform no runtime symmetry check; the legacy + // sym_check_tol overloads retain the check (covered under MARK: LEGACY + // OVERLOADS below). Error paths for the Uplo overloads (leading dims, + // submatrix window bounds) are covered under MARK: ERROR PATHS. + + // ============================================================================= + // 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). With materialize=false, S is handed to + // sketch_symmetric unmaterialized (nnz < 0), exercising the + // submatrix_as_coo window-sampling path; the reference is always built + // from a materialized twin (same dist and seed give identical samples). 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; + 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, + bool materialize = true + ) { + 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); + if (materialize) + RandBLAS::fill_sparse(S); + SparseSkOp S_ref(DS, seed_skop); + RandBLAS::fill_sparse(S_ref); + + // Densify the reference twin into a buffer matching the requested + // layout, for the SYMM reference. lds is the major-axis leading dim + // (tight, so coo_to_dense's layout overload applies directly). + 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_ref); + RandBLAS::sparse_data::coo::coo_to_dense(Scoo, layout, S_dense.data()); + + 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 + // column-driven sparse kernel accumulates in a different 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, + __RANDBLAS_PRETTY_FUNCTION__, __FILE__, __LINE__, atol, rtol + ); + if (!msg.empty()) FAIL() << msg; } }; @@ -391,9 +437,155 @@ 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(); +// 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 symmetric-read resolution). + 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); +} + + +TEST_F(TestSketchSymmetric, sparse_skop_unmaterialized) { + // S is handed to sketch_symmetric with nnz < 0: the wrapper samples only + // the requested window via submatrix_as_coo instead of materializing S. + test_sparse_skop(Layout::ColMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, 0.0, blas::Side::Left, 2, Uplo::Upper, /*materialize=*/false); + test_sparse_skop(Layout::RowMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, 0.0, blas::Side::Left, 2, Uplo::Lower, /*materialize=*/false); + test_sparse_skop(Layout::ColMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, -1.0, blas::Side::Right, 2, Uplo::Upper, /*materialize=*/false); +} + + +// MARK: ERROR PATHS + +TEST_F(TestSketchSymmetric, sparse_skop_bad_arguments_throw) { + int64_t d = 3, n = 10; + SparseDist DS(d, n, 2, Axis::Short); + SparseSkOp S(DS, 11); + RandBLAS::fill_sparse(S); + std::vector A(n * n, 1.0); + std::vector B(static_cast(d) * n, 0.0); + + // ldb below its ColMajor lower bound (B is d-by-n, so ldb >= d). + ASSERT_THROW( + RandBLAS::sketch_symmetric(Layout::ColMajor, Uplo::Upper, d, n, + 1.0, S, 0, 0, A.data(), n, 0.0, B.data(), d - 1), + RandBLAS::Error); + // lda below its lower bound (A is n-by-n). + ASSERT_THROW( + RandBLAS::sketch_symmetric(Layout::ColMajor, Uplo::Upper, d, n, + 1.0, S, 0, 0, A.data(), n - 1, 0.0, B.data(), d), + RandBLAS::Error); + // Submatrix window exceeding the operator: ro_s + d > S.n_rows. The + // dense-SkOp branch has always thrown here; the sparse branch must too + // (a silent filter would return a sketch with missing rows). + ASSERT_THROW( + RandBLAS::sketch_symmetric(Layout::ColMajor, Uplo::Upper, d, n, + 1.0, S, 1, 0, A.data(), n, 0.0, B.data(), d), + RandBLAS::Error); +} + + +// MARK: LEGACY OVERLOADS + +// Restores the pre-Uplo API test: the legacy sym_check_tol overloads must +// reject an asymmetric matrix at runtime, exactly as before. +TEST_F(TestSketchSymmetric, legacy_symmetry_check_fails_for_asymmetric_matrix) { + int64_t n = 3, d = 2; + // Column-major 3x3, symmetric except A(0,1)=5 vs A(1,0)=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, 1.0, S, 0, 0, A.data(), n, 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; + } +} + +// The legacy overloads (both-triangles A, runtime check, sketch_general +// forwarding) must agree with the Uplo overloads on symmetric input. +TEST_F(TestSketchSymmetric, legacy_overloads_agree_with_uplo_overloads) { + int64_t n = 10, d = 3; + std::vector A(n * n, 0.0); + random_symmetric_mat(n, A.data(), n, RNGState(7)); + DenseDist D(d, n, ScalarDist::Uniform); + DenseSkOp S(D, 21); + RandBLAS::fill_dense(S); + + std::vector B_legacy(d * n, 0.5), B_uplo(d * n, 0.5); + // Left sketch, SUBMAT form, beta != 0. + RandBLAS::sketch_symmetric(Layout::ColMajor, d, n, 2.0, S, 0, 0, A.data(), n, -1.0, B_legacy.data(), d); + RandBLAS::sketch_symmetric(Layout::ColMajor, Uplo::Upper, d, n, 2.0, S, 0, 0, A.data(), n, -1.0, B_uplo.data(), d); + auto msg = RandBLAS::testing::matrices_approx_equal( + Layout::ColMajor, blas::Op::NoTrans, d, n, + B_legacy.data(), d, B_uplo.data(), d, + __RANDBLAS_PRETTY_FUNCTION__, __FILE__, __LINE__, + 100 * std::numeric_limits::epsilon(), + 10 * std::numeric_limits::epsilon() + ); + if (!msg.empty()) FAIL() << msg; + + // Right sketch, FULL form. + DenseDist DR(n, d, ScalarDist::Uniform); + DenseSkOp SR(DR, 23); + RandBLAS::fill_dense(SR); + std::vector C_legacy(n * d, 0.0), C_uplo(n * d, 0.0); + RandBLAS::sketch_symmetric(Layout::ColMajor, 1.0, A.data(), n, SR, 0.0, C_legacy.data(), n); + RandBLAS::sketch_symmetric(Layout::ColMajor, Uplo::Upper, 1.0, A.data(), n, SR, 0.0, C_uplo.data(), n); + auto msg2 = RandBLAS::testing::matrices_approx_equal( + Layout::ColMajor, blas::Op::NoTrans, n, d, + C_legacy.data(), n, C_uplo.data(), n, + __RANDBLAS_PRETTY_FUNCTION__, __FILE__, __LINE__, + 100 * std::numeric_limits::epsilon(), + 10 * std::numeric_limits::epsilon() + ); + if (!msg2.empty()) FAIL() << msg2; } diff --git a/test/linops/test_spsymm.cc b/test/linops/test_spsymm.cc new file mode 100644 index 00000000..dea00063 --- /dev/null +++ b/test/linops/test_spsymm.cc @@ -0,0 +1,417 @@ +// 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 +#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 (k=1 skips the + // diagonal), leaving only the structurally-stored side + diagonal. + auto other = (uplo == Uplo::Upper) ? Uplo::Lower : Uplo::Upper; + RandBLAS::overwrite_triangle(Layout::ColMajor, other, n, 1, A, lda); + } + + 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, + bool route_via_wrapper = false + ) { + // For side=Left: C = alpha*A*B + beta*C, A is n_A x n_A, B and C are n_A x d. + // For side=Right: C = alpha*B*A + beta*C, B and C are d x n_A, A is n_A x n_A. + int64_t m_BC, n_BC; + if (side == Side::Left) { m_BC = n_A; n_BC = d; } + else { m_BC = d; n_BC = 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 C (for beta != 0 to be non-trivial). + int64_t ldb = (layout == Layout::ColMajor) ? m_BC : n_BC; + int64_t ldc = ldb; + std::vector B(m_BC * n_BC); + DenseDist DB(m_BC, n_BC, ScalarDist::Uniform); + RandBLAS::fill_dense_unpacked(layout, DB, m_BC, n_BC, 0, 0, B.data(), RNGState(seed_B)); + + std::vector C_actual(m_BC * n_BC); + RandBLAS::fill_dense_unpacked(layout, DB, m_BC, n_BC, 0, 0, C_actual.data(), RNGState(seed_B + 7)); + std::vector C_expect = C_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_BC, n_BC, + alpha, A_full.data(), lda, B.data(), ldb, + beta, C_expect.data(), ldc); + + // 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::spmm(Symmetric) + // overload to exercise the wrapper-routing path. side=Left only for + // the wrapper path since the public overload has no side parameter + // (it is unconditionally side=Left). + if (route_via_wrapper) { + randblas_require(side == Side::Left); + auto A_sym = RandBLAS::as_symmetric(A_sparse, uplo); + RandBLAS::spmm(layout, n_BC, + alpha, A_sym, B.data(), ldb, + beta, C_actual.data(), ldc); + } else { + RandBLAS::sparse_data::spsymm(layout, side, uplo, m_BC, n_BC, + alpha, A_sparse, B.data(), ldb, + beta, C_actual.data(), ldc); + } + + // 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_BC, n_BC, + C_actual.data(), ldc, C_expect.data(), ldc, + __RANDBLAS_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); + } + } + } + + // Case D: sparse-symmetric A times sparse B -> dense C. 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_BC, n_BC; + if (side == Side::Left) { m_BC = n_A; n_BC = d; } + else { m_BC = d; n_BC = 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_BC * n_BC, T(0)); + { + // std::mt19937_64 rather than a RandBLAS sampler: B here is + // arbitrary fixed test data, not a sketching operator. + 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_BC; ++j) { + for (int64_t i = 0; i < m_BC; ++i) { + if (uni01(rng) < density_B) { + B_dense[i + j * m_BC] = static_cast(univ(rng)); + } + } + } + } + SpMatB B_sparse(m_BC, n_BC); + dense_to_sparse_format(Layout::ColMajor, B_dense.data(), T(0), B_sparse); + + int64_t ldb = (layout == Layout::ColMajor) ? m_BC : n_BC; + int64_t ldc = ldb; + // The dense reference call to blas::symm uses the requested layout. + std::vector B_dense_layout(m_BC * n_BC); + if (layout == Layout::ColMajor) { + std::copy(B_dense.begin(), B_dense.end(), B_dense_layout.begin()); + } else { + for (int64_t i = 0; i < m_BC; ++i) + for (int64_t j = 0; j < n_BC; ++j) + B_dense_layout[i * ldb + j] = B_dense[i + j * m_BC]; + } + + std::vector C_actual(m_BC * n_BC); + DenseDist DC(m_BC, n_BC, ScalarDist::Uniform); + RandBLAS::fill_dense_unpacked(layout, DC, m_BC, n_BC, 0, 0, C_actual.data(), RNGState(seed_B + 13)); + std::vector C_expect = C_actual; + + // Reference: dense blas::symm on full-storage A and dense B. + blas::symm(layout, side, uplo, m_BC, n_BC, + alpha, A_full.data(), lda, B_dense_layout.data(), ldb, + beta, C_expect.data(), ldc); + + // Under test: sparse-symm A times sparse B via Case D. + RandBLAS::sparse_data::spsymm(layout, side, uplo, m_BC, n_BC, + alpha, A_sparse, B_sparse, + beta, C_actual.data(), ldc); + + 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_BC, n_BC, + C_actual.data(), ldc, C_expect.data(), ldc, + __RANDBLAS_PRETTY_FUNCTION__, __FILE__, __LINE__, atol, rtol + ); + if (!msg.empty()) FAIL() << msg; + } +}; + + +// 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_matches_dense_symm_reference) { sweep_layout_uplo>(Side::Left, 10, 4, 1.5, -0.5); } +TEST_F(TestSpsymm, csr_right_matches_dense_symm_reference) { sweep_layout_uplo>(Side::Right, 10, 4, 1.5, -0.5); } +TEST_F(TestSpsymm, csc_left_matches_dense_symm_reference) { sweep_layout_uplo>(Side::Left, 10, 4, 1.5, -0.5); } +TEST_F(TestSpsymm, csc_right_matches_dense_symm_reference) { sweep_layout_uplo>(Side::Right, 10, 4, 1.5, -0.5); } +TEST_F(TestSpsymm, coo_left_matches_dense_symm_reference) { sweep_layout_uplo>(Side::Left, 10, 4, 1.5, -0.5); } +TEST_F(TestSpsymm, coo_right_matches_dense_symm_reference) { sweep_layout_uplo>(Side::Right, 10, 4, 1.5, -0.5); } + +// Float coverage for one representative format +TEST_F(TestSpsymm, csr_left_float_matches_dense_symm_reference) { sweep_layout_uplo>(Side::Left, 10, 4, 1.5, -0.5); } + +// beta=0 (init-from-zero) edge case +TEST_F(TestSpsymm, csr_left_beta_zero_overwrites_output) { + 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_alpha_zero_scales_by_beta) { + run_case>(Layout::ColMajor, Side::Left, Uplo::Upper, 10, 4, 0.0, 0.5, 0, 3); +} + + +// 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, sparse_times_sparse_csr_csr_left_matches_reference) { + 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, sparse_times_sparse_csc_csc_left_matches_reference) { + 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, sparse_times_sparse_coo_coo_left_matches_reference) { + 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, sparse_times_sparse_mixed_formats_match_reference) { + // All six mixed (A, B) format pairings; the three same-format pairings + // are covered by the tests above. + 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); + run_case_d, COOMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 1.0, 0.0, 0, 1); + run_case_d, CSCMatrix>(Layout::ColMajor, Side::Left, Uplo::Lower, 8, 3, 1.0, 0.0, 0, 1); + run_case_d, COOMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 1.0, 0.0, 0, 1); +} +TEST_F(TestSpsymm, sparse_times_sparse_right_side_matches_reference) { + 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); + run_case_d, CSCMatrix>(Layout::ColMajor, Side::Right, Uplo::Upper, 8, 3, 1.5, -0.5, 0, 1); +} +TEST_F(TestSpsymm, sparse_times_sparse_float_matches_reference) { + run_case_d, CSRMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 1.5f, -0.5f, 0, 1); +} +TEST_F(TestSpsymm, sparse_times_sparse_alpha_zero_scales_by_beta) { + // alpha=0 path: just beta-scales C, 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::spmm(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, symmetric_wrapper_routes_to_same_result) { + 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 + ); +} + + +// MARK: ERROR PATHS + +TEST_F(TestSpsymm, leading_dim_too_small_throws) { + int64_t n_A = 6, d = 3; + std::vector A_tri(n_A * n_A, 0.0); + fill_sym_dense(n_A, A_tri.data(), n_A, 99); + zero_other_triangle(n_A, A_tri.data(), n_A, Uplo::Upper); + CSRMatrix A_sparse(n_A, n_A); + dense_to_sparse_format, double>(Layout::ColMajor, A_tri.data(), 0.0, A_sparse); + + std::vector B(n_A * d, 1.0), C(n_A * d, 0.0); + // side=Left, ColMajor: B and C are n_A-by-d, so ldb, ldc >= n_A. + ASSERT_THROW( + RandBLAS::sparse_data::spsymm(Layout::ColMajor, Side::Left, Uplo::Upper, n_A, d, + 1.0, A_sparse, B.data(), n_A, 0.0, C.data(), n_A - 1), + RandBLAS::Error); + ASSERT_THROW( + RandBLAS::sparse_data::spsymm(Layout::ColMajor, Side::Left, Uplo::Upper, n_A, d, + 1.0, A_sparse, B.data(), n_A - 1, 0.0, C.data(), n_A), + RandBLAS::Error); +} + +TEST_F(TestSpsymm, one_based_indices_throw) { + // The COOMatrix expert constructor accepts IndexBase::One, but spsymm + // requires zero-based indices (as left_spmm does). + double vals[] = {2.0, 1.0, 3.0}; + int64_t rows[] = {1, 1, 2}; + int64_t cols[] = {1, 2, 2}; + COOMatrix A(2, 2, 3, vals, rows, cols, true, + RandBLAS::sparse_data::IndexBase::One); + std::vector B(2 * 2, 1.0), C(2 * 2, 0.0); + ASSERT_THROW( + RandBLAS::sparse_data::spsymm(Layout::ColMajor, Side::Left, Uplo::Upper, 2, 2, + 1.0, A, B.data(), 2, 0.0, C.data(), 2), + RandBLAS::Error); +} + +// Case D with int32 indices: exercises whichever branch the build selects +// (expand-A + spgemm when the index width matches MKL_INT, the densify-B +// composition otherwise). Both must produce the same answer. +TEST_F(TestSpsymm, sparse_times_sparse_int32_indices_match_reference) { + // A_sym = [[2, 1, 0], [1, 3, 0], [0, 0, 4]], upper triangle stored. + double a_vals[] = {2.0, 1.0, 3.0, 4.0}; + int32_t a_rows[] = {0, 0, 1, 2}; + int32_t a_cols[] = {0, 1, 1, 2}; + COOMatrix A(3, 3, 4, a_vals, a_rows, a_cols); + // B = [[1, 0], [0, 5], [2, 0]]. + double b_vals[] = {1.0, 5.0, 2.0}; + int32_t b_rows[] = {0, 1, 2}; + int32_t b_cols[] = {0, 1, 0}; + COOMatrix B(3, 2, 3, b_vals, b_rows, b_cols); + + // A_sym * B = [[2, 5], [1, 15], [8, 0]]. + std::vector C(3 * 2, 0.0); + RandBLAS::sparse_data::spsymm(Layout::ColMajor, Side::Left, Uplo::Upper, 3, 2, + 1.0, A, B, 0.0, C.data(), 3); + std::vector expect = {2.0, 1.0, 8.0, 5.0, 15.0, 0.0}; + for (size_t i = 0; i < expect.size(); ++i) + EXPECT_NEAR(C[i], expect[i], 1e-14) << "mismatch at flat index " << i; +} + +// Empty operands leave beta * C, matching the left_spmm contract from #196 +// (MKL rejects some valid empty sparse matrices at handle creation, so the +// dispatcher must not reach it). +TEST_F(TestSpsymm, empty_sparse_operands_leave_beta_scaled_output) { + int64_t n_A = 4, d = 2; + CSRMatrix A_empty(n_A, n_A); // nnz == 0 + std::vector B(n_A * d, 1.0); + std::vector C(n_A * d, 2.0); + // Case C with structurally empty A: C <- 0.5 * C. + RandBLAS::sparse_data::spsymm(Layout::ColMajor, Side::Left, Uplo::Upper, n_A, d, + 1.0, A_empty, B.data(), n_A, 0.5, C.data(), n_A); + for (auto c : C) EXPECT_DOUBLE_EQ(c, 1.0); + + // Case D with structurally empty B: C <- 0.5 * C again. + COOMatrix B_empty(n_A, d); // nnz == 0 + double a_vals[] = {1.0}; + int64_t a_rows[] = {0}, a_cols[] = {0}; + COOMatrix A_one(n_A, n_A, 1, a_vals, a_rows, a_cols); + RandBLAS::sparse_data::spsymm(Layout::ColMajor, Side::Left, Uplo::Upper, n_A, d, + 1.0, A_one, B_empty, 0.5, C.data(), n_A); + for (auto c : C) EXPECT_DOUBLE_EQ(c, 0.5); +}