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..68ae37a1 100644 --- a/RandBLAS/sksy.hh +++ b/RandBLAS/sksy.hh @@ -32,162 +32,275 @@ #include "RandBLAS/util.hh" #include "RandBLAS/base.hh" #include "RandBLAS/skge.hh" +#include "RandBLAS/sparse_data/coo_sksys_impl.hh" -namespace RandBLAS { -using namespace RandBLAS::dense; -using namespace RandBLAS::sparse; +// ============================================================================= +// Symmetric sketching helpers (SYMM-backed). See project-plans/randblas-symm-plan.md +// for the four-case API design and the implementation status of each case. +// - lsksy3, rsksy3: Case A (dense-symm A x dense Omega), via blas::symm. +// - lsksys, rsksys: Case B (dense-symm A x sparse SkOp). Thin wrappers +// handling the SparseSkOp materialization-and-recurse; the actual COO +// walk + two-axpy scatter lives in sparse_data/coo_sksys_impl.hh +// as coo_lsksys / coo_rsksys. +// ============================================================================= +namespace RandBLAS::dense { + +// ============================================================================= +/// 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, + blas::Uplo uplo, + int64_t d, // B is d-by-n + int64_t n, // A is n-by-n, S is d-by-n (after submat view) + T alpha, + const DenseSkOp &S, + int64_t ro_s, + int64_t co_s, + const T *A, + int64_t lda, + T beta, + T *B, + int64_t ldb +){ + constexpr bool maybe_denseskop = !std::is_same_v, BLASFriendlyOperator>; + if constexpr (maybe_denseskop) { + if (!S.buff) { + auto submat_S = submatrix_as_blackbox>(S, d, n, ro_s, co_s); + lsksy3(layout, uplo, d, n, alpha, submat_S, 0, 0, A, lda, beta, B, ldb); + return; + } + } + randblas_require( S.buff != nullptr ); + randblas_require( S.n_rows >= d + ro_s ); + randblas_require( S.n_cols >= n + co_s ); + if (layout == blas::Layout::ColMajor) { + randblas_require(lda >= n); + randblas_require(ldb >= d); + } else { + randblas_require(lda >= n); + randblas_require(ldb >= n); + } + + auto [pos, lds] = offset_and_ldim(S.layout, S.n_rows, S.n_cols, ro_s, co_s); + T* S_ptr = &S.buff[pos]; + + if (S.layout == layout) { + // Fast path: SYMM directly. + blas::symm(layout, blas::Side::Right, uplo, d, n, + alpha, A, lda, S_ptr, lds, beta, B, ldb); + } else { + // Layout mismatch: transpose-copy S into a tight buffer in the + // caller's layout, then SYMM. Costs O(d*n) for the copy, but keeps + // the SYMM speedup over GEMM on the d*n*n_A matvec. + int64_t lds_new = (layout == blas::Layout::ColMajor) ? d : n; + std::vector S_copy(static_cast(d) * static_cast(n)); + auto [irs_in, ics_in] = layout_to_strides(S.layout, lds); + auto [irs_out, ics_out] = layout_to_strides(layout, lds_new); + util::omatcopy(d, n, S_ptr, irs_in, ics_in, + S_copy.data(), irs_out, ics_out); + blas::symm(layout, blas::Side::Right, uplo, d, n, + alpha, A, lda, S_copy.data(), lds_new, beta, B, ldb); + } + return; +} -// MARK: SUBMAT(S) // ============================================================================= -/// \fn sketch_symmetric(blas::Layout layout, int64_t n, -/// int64_t d, T alpha, const T *A, int64_t lda, -/// const SKOP &S, int64_t ro_s, int64_t co_s, -/// T beta, T *B, int64_t ldb, T sym_check_tol = 0 -/// ) -/// @verbatim embed:rst:leading-slashes -/// Check that :math:`\mat(A)` is symmetric up to tolerance :math:`\texttt{sym_check_tol}`, then sketch from the right in a SYMM-like operation -/// -/// .. math:: -/// \mat(B) = \alpha \cdot \underbrace{\mat(A)}_{n \times n} \cdot \underbrace{\submat(\mtxS)}_{n \times d} + \beta \cdot \underbrace{\mat(B)}_{n \times d}, \tag{$\star$} -/// -/// where :math:`\alpha` and :math:`\beta` are real scalars and :math:`\mtxS` is a sketching operator. -/// -/// .. dropdown:: FAQ -/// :animate: fade-in-slide-down -/// -/// **What's** :math:`\mat(A)?` -/// -/// It's a symmetric matrix of order :math:`n`. Its precise contents depend on :math:`(A, \lda)`, -/// according to -/// -/// .. math:: -/// \mat(A)_{ij} = A[i + j \cdot \lda] = A[i \cdot \lda + j]. -/// -/// Note that the the "layout" parameter passed to this function is not used here. -/// That's because this function requires :math:`\mat(A)` to be stored in the format -/// of a general matrix (with both upper and lower triangles). -/// -/// This function's default behavior is to check that :math:`\mat(A)` is symmetric before -/// attempting sketching. That check can be skipped (at your own peril!) by calling this -/// function with sym_check_tol < 0. -/// -/// **What's** :math:`\mat(B)?` -/// -/// It's an :math:`n \times d` matrix. Its precise contents depend on :math:`(B,\ldb)` and "layout." -/// -/// If layout == ColMajor, then -/// -/// .. math:: -/// \mat(B)_{ij} = B[i + j \cdot \ldb]. -/// -/// In this case, :math:`\ldb` must be :math:`\geq n.` -/// -/// If layout == RowMajor, then -/// -/// .. math:: -/// \mat(B)_{ij} = B[i \cdot \ldb + j]. -/// -/// In this case, :math:`\ldb` must be :math:`\geq d.` -/// -/// **What is** :math:`\submat(\mtxS)` **?** -/// -/// It's the :math:`n \times d` submatrix of :math:`{\mtxS}` whose upper-left corner appears -/// at index :math:`(\texttt{ro_s}, \texttt{co_s})` of :math:`{\mtxS}.` -/// -/// .. dropdown:: Full parameter descriptions -/// :animate: fade-in-slide-down -/// -/// layout - [in] -/// * Either Layout::ColMajor or Layout::RowMajor -/// * Matrix storage for :math:`\mat(B).` -/// -/// n - [in] -/// * A nonnegative integer. -/// * The number of rows in :math:`\mat(B).` -/// * The number of rows and columns in :math:`\mat(A).` -/// -/// d - [in] -/// * A nonnegative integer. -/// * The number of columns in :math:`\mat(B)` and :math:`\submat(\mtxS).` -/// -/// alpha - [in] -/// * A real scalar. -/// * If zero, then :math:`A` is not accessed. -/// -/// A - [in] -/// * Pointer to a 1D array of real scalars. -/// * Defines :math:`\mat(A).` -/// -/// lda - [in] -/// * A nonnegative integer. -/// * Leading dimension of :math:`\mat(A)` when reading from :math:`A.` -/// -/// S - [in] -/// * A DenseSkOp or SparseSkOp object. -/// * Defines :math:`\submat(\mtxS).` -/// -/// ro_s - [in] -/// * A nonnegative integer. -/// * The rows of :math:`\submat(\mtxS)` are a contiguous subset of rows of :math:`S.` -/// * The rows of :math:`\submat(\mtxS)` start at :math:`S[\texttt{ro_s}, :].` -/// -/// co_s - [in] -/// * A nonnnegative integer. -/// * The columns of :math:`\submat(\mtxS)` are a contiguous subset of columns of :math:`S.` -/// * The columns of :math:`\submat(\mtxS)` start at :math:`S[:,\texttt{co_s}].` -/// -/// beta - [in] -/// * A real scalar. -/// * If zero, then :math:`B` need not be set on input. -/// -/// B - [in,out] -/// * Pointer to 1D array of real scalars. -/// * On entry, defines :math:`\mat(B)` -/// on the RIGHT-hand side of :math:`(\star).` -/// * On exit, defines :math:`\mat(B)` -/// on the LEFT-hand side of :math:`(\star).` +/// RSKSY3: SYMM-backed right-sketch with a symmetric matrix A. /// -/// ldb - [in] -/// * A nonnegative integer. -/// * Leading dimension of :math:`\mat(B)` when reading from :math:`B.` +/// Computes B = alpha * mat(A) * submat(S) + beta * B, where: +/// - mat(A) is n-by-n symmetric. Only the triangle named by `uplo` is read. +/// - submat(S) is the n-by-d view of S at (ro_s, co_s). +/// - mat(B) is n-by-d. /// -/// @endverbatim -template -inline void sketch_symmetric( - // B = alpha*A*S + beta*B, where A is a symmetric matrix stored in the format of a general matrix. +/// Same materialization and layout-mismatch fallback semantics as `lsksy3`. +/// Final call (matching layout): `blas::symm` with `side = Left`. +template +void rsksy3( blas::Layout layout, - int64_t n, // number of rows in B - int64_t d, // number of columns in B + blas::Uplo uplo, + int64_t n, // A is n-by-n, S is n-by-d (after submat view), B is n-by-d + int64_t d, T alpha, - const T* A, + const T *A, int64_t lda, - const SKOP &S, + const DenseSkOp &S, int64_t ro_s, int64_t co_s, T beta, - T* B, - int64_t ldb, - T sym_check_tol = 0 + T *B, + int64_t ldb +){ + constexpr bool maybe_denseskop = !std::is_same_v, BLASFriendlyOperator>; + if constexpr (maybe_denseskop) { + if (!S.buff) { + auto submat_S = submatrix_as_blackbox>(S, n, d, ro_s, co_s); + rsksy3(layout, uplo, n, d, alpha, A, lda, submat_S, 0, 0, beta, B, ldb); + return; + } + } + randblas_require( S.buff != nullptr ); + randblas_require( S.n_rows >= n + ro_s ); + randblas_require( S.n_cols >= d + co_s ); + if (layout == blas::Layout::ColMajor) { + randblas_require(lda >= n); + randblas_require(ldb >= n); + } else { + randblas_require(lda >= n); + randblas_require(ldb >= d); + } + + auto [pos, lds] = offset_and_ldim(S.layout, S.n_rows, S.n_cols, ro_s, co_s); + T* S_ptr = &S.buff[pos]; + + if (S.layout == layout) { + blas::symm(layout, blas::Side::Left, uplo, n, d, + alpha, A, lda, S_ptr, lds, beta, B, ldb); + } else { + // Layout mismatch: transpose-copy S into a tight matching-layout + // buffer, then SYMM. See lsksy3 for the trade-off discussion. + int64_t lds_new = (layout == blas::Layout::ColMajor) ? n : d; + std::vector S_copy(static_cast(n) * static_cast(d)); + auto [irs_in, ics_in] = layout_to_strides(S.layout, lds); + auto [irs_out, ics_out] = layout_to_strides(layout, lds_new); + util::omatcopy(n, d, S_ptr, irs_in, ics_in, + S_copy.data(), irs_out, ics_out); + blas::symm(layout, blas::Side::Left, uplo, n, d, + alpha, A, lda, S_copy.data(), lds_new, beta, B, ldb); + } + return; +} + +} // end namespace RandBLAS::dense + + +namespace RandBLAS::sparse { + +// ============================================================================= +/// LSKSYS: dense symmetric A on the right of a SparseSkOp. +/// +/// Computes B = alpha * submat(S) * mat(A) + beta * B, where: +/// - mat(A) is n-by-n dense symmetric, with only the `uplo` triangle stored. +/// - submat(S) is the d-by-n view of S at (ro_s, co_s); S is a SparseSkOp. +/// - mat(B) is d-by-n dense. +/// +/// Inner loop: for each stored nonzero (row_S, col_S, v) of S inside the +/// submatrix window, contribute alpha*v * (row col_S of A) to row (row_S - ro_s) +/// of B. The row of symmetric A is assembled from the stored triangle as two +/// blas::axpy calls (one along the stored column, one along the stored row past +/// the diagonal). +template +void lsksys( + blas::Layout layout, + blas::Uplo uplo, + int64_t d, + int64_t n, + T alpha, + const SparseSkOp &S, + int64_t ro_s, + int64_t co_s, + const T *A, + int64_t lda, + T beta, + T *B, + int64_t ldb ) { - 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); + if (S.nnz < 0) { + SparseSkOp shallowcopy(S.dist, S.seed_state); + fill_sparse(shallowcopy); + lsksys(layout, uplo, d, n, alpha, shallowcopy, ro_s, co_s, A, lda, beta, B, ldb); + return; + } + + util::lascl(layout, d, n, beta, B, ldb); + 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 + ); } // ============================================================================= -/// \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 -/// ) +/// RSKSYS: dense symmetric A on the left of a SparseSkOp. +/// +/// Computes B = alpha * mat(A) * submat(S) + beta * B, where: +/// - mat(A) is n-by-n dense symmetric, with only the `uplo` triangle stored. +/// - submat(S) is the n-by-d view of S at (ro_s, co_s); S is a SparseSkOp. +/// - mat(B) is n-by-d dense. +/// +/// Same two-axpy scatter pattern as lsksys, but on columns of A (the column +/// of A indexed by row_S of S contributes into the column of B indexed by +/// col_S - co_s). +template +void rsksys( + blas::Layout layout, + blas::Uplo uplo, + int64_t n, + int64_t d, + T alpha, + const T *A, + int64_t lda, + const SparseSkOp &S, + int64_t ro_s, + int64_t co_s, + T beta, + T *B, + int64_t ldb +) { + if (S.nnz < 0) { + SparseSkOp shallowcopy(S.dist, S.seed_state); + fill_sparse(shallowcopy); + rsksys(layout, uplo, n, d, alpha, A, lda, shallowcopy, ro_s, co_s, beta, B, ldb); + return; + } + + util::lascl(layout, n, d, beta, B, ldb); + 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 +308,28 @@ inline void sketch_symmetric( /// /// **What's** :math:`\mat(A)?` /// -/// It's a symmetric matrix of order :math:`n`. Its precise contents depend on :math:`(A, \lda)`, -/// according to -/// -/// .. math:: -/// \mat(A)_{ij} = A[i + j \cdot \lda] = A[i \cdot \lda + j]. -/// -/// Note that the the "layout" parameter passed to this function is not used here. -/// That's because this function requires :math:`\mat(A)` to be stored in the format -/// of a general matrix (with both upper and lower triangles). -/// -/// This function's default behavior is to check that :math:`\mat(A)` is symmetric before -/// attempting sketching. That check can be skipped (at your own peril!) by calling this -/// function with sym_check_tol < 0. -/// -/// **What's** :math:`\mat(B)?` -/// -/// It's a :math:`d \times n` matrix. Its precise contents depend on :math:`(B,\ldb)` and "layout." +/// 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. /// -/// If layout == ColMajor, then +/// The :math:`\texttt{layout}` parameter governs the indexing convention into :math:`A`: /// /// .. 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.` +/// \mat(A)_{ij} = A[i + j \cdot \lda] \quad (\text{ColMajor}) +/// = A[i \cdot \lda + j] \quad (\text{RowMajor}). /// -/// **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 +337,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 hand-rolled +/// per-stored-nonzero scatter that reads only the named triangle of +/// :math:`A` (Case B). See `project-plans/randblas-symm-plan.md`. /// /// ro_s - [in] /// * A nonnegative integer. -/// * The rows of :math:`\submat(\mtxS)` are a contiguous subset of rows of :math:`S.` -/// * The rows of :math:`\submat(\mtxS)` start at :math:`S[\texttt{ro_s}, :].` +/// * The rows of :math:`\submat(\mtxS)` start at :math:`\mtxS[\texttt{ro_s}, :].` /// /// co_s - [in] -/// * A nonnnegative integer. -/// * The columns of :math:`\submat(\mtxS)` are a contiguous subset of columns of :math:`S.` -/// * The columns of :math:`\submat(\mtxS)` start at :math:`S[:,\texttt{co_s}].` +/// * A nonnegative integer. +/// * The columns of :math:`\submat(\mtxS)` start at :math:`\mtxS[:, \texttt{co_s}].` /// /// A - [in] /// * Pointer to a 1D array of real scalars. @@ -278,255 +371,227 @@ inline void sketch_symmetric( /// * A real scalar. /// * If zero, then :math:`B` need not be set on input. /// -/// B - [in,out] +/// B - [in, out] /// * Pointer to 1D array of real scalars. -/// * On entry, defines :math:`\mat(B)` -/// on the RIGHT-hand side of :math:`(\star).` -/// * On exit, defines :math:`\mat(B)` -/// on the LEFT-hand side of :math:`(\star).` +/// * On entry, defines :math:`\mat(B)` on the RIGHT-hand side of :math:`(\star).` +/// * On exit, defines :math:`\mat(B)` on the LEFT-hand side of :math:`(\star).` /// /// ldb - [in] /// * A nonnegative integer. /// * Leading dimension of :math:`\mat(B)` when reading from :math:`B.` /// /// @endverbatim -template +template inline void sketch_symmetric( - // B = alpha*S*A + beta*B - blas::Layout layout, - int64_t d, // number of rows in B - int64_t n, // number of columns in B + blas::Layout layout, blas::Uplo uplo, + int64_t d, int64_t n, T alpha, - const SKOP &S, - int64_t ro_s, - int64_t co_s, - const T* A, - int64_t lda, + const SKOP &S, int64_t ro_s, int64_t co_s, + const T* A, int64_t lda, T beta, - T* B, - int64_t ldb, - T sym_check_tol = 0 + T* B, int64_t ldb +); + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + int64_t d, int64_t n, + T alpha, + const DenseSkOp &S, int64_t ro_s, int64_t co_s, + const T* A, int64_t lda, + T beta, + T* B, int64_t ldb ) { - RandBLAS::util::require_symmetric(layout, A, n, lda, sym_check_tol); - sketch_general(layout, blas::Op::NoTrans, blas::Op::NoTrans, d, n, n, alpha, S, ro_s, co_s, A, lda, beta, B, ldb); + RandBLAS::dense::lsksy3(layout, uplo, d, n, alpha, S, ro_s, co_s, A, lda, beta, B, ldb); +} + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + int64_t d, int64_t n, + T alpha, + const SparseSkOp &S, int64_t ro_s, int64_t co_s, + const T* A, int64_t lda, + T beta, + T* B, int64_t ldb +) { + RandBLAS::sparse::lsksys(layout, uplo, d, n, alpha, S, ro_s, co_s, A, lda, beta, B, ldb); } -// 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). -/// -/// 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.` +/// Sketch from the right in a SYMM-like operation /// -/// If layout == RowMajor, then -/// -/// .. math:: -/// \mat(B)_{ij} = B[i \cdot \ldb + j]. -/// -/// In this case, :math:`\ldb` must be :math:`\geq d.` -/// -/// .. dropdown:: Full parameter descriptions -/// :animate: fade-in-slide-down -/// -/// layout - [in] -/// * Either Layout::ColMajor or Layout::RowMajor -/// * Matrix storage for :math:`\mat(B).` -/// -/// alpha - [in] -/// * A real scalar. -/// * If zero, then :math:`A` is not accessed. +/// .. math:: +/// \mat(B) = \alpha \cdot \underbrace{\mat(A)}_{n \times n} \cdot \underbrace{\submat(\mtxS)}_{n \times d} + \beta \cdot \underbrace{\mat(B)}_{n \times d}, \tag{$\star$} /// -/// A - [in] -/// * Pointer to a 1D array of real scalars. -/// * Defines :math:`\mat(A).` +/// where :math:`\alpha` and :math:`\beta` are real scalars and :math:`\mtxS` is a sketching operator. /// -/// lda - [in] -/// * A nonnegative integer. -/// * Leading dimension of :math:`\mat(A)` when reading from :math:`A.` +/// See the left-side overload above for the meaning of :math:`\mat(A)` and its +/// :math:`\texttt{uplo}` storage convention. The roles of :math:`d` and :math:`n` are mirrored: +/// :math:`d` is the embedding dimension (cols of :math:`\submat(\mtxS)`) and :math:`n` is the +/// order of :math:`\mat(A).` /// -/// S - [in] -/// * A DenseSkOp or SparseSkOp object. +/// @endverbatim +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + int64_t n, int64_t d, + T alpha, + const T* A, int64_t lda, + const SKOP &S, int64_t ro_s, int64_t co_s, + T beta, + T* B, int64_t ldb +); + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + int64_t n, int64_t d, + T alpha, + const T* A, int64_t lda, + const DenseSkOp &S, int64_t ro_s, int64_t co_s, + T beta, + T* B, int64_t ldb +) { + RandBLAS::dense::rsksy3(layout, uplo, n, d, alpha, A, lda, S, ro_s, co_s, beta, B, ldb); +} + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + int64_t n, int64_t d, + T alpha, + const T* A, int64_t lda, + const SparseSkOp &S, int64_t ro_s, int64_t co_s, + T beta, + T* B, int64_t ldb +) { + RandBLAS::sparse::rsksys(layout, uplo, n, d, alpha, A, lda, S, ro_s, co_s, beta, B, ldb); +} + + +// MARK: FULL(S) + +// ============================================================================= +/// \fn sketch_symmetric(blas::Layout layout, blas::Uplo uplo, T alpha, +/// const SKOP &S, const T *A, int64_t lda, +/// T beta, T *B, int64_t ldb +/// ) +/// @verbatim embed:rst:leading-slashes +/// Sketch from the left in a SYMM-like operation, with :math:`\mtxS` used in full +/// (no submatrix offsets). /// -/// beta - [in] -/// * A real scalar. -/// * If zero, then :math:`B` need not be set on input. +/// .. math:: +/// \mat(B) = \alpha \cdot \underbrace{\mtxS}_{d \times n} \cdot \underbrace{\mat(A)}_{n \times n} + \beta \cdot \underbrace{\mat(B)}_{d \times n}, \tag{$\star$} /// -/// B - [in,out] -/// * Pointer to 1D array of real scalars. -/// * On entry, defines :math:`\mat(B)` -/// on the RIGHT-hand side of :math:`(\star).` -/// * On exit, defines :math:`\mat(B)` -/// on the LEFT-hand side of :math:`(\star).` +/// The dimensions :math:`d` and :math:`n` are taken from :math:`\mtxS` directly +/// (:math:`d = \mtxS.\ttt{dist.n\_rows}`, :math:`n = \mtxS.\ttt{dist.n\_cols}`). /// -/// ldb - [in] -/// * A nonnegative integer. -/// * Leading dimension of :math:`\mat(B)` when reading from :math:`B.` +/// See the submatrix overload above for the meaning of :math:`\mat(A)` and +/// :math:`\texttt{uplo}`. /// /// @endverbatim -template +template inline void sketch_symmetric( - // B = alpha*A*S + beta*B, where A is a symmetric matrix stored in the format of a general matrix. - blas::Layout layout, + blas::Layout layout, blas::Uplo uplo, T alpha, - const T* A, - int64_t lda, const SKOP &S, + const T* A, int64_t lda, T beta, - T* B, - int64_t ldb, - T sym_check_tol = 0 + T* B, int64_t ldb +); + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + T alpha, + const DenseSkOp &S, + const T* A, int64_t lda, + T beta, + T* B, int64_t ldb ) { - int64_t n = S.dist.n_rows; - int64_t d = S.dist.n_cols; - RandBLAS::util::require_symmetric(layout, A, n, lda, sym_check_tol); - sketch_general(layout, blas::Op::NoTrans, blas::Op::NoTrans, n, d, n, alpha, A, lda, S, 0, 0, beta, B, ldb); + int64_t d = S.dist.n_rows; + int64_t n = S.dist.n_cols; + RandBLAS::dense::lsksy3(layout, uplo, d, n, alpha, S, 0, 0, A, lda, beta, B, ldb); +} + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + T alpha, + const SparseSkOp &S, + const T* A, int64_t lda, + T beta, + T* B, int64_t ldb +) { + int64_t d = S.dist.n_rows; + int64_t n = S.dist.n_cols; + RandBLAS::sparse::lsksys(layout, uplo, d, n, alpha, S, 0, 0, A, lda, beta, B, ldb); } // ============================================================================= -/// \fn sketch_symmetric(blas::Layout layout, T alpha, const SKOP &S, -/// const T *A, int64_t lda, T beta, T *B, int64_t ldb, T sym_check_tol = 0 -/// ) +/// \fn sketch_symmetric(blas::Layout layout, blas::Uplo uplo, T alpha, +/// const T *A, int64_t lda, const SKOP &S, +/// T beta, T *B, int64_t ldb +/// ) /// @verbatim embed:rst:leading-slashes -/// Check that :math:`\mat(A)` is symmetric up to tolerance :math:`\texttt{sym_check_tol}`, then sketch from the left in a SYMM-like operation -/// -/// .. math:: -/// \mat(B) = \alpha \cdot \mtxS \cdot \underbrace{\mat(A)}_{n \times n} + \beta \cdot \underbrace{\mat(B)}_{d \times n}, \tag{$\star$} -/// -/// where :math:`\alpha` and :math:`\beta` are real scalars and :math:`\mtxS` is a :math:`d \times n` sketching operator. -/// -/// .. dropdown:: FAQ -/// :animate: fade-in-slide-down -/// -/// **What's** :math:`\mat(A)?` -/// -/// It's a symmetric matrix of order :math:`n`. Its precise contents depend on :math:`(A, \lda)`, -/// according to -/// -/// .. math:: -/// \mat(A)_{ij} = A[i + j \cdot \lda] = A[i \cdot \lda + j]. -/// -/// Note that the the "layout" parameter passed to this function is not used here. -/// That's because this function requires :math:`\mat(A)` to be stored in the format -/// of a general matrix (with both upper and lower triangles). -/// -/// This function's default behavior is to check that :math:`\mat(A)` is symmetric before -/// attempting sketching. That check can be skipped (at your own peril!) by calling this -/// function with sym_check_tol < 0. -/// -/// **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).` +/// Sketch from the right in a SYMM-like operation, with :math:`\mtxS` used in full. /// -/// alpha - [in] -/// * A real scalar. -/// * If zero, then :math:`A` is not accessed. -/// -/// S - [in] -/// * A DenseSkOp or SparseSkOp object. -/// -/// A - [in] -/// * Pointer to a 1D array of real scalars. -/// * Defines :math:`\mat(A).` -/// -/// lda - [in] -/// * A nonnegative integer. -/// * Leading dimension of :math:`\mat(A)` when reading from :math:`A.` +/// .. math:: +/// \mat(B) = \alpha \cdot \underbrace{\mat(A)}_{n \times n} \cdot \underbrace{\mtxS}_{n \times d} + \beta \cdot \underbrace{\mat(B)}_{n \times d}, \tag{$\star$} /// -/// beta - [in] -/// * A real scalar. -/// * If zero, then :math:`B` need not be set on input. +/// The dimensions :math:`n` and :math:`d` are taken from :math:`\mtxS` directly +/// (:math:`n = \mtxS.\ttt{dist.n\_rows}`, :math:`d = \mtxS.\ttt{dist.n\_cols}`). /// -/// B - [in,out] -/// * Pointer to 1D array of real scalars. -/// * On entry, defines :math:`\mat(B)` -/// on the RIGHT-hand side of :math:`(\star).` -/// * On exit, defines :math:`\mat(B)` -/// on the LEFT-hand side of :math:`(\star).` -/// -/// ldb - [in] -/// * A nonnegative integer. -/// * Leading dimension of :math:`\mat(B)` when reading from :math:`B.` +/// See the submatrix overload above for the meaning of :math:`\mat(A)` and +/// :math:`\texttt{uplo}`. /// /// @endverbatim -template +template inline void sketch_symmetric( - // B = alpha*S*A + beta*B - blas::Layout layout, + blas::Layout layout, blas::Uplo uplo, T alpha, + const T* A, int64_t lda, const SKOP &S, - const T* A, - int64_t lda, T beta, - T* B, - int64_t ldb, - T sym_check_tol = 0 + T* B, int64_t ldb +); + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + T alpha, + const T* A, int64_t lda, + const DenseSkOp &S, + T beta, + T* B, int64_t ldb ) { - int64_t d = S.dist.n_rows; - int64_t n = S.dist.n_cols; - RandBLAS::util::require_symmetric(layout, A, n, lda, sym_check_tol); - sketch_general(layout, blas::Op::NoTrans, blas::Op::NoTrans, d, n, n, alpha, S, 0, 0, A, lda, beta, B, ldb); + int64_t n = S.dist.n_rows; + int64_t d = S.dist.n_cols; + RandBLAS::dense::rsksy3(layout, uplo, n, d, alpha, A, lda, S, 0, 0, beta, B, ldb); +} + +template +inline void sketch_symmetric( + blas::Layout layout, blas::Uplo uplo, + T alpha, + const T* A, int64_t lda, + const SparseSkOp &S, + T beta, + T* B, int64_t ldb +) { + int64_t n = S.dist.n_rows; + int64_t d = S.dist.n_cols; + RandBLAS::sparse::rsksys(layout, uplo, n, d, alpha, A, lda, S, 0, 0, beta, B, ldb); } } // end namespace RandBLAS diff --git a/RandBLAS/sparse_data/DevNotes.md b/RandBLAS/sparse_data/DevNotes.md index 24d52336..3ba86cb7 100644 --- a/RandBLAS/sparse_data/DevNotes.md +++ b/RandBLAS/sparse_data/DevNotes.md @@ -54,6 +54,118 @@ either ``lskges`` or ``rskges``. Here's what would happen after we entered one o (if inside ``lskges``) or ``right_spmm`` (if inside ``rskges``). +## SYMM-shaped kernels (spsymm) + +RandBLAS exposes a SYMM-style API for sparse symmetric matrices via the +``spsymm`` family. The design covers four cases based on the structure of the +two operands (``A`` symmetric vs. the second factor ``B``): + +| Tag | Operation | A storage | B storage | Status this PR | +|-----|---------------------------------|---------------------|-----------|-------------------------------------------------| +| A | dense-symm × dense | dense, one triangle | dense | Implemented via ``blas::symm`` in ``sksy.hh``. | +| B | dense-symm × sparse | dense, one triangle | sparse | Implemented via hand-rolled ``lsksys`` / ``rsksys`` wrappers in ``sksy.hh`` that handle the SparseSkOp materialization-and-recurse, and the ``coo_lsksys`` / ``coo_rsksys`` COO kernels in ``sparse_data/coo_sksys_impl.hh`` that do the two-axpy scatter per stored nonzero of S (reading only the named triangle of A). | +| C | sparse-symm × dense (→ dense) | sparse, one triangle | dense | Implemented in ``spsymm_dispatch.hh`` (MKL fast path + per-format fallbacks). | +| D | sparse-symm × sparse → dense | sparse, one triangle | sparse | Implemented via densify-B + Case-C composition in ``spsymm_dispatch.hh``. | + +### 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 (mostly) | ``mkl_sparse_d_mm`` with ``descr.type = SPARSE_MATRIX_TYPE_SYMMETRIC``. RandBLAS falls back to a hand kernel for side=Right (MKL has no Side parameter) and for CSC (``mkl_sparse_d_mm`` returns NOT_SUPPORTED on CSC). | +| D | No | ``mkl_sparse_sp2m`` returns ``SPARSE_STATUS_NOT_SUPPORTED`` when ``descrA.type == SPARSE_MATRIX_TYPE_SYMMETRIC`` (only ``GENERAL`` is accepted there); ``mkl_sparse_d_spmmd`` takes no descriptor at all. Symmetric expansion has to happen on the RandBLAS side, so we don't gain anything by routing through MKL. | + +### Case C dispatch (``spsymm_dispatch.hh``) + +``RandBLAS::sparse_data::spsymm(layout, side, uplo, m, n, alpha, A, B, ldb, beta, Y, ldy)`` +dispatches as follows: + +1. Validate: A is square (``A.n_rows == A.n_cols``), and matches the side + convention (``A.n_rows == m`` for side=Left, ``A.n_rows == n`` for side=Right). +2. If RandBLAS was built with MKL and the index width matches ``MKL_INT``, + try ``mkl::mkl_spsymm``. It applies the ``SPARSE_MATRIX_TYPE_SYMMETRIC`` + descriptor and calls ``mkl_sparse_d_mm`` directly. Returns false for + side=Right and CSC; control falls through to step 3 in either case. +3. Format-specific fallback: ``csr_spsymm`` / ``csc_spsymm`` / ``coo_spsymm``. + Each iterates the named triangle once. For each stored entry ``A(i,j) = v``, + it emits one ``blas::axpy`` for the structural location and a second one + for the implied symmetric counterpart (when ``i != j``). Diagonal entries + contribute once. Entries outside the named triangle are silently skipped, + so a caller that mistakenly stored both triangles still gets the correct + answer (the kernel just behaves as if the "extra" entries were absent). + +A shared ``internal::apply_beta_scale`` helper (defined in +``csr_spsymm_impl.hh``, re-included by the other two format files) handles +the ``Y <- beta * Y`` pass on entry. + +The public-facing wrappers in the top-level ``RandBLAS::`` namespace are: + + - ``spsymm(layout, uplo, m, n, alpha, A, B, ldb, beta, Y, ldy)``, + convenience for side=Left. + - ``spsymm(layout, m, n, alpha, Symmetric A_sym, B, ldb, beta, Y, ldy)``, + routes via the ``Symmetric`` carrier so the uplo annotation + travels with the matrix. + +### Case B: hand-rolled COO kernels in ``coo_sksys_impl.hh`` + +The COO walk and scatter live in ``sparse_data/coo_sksys_impl.hh`` as +``coo_lsksys`` / ``coo_rsksys``, taking a ``COOMatrix`` plus +submatrix offsets and writing into a dense buffer. Thin wrappers +``lsksys`` / ``rsksys`` in ``sksy.hh`` handle the SparseSkOp +materialization-and-recurse pattern, beta-scale ``B`` via ``util::lascl``, +unpack the SparseSkOp into its COO view, and forward into the kernel. +The split keeps ``sksy.hh`` focused on SkOp dispatch and puts the +format-specific work next to the other COO kernels. + +For each stored nonzero ``(i_S, j_S, v)`` of the SparseSkOp's COO view (with +submatrix offsets ``(ro_s, co_s)`` filtered inline), the kernel applies +``alpha * v`` to one row or column of the symmetric dense A and accumulates +into the corresponding row or column of the output. Reading "row j of A" (or +"col i of A") splits into two ranges based on the diagonal: + + - For ``Uplo::Upper``: the part above the diagonal walks A's stored row / + column directly; the part below comes from the symmetric reflection + (reading the transposed-position entry from the stored triangle). + - For ``Uplo::Lower``: the roles swap. + +Each range becomes a single ``blas::axpy`` with the appropriate stride, so +the inner-loop body is exactly two AXPY calls per stored nonzero of S +(plus a uniform layout / uplo branch). No special handling for the diagonal +beyond consistent inclusion in one of the two ranges. SparseSkOp is COO +internally, so submatrix filtering is a direct ``if (row < ro_s ...) continue`` +on the COO triples. + +### Case D: densify-B + Case-C composition + +Lives in ``spsymm_dispatch.hh`` next to the Case-C dispatcher. The +overload taking two ``SparseMatrix`` operands allocates an ``m`` by +``n`` ``std::vector`` (tight leading dim in the caller's +``layout``), fills it via the format-specific ``coo_to_dense`` / +``csr_to_dense`` / ``csc_to_dense`` helper picked by ``if constexpr``, +then calls the existing Case-C ``spsymm`` overload on the densified +buffer. Works in any build (the MKL fast path or the per-format hand +kernel both apply, depending on the build), and across all 3 × 3 = 9 +sparse-format pairings for ``(A, B)`` since the densification picks +the right format-specific helper. + +Why this composition rather than a single MKL ``sp2m`` call: + + - ``mkl_sparse_sp2m`` returns ``SPARSE_STATUS_NOT_SUPPORTED`` when + the ``matrix_descr`` on either operand is + ``SPARSE_MATRIX_TYPE_SYMMETRIC``; only ``GENERAL`` is accepted + there. + - ``mkl_sparse_d_spmmd`` (which writes directly to dense ``C``) + accepts no descriptor at all. + +So the symmetric expansion has to happen on the RandBLAS side either +way. Composing through Case C gets it for free at the cost of a +temporary dense buffer for ``B`` (cost ``O(m*n)``), and for the +typical RandNLA workload where ``B`` is a sketching operator with +``nnz(B) << m*n`` the buffer cost is small relative to the work that +would have to happen anyway. ``Y`` itself is never touched until the +Case-C call. + ## Sketching sparse data with dense operators 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 7eb4df2a..8926b81d 100644 --- a/RandBLAS/sparse_data/base.hh +++ b/RandBLAS/sparse_data/base.hh @@ -31,6 +31,7 @@ #include "RandBLAS/config.h" #include "RandBLAS/base.hh" #include +#include // std::iota, used below #ifdef __cpp_concepts #include diff --git a/RandBLAS/sparse_data/coo_sksys_impl.hh b/RandBLAS/sparse_data/coo_sksys_impl.hh new file mode 100644 index 00000000..7aa82dd2 --- /dev/null +++ b/RandBLAS/sparse_data/coo_sksys_impl.hh @@ -0,0 +1,184 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include "RandBLAS/base.hh" +#include "RandBLAS/util.hh" +#include "RandBLAS/sparse_data/base.hh" +#include "RandBLAS/sparse_data/coo_matrix.hh" + +#include + + +namespace RandBLAS::sparse_data { + +// ============================================================================= +// COO sparse-from-left times dense symmetric A into dense B: +// B = alpha * submat(Scoo) * sym(A, uplo) + beta * B +// +// - sym(A, uplo) is n-by-n dense symmetric, only the `uplo` triangle stored. +// - submat(Scoo) is the d-by-n window of Scoo at (ro_s, co_s). +// - B is d-by-n dense, layout-matched. +// +// For each stored nonzero (row_S, col_S, v) of Scoo inside the window, the +// kernel contributes alpha*v * (row col_S of A) to row (row_S - ro_s) of B. +// The row of symmetric A is assembled from the stored triangle as two +// blas::axpy calls: one along the stored column up to (and including) the +// diagonal, one along the stored row past the diagonal. The split is chosen +// by `uplo` and the matrix layout. +// +// Beta-scaling of B and the alpha==0 short-circuit are the caller's +// responsibility (so that this kernel can be composed with other accumulators +// without redundant scaling). +// ============================================================================= +template +void coo_lsksys( + blas::Layout layout, + blas::Uplo uplo, + int64_t d, + int64_t n, + T alpha, + const COOMatrix &Scoo, + int64_t ro_s, + int64_t co_s, + const T *A, + int64_t lda, + T *B, + int64_t ldb +) { + if (alpha == T(0)) return; + + const bool col_major = (layout == blas::Layout::ColMajor); + + for (int64_t p = 0; p < Scoo.nnz; ++p) { + sint_t row_S = Scoo.rows[p]; + sint_t col_S = Scoo.cols[p]; + if (row_S < ro_s || row_S >= ro_s + d) continue; + if (col_S < co_s || col_S >= co_s + n) continue; + int64_t i = static_cast(row_S) - ro_s; // B row + int64_t j = static_cast(col_S) - co_s; // A row index (sym A read as row j) + T av = alpha * Scoo.vals[p]; + + // B[i, :] += av * row_j(A_sym), split by uplo into two contiguous + // ranges of A so each becomes a single axpy. + if (uplo == blas::Uplo::Upper) { + // row j of A: + // c in [0, j]: read A[c, j] (column j of stored Upper, rows 0..j) + // c in (j, n-1]: read A[j, c] (row j of stored Upper, cols j+1..n-1) + if (col_major) { + blas::axpy(j + 1, av, &A[j * lda], 1, &B[i], ldb); + blas::axpy(n - j - 1, av, &A[j + (j + 1) * lda], lda, &B[i + (j + 1) * ldb], ldb); + } else { + blas::axpy(j + 1, av, &A[j], lda, &B[i * ldb], 1); + blas::axpy(n - j - 1, av, &A[j * lda + j + 1], 1, &B[i * ldb + j + 1], 1); + } + } else { + // Lower-stored: + // c in [0, j]: read A[j, c] (row j of stored Lower, cols 0..j) + // c in (j, n-1]: read A[c, j] (column j of stored Lower, rows j+1..n-1) + if (col_major) { + blas::axpy(j + 1, av, &A[j], lda, &B[i], ldb); + blas::axpy(n - j - 1, av, &A[(j + 1) + j * lda], 1, &B[i + (j + 1) * ldb], ldb); + } else { + blas::axpy(j + 1, av, &A[j * lda], 1, &B[i * ldb], 1); + blas::axpy(n - j - 1, av, &A[(j + 1) * lda + j], lda, &B[i * ldb + j + 1], 1); + } + } + } +} + + +// ============================================================================= +// COO sparse-from-right times dense symmetric A into dense B: +// B = alpha * sym(A, uplo) * submat(Scoo) + beta * B +// +// Same two-axpy scatter pattern as coo_lsksys, but on columns of A: the +// column of A indexed by row_S of S contributes into the column of B indexed +// by col_S - co_s. +// +// Beta-scaling of B and the alpha==0 short-circuit are the caller's +// responsibility. +// ============================================================================= +template +void coo_rsksys( + blas::Layout layout, + blas::Uplo uplo, + int64_t n, + int64_t d, + T alpha, + const T *A, + int64_t lda, + const COOMatrix &Scoo, + int64_t ro_s, + int64_t co_s, + T *B, + int64_t ldb +) { + if (alpha == T(0)) return; + + const bool col_major = (layout == blas::Layout::ColMajor); + + for (int64_t p = 0; p < Scoo.nnz; ++p) { + sint_t row_S = Scoo.rows[p]; + sint_t col_S = Scoo.cols[p]; + if (row_S < ro_s || row_S >= ro_s + n) continue; + if (col_S < co_s || col_S >= co_s + d) continue; + int64_t i = static_cast(row_S) - ro_s; // A column index + int64_t j = static_cast(col_S) - co_s; // B column + T av = alpha * Scoo.vals[p]; + + // B[:, j] += av * col_i(A_sym), split by uplo: + if (uplo == blas::Uplo::Upper) { + // col i of A: + // r in [0, i]: read A[r, i] (column i of stored Upper, rows 0..i) + // r in (i, n-1]: read A[i, r] (row i of stored Upper, cols i+1..n-1) + if (col_major) { + blas::axpy(i + 1, av, &A[i * lda], 1, &B[j * ldb], 1); + blas::axpy(n - i - 1, av, &A[i + (i + 1) * lda], lda, &B[(i + 1) + j * ldb], 1); + } else { + blas::axpy(i + 1, av, &A[i], lda, &B[j], ldb); + blas::axpy(n - i - 1, av, &A[i * lda + i + 1], 1, &B[(i + 1) * ldb + j], ldb); + } + } else { + // Lower-stored col i: + // r in [0, i-1]: read A[i, r] (row i of stored Lower, cols 0..i-1) + // r in [i, n-1]: read A[r, i] (column i of stored Lower, rows i..n-1) + if (col_major) { + blas::axpy(i, av, &A[i], lda, &B[j * ldb], 1); + blas::axpy(n - i, av, &A[i + i * lda], 1, &B[i + j * ldb], 1); + } else { + blas::axpy(i, av, &A[i * lda], 1, &B[j], ldb); + blas::axpy(n - i, av, &A[i * lda + i], lda, &B[i * ldb + j], ldb); + } + } + } +} + +} // end namespace RandBLAS::sparse_data diff --git a/RandBLAS/sparse_data/coo_spsymm_impl.hh b/RandBLAS/sparse_data/coo_spsymm_impl.hh new file mode 100644 index 00000000..31f5a7a5 --- /dev/null +++ b/RandBLAS/sparse_data/coo_spsymm_impl.hh @@ -0,0 +1,80 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include "RandBLAS/exceptions.hh" +#include "RandBLAS/util.hh" +#include "RandBLAS/sparse_data/base.hh" +#include "RandBLAS/sparse_data/coo_matrix.hh" +#include "RandBLAS/sparse_data/spsymm_internal.hh" +#include + +namespace RandBLAS::sparse_data { + +// ============================================================================= +/// COO fallback for symmetric sparse-times-dense. +/// +/// Iterates over the (row, col, val) triples directly. For uplo=Upper, +/// structurally populated entries satisfy row <= col; for uplo=Lower, +/// row >= col. Entries outside the named triangle are silently skipped. +/// No assumption on the order of the COO entries. +template +void coo_spsymm( + blas::Layout layout, + blas::Side side, + blas::Uplo uplo, + int64_t m, int64_t n, + T alpha, + const COOMatrix& A, + const T* B, int64_t ldb, + T beta, + T* Y, int64_t ldy +) { + randblas_require(A.n_rows == A.n_cols); + int64_t k = (side == blas::Side::Left) ? m : n; + randblas_require(A.n_rows == k); + + RandBLAS::util::lascl(layout, m, n, beta, Y, ldy); + if (alpha == T(0)) return; + + for (int64_t p = 0; p < A.nnz; ++p) { + sint_t i = A.rows[p]; + sint_t j = A.cols[p]; + if (uplo == blas::Uplo::Upper && j < i) continue; + if (uplo == blas::Uplo::Lower && j > i) continue; + T av = alpha * A.vals[p]; + if (side == blas::Side::Left) + internal::spsymm_scatter_left(layout, n, av, i, j, B, ldb, Y, ldy); + else + internal::spsymm_scatter_right(layout, m, av, i, j, B, ldb, Y, ldy); + } +} + +} // namespace RandBLAS::sparse_data diff --git a/RandBLAS/sparse_data/csc_spsymm_impl.hh b/RandBLAS/sparse_data/csc_spsymm_impl.hh new file mode 100644 index 00000000..f5fef594 --- /dev/null +++ b/RandBLAS/sparse_data/csc_spsymm_impl.hh @@ -0,0 +1,92 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include "RandBLAS/exceptions.hh" +#include "RandBLAS/util.hh" +#include "RandBLAS/sparse_data/base.hh" +#include "RandBLAS/sparse_data/csc_matrix.hh" +#include "RandBLAS/sparse_data/spsymm_internal.hh" +#include + +namespace RandBLAS::sparse_data { + +// ============================================================================= +/// CSC fallback for symmetric sparse-times-dense. +/// +/// Same semantics as the CSR variant; iteration is column-driven via +/// colptr/rowidxs/vals. For CSC stored with uplo=Upper, structurally +/// populated entries satisfy row <= col; for uplo=Lower, row >= col. +/// Entries outside the named triangle are silently skipped. +template +void csc_spsymm( + blas::Layout layout, + blas::Side side, + blas::Uplo uplo, + int64_t m, int64_t n, + T alpha, + const CSCMatrix& A, + const T* B, int64_t ldb, + T beta, + T* Y, int64_t ldy +) { + randblas_require(A.n_rows == A.n_cols); + int64_t k = (side == blas::Side::Left) ? m : n; + randblas_require(A.n_rows == k); + + RandBLAS::util::lascl(layout, m, n, beta, Y, ldy); + if (alpha == T(0)) return; + + if (side == blas::Side::Left) { + // Y = alpha * A * B + ... (A is m-by-m) + for (int64_t j = 0; j < m; ++j) { + for (int64_t p = A.colptr[j]; p < A.colptr[j+1]; ++p) { + sint_t i = A.rowidxs[p]; + if (uplo == blas::Uplo::Upper && i > j) continue; + if (uplo == blas::Uplo::Lower && i < j) continue; + T av = alpha * A.vals[p]; + internal::spsymm_scatter_left(layout, n, av, i, j, B, ldb, Y, ldy); + } + } + } else { + // Y = alpha * B * A + ... (A is n-by-n) + for (int64_t j = 0; j < n; ++j) { + for (int64_t p = A.colptr[j]; p < A.colptr[j+1]; ++p) { + sint_t i = A.rowidxs[p]; + if (uplo == blas::Uplo::Upper && i > j) continue; + if (uplo == blas::Uplo::Lower && i < j) continue; + T av = alpha * A.vals[p]; + internal::spsymm_scatter_right(layout, m, av, i, j, B, ldb, Y, ldy); + } + } + } +} + +} // namespace RandBLAS::sparse_data diff --git a/RandBLAS/sparse_data/csr_spsymm_impl.hh b/RandBLAS/sparse_data/csr_spsymm_impl.hh new file mode 100644 index 00000000..f088da96 --- /dev/null +++ b/RandBLAS/sparse_data/csr_spsymm_impl.hh @@ -0,0 +1,99 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include "RandBLAS/exceptions.hh" +#include "RandBLAS/util.hh" +#include "RandBLAS/sparse_data/base.hh" +#include "RandBLAS/sparse_data/csr_matrix.hh" +#include "RandBLAS/sparse_data/spsymm_internal.hh" +#include + +namespace RandBLAS::sparse_data { + + +// ============================================================================= +/// CSR fallback for the symmetric sparse-times-dense kernel. +/// +/// Computes Y = alpha * op(A) * B + beta * Y (side=Left) or +/// Y = alpha * B * op(A) + beta * Y (side=Right), +/// where A is symmetric and only the triangle named by uplo is structurally +/// stored. The off-diagonal entries contribute twice (the (i,j) and the +/// implied (j,i)); diagonal entries contribute once. Entries that fall +/// outside the named triangle are silently skipped, so the kernel is robust +/// against callers who store both triangles by mistake (it just behaves like +/// the "correctly stored" case). +/// +/// Inner-loop updates use blas::axpy on slices of B and Y. +template +void csr_spsymm( + blas::Layout layout, + blas::Side side, + blas::Uplo uplo, + int64_t m, int64_t n, + T alpha, + const CSRMatrix& A, + const T* B, int64_t ldb, + T beta, + T* Y, int64_t ldy +) { + randblas_require(A.n_rows == A.n_cols); + int64_t k = (side == blas::Side::Left) ? m : n; + randblas_require(A.n_rows == k); + + RandBLAS::util::lascl(layout, m, n, beta, Y, ldy); + if (alpha == T(0)) return; + + if (side == blas::Side::Left) { + // Y = alpha * A * B + ... (A is m-by-m) + for (int64_t i = 0; i < m; ++i) { + for (int64_t p = A.rowptr[i]; p < A.rowptr[i+1]; ++p) { + sint_t j = A.colidxs[p]; + if (uplo == blas::Uplo::Upper && j < i) continue; + if (uplo == blas::Uplo::Lower && j > i) continue; + T av = alpha * A.vals[p]; + internal::spsymm_scatter_left(layout, n, av, i, j, B, ldb, Y, ldy); + } + } + } else { + // Y = alpha * B * A + ... (A is n-by-n) + for (int64_t i = 0; i < n; ++i) { + for (int64_t p = A.rowptr[i]; p < A.rowptr[i+1]; ++p) { + sint_t j = A.colidxs[p]; + if (uplo == blas::Uplo::Upper && j < i) continue; + if (uplo == blas::Uplo::Lower && j > i) continue; + T av = alpha * A.vals[p]; + internal::spsymm_scatter_right(layout, m, av, i, j, B, ldb, Y, ldy); + } + } + } +} + +} // namespace RandBLAS::sparse_data diff --git a/RandBLAS/sparse_data/mkl_spmm_impl.hh b/RandBLAS/sparse_data/mkl_spmm_impl.hh index 6d2af501..1a3e0c27 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" @@ -256,6 +257,38 @@ MKLSparseHandle make_mkl_handle(const SpMat& A) { } } +// ============================================================================ +// Type-dispatched wrapper around mkl_sparse_d_mm / mkl_sparse_s_mm. +// Shared between mkl_left_spmm (general A) and mkl_spsymm (symmetric A); +// the caller controls the matrix_descr (general vs symmetric) and the +// `op` flag, plus the post-call status interpretation. +// ============================================================================ +template +inline sparse_status_t mkl_sparse_mm_call( + sparse_operation_t op, T alpha, sparse_matrix_t A_handle, + const struct matrix_descr& descr, + sparse_layout_t mkl_layout, + const T* B, int64_t n_rhs, int64_t ldb, + T beta, T* C, int64_t ldc +) { + if constexpr (std::is_same_v) { + return mkl_sparse_d_mm( + op, alpha, A_handle, descr, mkl_layout, + B, (MKL_INT)n_rhs, (MKL_INT)ldb, + beta, C, (MKL_INT)ldc + ); + } else if constexpr (std::is_same_v) { + return mkl_sparse_s_mm( + op, alpha, A_handle, descr, mkl_layout, + B, (MKL_INT)n_rhs, (MKL_INT)ldb, + beta, C, (MKL_INT)ldc + ); + } else { + static_assert(sizeof(T) == 0, "MKL sparse BLAS only supports float and double."); + // see GitHub PR #155 for why we don't use static_assert(false, ...). + } +} + // ============================================================================ // MKL-accelerated left_spmm: C = alpha * op(A) * op(B) + beta * C // where A is sparse, B and C are dense. @@ -312,22 +345,11 @@ bool mkl_left_spmm( struct matrix_descr descr; descr.type = SPARSE_MATRIX_TYPE_GENERAL; - sparse_status_t status; - if constexpr (std::is_same_v) { - status = mkl_sparse_d_mm( - to_mkl_op(opA), alpha, h.handle, descr, - to_mkl_layout(layout), - B, (MKL_INT)n, (MKL_INT)ldb, - beta, C, (MKL_INT)ldc - ); - } else if constexpr (std::is_same_v) { - status = mkl_sparse_s_mm( - to_mkl_op(opA), alpha, h.handle, descr, - to_mkl_layout(layout), - B, (MKL_INT)n, (MKL_INT)ldb, - beta, C, (MKL_INT)ldc - ); - } + sparse_status_t status = mkl_sparse_mm_call( + to_mkl_op(opA), alpha, h.handle, descr, + to_mkl_layout(layout), + B, n, ldb, beta, C, ldc + ); check_mkl_status(status, "mkl_sparse_mm"); return true; // signal: MKL handled it } @@ -369,20 +391,8 @@ void mkl_spgemm_to_dense( // 3. C = alpha * temp + C (if needed) if (alpha == (T)0) { - // Just scale C by beta - if (beta == (T)0) { - int64_t total = (layout == blas::Layout::ColMajor) ? ldc * n : ldc * m; - std::fill(C, C + total, (T)0); - } else if (beta != (T)1) { - // Scale each column/row of C - if (layout == blas::Layout::ColMajor) { - for (int64_t j = 0; j < n; ++j) - blas::scal(m, beta, &C[j * ldc], 1); - } else { - for (int64_t i = 0; i < m; ++i) - blas::scal(n, beta, &C[i * ldc], 1); - } - } + // Just scale C by beta. + RandBLAS::util::lascl(layout, m, n, beta, C, ldc); return; } @@ -412,17 +422,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..9848f73c --- /dev/null +++ b/RandBLAS/sparse_data/mkl_spsymm_impl.hh @@ -0,0 +1,150 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include "RandBLAS/config.h" + +#if defined(RandBLAS_HAS_MKL) + +#if defined(BLAS_ILP64) && !defined(MKL_ILP64) +#define MKL_ILP64 +#endif + +#include +#include + +#include "RandBLAS/exceptions.hh" +#include "RandBLAS/sparse_data/base.hh" +#include "RandBLAS/sparse_data/coo_matrix.hh" +#include "RandBLAS/sparse_data/csr_matrix.hh" +#include "RandBLAS/sparse_data/csc_matrix.hh" +#include "RandBLAS/sparse_data/mkl_spmm_impl.hh" // reuse make_mkl_handle, to_mkl_layout, check_mkl_status + +namespace RandBLAS::sparse_data::mkl { + +// ============================================================================ +// MKL-accelerated symmetric SpMM: Y = alpha * A * B + beta * Y (side=Left only). +// A is symmetric sparse (one triangle stored, named by uplo); B and Y dense. +// +// Returns true if MKL handled the call, false to signal fallback to the +// hand-rolled per-format kernel. MKL applies alpha and beta internally; the +// caller therefore passes them through unchanged. +// +// Known limitation that triggers a fallback: +// - Index type mismatched with MKL_INT. +// +// Both Side values are handled: +// - side=Left (Y = alpha*A*B + beta*Y): direct call to mkl_sparse_d_mm +// with the symmetric descriptor. +// - side=Right (Y = alpha*B*A + beta*Y): A is symmetric so A == A^T, +// and (B*A)^T = A^T * B^T = A * B^T. We tell MKL to compute A*B^T +// into Y^T by flipping the MKL layout flag. Concretely, the same +// user buffers for B and Y are reinterpreted in the opposite layout +// (ColMajor <-> RowMajor); ldb and ldy carry through unchanged +// because the reinterpretation has the same leading-dim semantics +// (rows-of-RowMajor have the same stride as cols-of-ColMajor). +// +// CSC handling: MKL's mkl_sparse_d_mm returns NOT_SUPPORTED for CSC even +// with a symmetric descriptor. We work around this by taking the +// CSC.transpose() view (a lightweight CSR view over the same buffers, +// since A is symmetric so A == A^T) and recursing. The triangle the user +// named in the CSC is in the *opposite* triangle of the CSR view (a CSC +// Upper entry at (i, j) with i <= j becomes a CSR-view entry at (j, i) +// with j >= i, i.e., Lower in the CSR view), so the recursive call flips +// uplo. +// ============================================================================ +template +bool mkl_spsymm( + blas::Layout layout, + blas::Side side, + blas::Uplo uplo, + int64_t m, int64_t n, + T alpha, + const SpMat &A, + const T *B, + int64_t ldb, + T beta, + T *Y, + int64_t ldy +) { + using sint_t = typename SpMat::index_t; + constexpr bool is_csc = std::is_same_v>; + + if constexpr (is_csc) { + // Symmetric A: A == A^T. The CSC->CSR view is lightweight (same + // buffers, reinterpreted) and gives us a CSR matrix MKL accepts; + // the structurally stored triangle moves from {Upper,Lower} to + // the opposite side in the CSR view. + auto At = A.transpose(); + blas::Uplo uplo_flipped = (uplo == blas::Uplo::Upper) + ? blas::Uplo::Lower + : blas::Uplo::Upper; + return mkl_spsymm(layout, side, uplo_flipped, m, n, + alpha, At, B, ldb, beta, Y, ldy); + } + + auto h = make_mkl_handle(A); + + 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; + + // For side=Right, reinterpret B and Y in the opposite layout to + // present them to MKL as B^T and Y^T; MKL then computes Y^T = A*B^T + // which is the transpose of Y = B*A. The number of MKL-side + // right-hand-side columns is m (the user's row count) rather than n + // (the user's col count, = A's order). + blas::Layout mkl_layout = (side == blas::Side::Left) + ? layout + : (layout == blas::Layout::ColMajor ? blas::Layout::RowMajor + : blas::Layout::ColMajor); + int64_t n_rhs = (side == blas::Side::Left) ? n : m; + + sparse_status_t status = mkl_sparse_mm_call( + SPARSE_OPERATION_NON_TRANSPOSE, alpha, h.handle, descr, + to_mkl_layout(mkl_layout), + B, n_rhs, ldb, beta, Y, ldy + ); + + // Some MKL versions return NOT_SUPPORTED for combinations we couldn't + // predict. Don't throw -- signal fallback. + if (status == SPARSE_STATUS_NOT_SUPPORTED) + return false; + + check_mkl_status(status, "mkl_sparse_mm (symmetric)"); + return true; +} + +} // namespace RandBLAS::sparse_data::mkl + +#endif // RandBLAS_HAS_MKL diff --git a/RandBLAS/sparse_data/spmm_dispatch.hh b/RandBLAS/sparse_data/spmm_dispatch.hh index c2edf112..192b6aee 100644 --- a/RandBLAS/sparse_data/spmm_dispatch.hh +++ b/RandBLAS/sparse_data/spmm_dispatch.hh @@ -31,6 +31,7 @@ #include "RandBLAS/base.hh" #include "RandBLAS/exceptions.hh" +#include "RandBLAS/util.hh" #include "RandBLAS/sparse_data/base.hh" #include "RandBLAS/sparse_data/coo_matrix.hh" #include "RandBLAS/sparse_data/csr_matrix.hh" @@ -125,14 +126,11 @@ void left_spmm( if (layout == Layout::ColMajor) { randblas_require(ldb >= rows_B); randblas_require(ldc >= d); - for (int64_t i = 0; i < n; ++i) - RandBLAS::util::safe_scal(d, beta, &C[i*ldc]); } else { randblas_require(ldc >= n); randblas_require(ldb >= cols_B); - for (int64_t i = 0; i < d; ++i) - RandBLAS::util::safe_scal(n, beta, &C[i*ldc]); } + RandBLAS::util::lascl(layout, d, n, beta, C, ldc); if (alpha == (T) 0) return; diff --git a/RandBLAS/sparse_data/spsymm_dispatch.hh b/RandBLAS/sparse_data/spsymm_dispatch.hh new file mode 100644 index 00000000..20354f85 --- /dev/null +++ b/RandBLAS/sparse_data/spsymm_dispatch.hh @@ -0,0 +1,243 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include "RandBLAS/base.hh" +#include "RandBLAS/exceptions.hh" +#include "RandBLAS/sparse_data/base.hh" +#include "RandBLAS/sparse_data/coo_matrix.hh" +#include "RandBLAS/sparse_data/csr_matrix.hh" +#include "RandBLAS/sparse_data/csc_matrix.hh" +#include "RandBLAS/sparse_data/coo_spsymm_impl.hh" +#include "RandBLAS/sparse_data/csr_spsymm_impl.hh" +#include "RandBLAS/sparse_data/csc_spsymm_impl.hh" +#include "RandBLAS/sparse_data/symmetric.hh" +#include "RandBLAS/config.h" + +#if defined(RandBLAS_HAS_MKL) +#include "RandBLAS/sparse_data/mkl_spsymm_impl.hh" +#endif + +#include + + +namespace RandBLAS::sparse_data { + +// ============================================================================= +/// Dispatched symmetric sparse-times-dense kernel (Case C in +/// project-plans/randblas-symm-plan.md). +/// +/// @verbatim embed:rst:leading-slashes +/// Computes +/// +/// .. math:: +/// \mat(Y) = \alpha \cdot \op_{\ttt{side}}(\mat(A), \mat(B)) + \beta \cdot \mat(Y), +/// +/// where: +/// +/// - :math:`\mat(A)` is a sparse symmetric matrix stored in COO, CSR, or +/// CSC format. Only the triangle named by :math:`\ttt{uplo}` is read. +/// The opposite triangle is implied by symmetry. The matrix is required +/// to be square. +/// - :math:`\mat(B)` and :math:`\mat(Y)` are dense, both m-by-n. +/// - For :math:`\ttt{side} = \ttt{Left}`, :math:`\mat(A)` is m-by-m and the +/// operation is :math:`\mat(Y) = \alpha A B + \beta Y`. +/// - For :math:`\ttt{side} = \ttt{Right}`, :math:`\mat(A)` is n-by-n and +/// the operation is :math:`\mat(Y) = \alpha B A + \beta Y`. +/// +/// Dispatch: +/// - If RandBLAS was built with MKL and the index type matches MKL_INT, +/// try the MKL fast path (mkl_sparse_d_mm with +/// ``SPARSE_MATRIX_TYPE_SYMMETRIC`` descriptor). The MKL path covers +/// side=Left for CSR and COO; it returns false (fall back) for CSC and +/// for side=Right. +/// - Otherwise (or on MKL fall back), call the format-specific fallback +/// kernel (``csr_spsymm`` / ``csc_spsymm`` / ``coo_spsymm``). +/// @endverbatim +template +void spsymm( + blas::Layout layout, + blas::Side side, + blas::Uplo uplo, + int64_t m, int64_t n, + T alpha, + const SpMat& A, + const T* B, int64_t ldb, + T beta, + T* Y, int64_t ldy +) { + using sint_t = typename SpMat::index_t; + constexpr bool is_coo = std::is_same_v>; + constexpr bool is_csr = std::is_same_v>; + constexpr bool is_csc = std::is_same_v>; + static_assert(is_coo || is_csr || is_csc, + "RandBLAS::sparse_data::spsymm requires COO, CSR, or CSC."); + + randblas_require(A.n_rows == A.n_cols); + int64_t k = (side == blas::Side::Left) ? m : n; + randblas_require(A.n_rows == k); + +#if defined(RandBLAS_HAS_MKL) + if constexpr (sizeof(sint_t) == sizeof(MKL_INT)) { + bool handled = mkl::mkl_spsymm( + layout, side, uplo, m, n, + alpha, A, B, ldb, beta, Y, ldy + ); + if (handled) return; + } +#endif + + // Fallback path + if constexpr (is_csr) { + csr_spsymm(layout, side, uplo, m, n, alpha, A, B, ldb, beta, Y, ldy); + } else if constexpr (is_csc) { + csc_spsymm(layout, side, uplo, m, n, alpha, A, B, ldb, beta, Y, ldy); + } else if constexpr (is_coo) { + coo_spsymm(layout, side, uplo, m, n, alpha, A, B, ldb, beta, Y, ldy); + } +} + +// ============================================================================= +/// Case D: sparse-symmetric A times sparse B, dense output. +/// +/// @verbatim embed:rst:leading-slashes +/// Implemented as a composition: densify ``B`` into a temporary dense buffer +/// matching the caller's ``layout``, then call the Case-C ``spsymm`` overload +/// (sparse-symm × dense) on the result. This works in any build (the MKL +/// fast path or the per-format hand kernel both apply, depending on the +/// build), and across all 3 × 3 = 9 sparse-format pairings for ``(A, B)``, +/// since the densification picks the right format-specific helper. +/// +/// Why this composition rather than a single MKL ``sp2m`` call: MKL's +/// ``mkl_sparse_sp2m`` returns ``SPARSE_STATUS_NOT_SUPPORTED`` when the +/// ``matrix_descr`` on either operand is ``SPARSE_MATRIX_TYPE_SYMMETRIC``; +/// only ``GENERAL`` is supported there. ``mkl_sparse_d_spmmd`` (which writes +/// directly to dense ``C``) accepts no descriptor at all. So the symmetric +/// expansion has to happen on the RandBLAS side. Composing through Case C +/// gets it for free at the cost of a temporary dense buffer for ``B`` +/// (cost: ``O(m * n)``), and for the typical RandNLA workload where ``B`` +/// is a sketching operator with ``nnz(B) << m*n`` the cost is small. +/// +/// The dimensions of the densified ``B`` are the same as the user's ``Y``: +/// ``m``-by-``n``, with the user's ``layout`` and a tight leading dim. The +/// densified buffer is freed when the temporary ``std::vector`` goes out +/// of scope. ``Y`` itself is never touched until the Case-C call. +/// @endverbatim +template +void spsymm( + blas::Layout layout, + blas::Side side, + blas::Uplo uplo, + int64_t m, int64_t n, + T alpha, + const SpMatA& A, + const SpMatB& B, + T beta, + T* Y, int64_t ldy +) { + static_assert(std::is_same_v, + "Case D: A and B must share scalar_t."); + + randblas_require(A.n_rows == A.n_cols); + int64_t k = (side == blas::Side::Left) ? m : n; + randblas_require(A.n_rows == k); + randblas_require(B.n_rows == m); + randblas_require(B.n_cols == n); + + // Densify B into a tight buffer in the caller's layout. + int64_t ldb_dense = (layout == blas::Layout::ColMajor) ? m : n; + std::vector B_dense(static_cast(m) * static_cast(n), T(0)); + + using sint_B = typename SpMatB::index_t; + if constexpr (std::is_same_v>) { + coo::coo_to_dense(B, layout, B_dense.data()); + } else if constexpr (std::is_same_v>) { + csr::csr_to_dense(B, layout, B_dense.data()); + } else if constexpr (std::is_same_v>) { + csc::csc_to_dense(B, layout, B_dense.data()); + } else { + static_assert(sizeof(SpMatB) == 0, + "RandBLAS::sparse_data::spsymm: SpMatB must be COO, CSR, or CSC."); + } + + // Compose into Case C with the densified B. + spsymm(layout, side, uplo, m, n, + alpha, A, B_dense.data(), ldb_dense, beta, Y, ldy); +} + +} // end namespace RandBLAS::sparse_data + + +namespace RandBLAS { + +// ============================================================================= +/// Convenience wrapper for symmetric sparse-times-dense matmul. +/// +/// Computes :math:`\ttt{Y} = \alpha A B + \beta Y` (side=Left default), where +/// :math:`A` is the symmetric sparse matrix and only the triangle named by +/// :math:`\ttt{uplo}` is read. +template +inline void spsymm( + blas::Layout layout, + blas::Uplo uplo, + int64_t m, int64_t n, + T alpha, + const SpMat& A, + const T* B, int64_t ldb, + T beta, + T* Y, int64_t ldy +) { + RandBLAS::sparse_data::spsymm( + layout, blas::Side::Left, uplo, m, n, + alpha, A, B, ldb, beta, Y, ldy + ); +} + +// ============================================================================= +/// Convenience overload taking a Symmetric wrapper. The uplo is read +/// from the wrapper. Defaults to side=Left. +template +inline void spsymm( + blas::Layout layout, + int64_t m, int64_t n, + T alpha, + const Symmetric& A_sym, + const T* B, int64_t ldb, + T beta, + T* Y, int64_t ldy +) { + RandBLAS::sparse_data::spsymm( + layout, blas::Side::Left, A_sym.uplo, m, n, + alpha, A_sym.A, B, ldb, beta, Y, ldy + ); +} + +} // end namespace RandBLAS diff --git a/RandBLAS/sparse_data/spsymm_internal.hh b/RandBLAS/sparse_data/spsymm_internal.hh new file mode 100644 index 00000000..4aea5678 --- /dev/null +++ b/RandBLAS/sparse_data/spsymm_internal.hh @@ -0,0 +1,87 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include +#include + +namespace RandBLAS::sparse_data::internal { + +// ============================================================================= +/// Add the contribution of one stored A(i, j) = v entry (and the implied +/// symmetric A(j, i) = v when i != j) to a side=Left spsymm output: +/// +/// Y[i, :] += av * B[j, :] +/// if (i != j) Y[j, :] += av * B[i, :] +/// +/// where av = alpha * v has already been folded by the caller. The inner +/// row updates use blas::axpy with strides determined by `layout`. +/// +/// Shared by all three spsymm fallback kernels (CSR / CSC / COO) -- they +/// differ only in their outer iteration order; the per-stored-entry scatter +/// is identical. +template +inline void spsymm_scatter_left(blas::Layout layout, int64_t n, T av, + int64_t i, int64_t j, + const T* B, int64_t ldb, + T* Y, int64_t ldy) { + if (layout == blas::Layout::ColMajor) { + blas::axpy(n, av, &B[j], ldb, &Y[i], ldy); + if (i != j) + blas::axpy(n, av, &B[i], ldb, &Y[j], ldy); + } else { + blas::axpy(n, av, &B[j*ldb], 1, &Y[i*ldy], 1); + if (i != j) + blas::axpy(n, av, &B[i*ldb], 1, &Y[j*ldy], 1); + } +} + +// ============================================================================= +/// Side=Right counterpart of spsymm_scatter_left. For stored A(i, j) = v: +/// +/// Y[:, j] += av * B[:, i] +/// if (i != j) Y[:, i] += av * B[:, j] +template +inline void spsymm_scatter_right(blas::Layout layout, int64_t m, T av, + int64_t i, int64_t j, + const T* B, int64_t ldb, + T* Y, int64_t ldy) { + if (layout == blas::Layout::ColMajor) { + blas::axpy(m, av, &B[i*ldb], 1, &Y[j*ldy], 1); + if (i != j) + blas::axpy(m, av, &B[j*ldb], 1, &Y[i*ldy], 1); + } else { + blas::axpy(m, av, &B[i], ldb, &Y[j], ldy); + if (i != j) + blas::axpy(m, av, &B[j], ldb, &Y[i], ldy); + } +} + +} // namespace RandBLAS::sparse_data::internal diff --git a/RandBLAS/sparse_data/symmetric.hh b/RandBLAS/sparse_data/symmetric.hh new file mode 100644 index 00000000..dee1a52a --- /dev/null +++ b/RandBLAS/sparse_data/symmetric.hh @@ -0,0 +1,97 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#pragma once + +#include "RandBLAS/sparse_data/base.hh" +#include "RandBLAS/exceptions.hh" + + +namespace RandBLAS::sparse_data { + +// ============================================================================= +/// Lightweight non-owning wrapper marking a SparseMatrix as symmetric with a +/// stored triangle. +/// +/// @verbatim embed:rst:leading-slashes +/// Holds: +/// - A const reference to the underlying SparseMatrix :math:`A.` +/// - :math:`\ttt{blas::Uplo uplo}`: names the triangle of :math:`A` that is +/// structurally populated. The opposite triangle is implied by symmetry. +/// +/// The wrapper performs **no validation** of the matrix contents --- it is the +/// caller's responsibility to guarantee that the named triangle is correctly +/// populated and that the opposite triangle is either structurally absent or +/// will be ignored by the SYMM-aware consumer (kernels in this library that +/// accept a ``Symmetric`` agree to consult ``uplo`` and respect this +/// contract). Construction does enforce that :math:`A` is square +/// (``A.n_rows == A.n_cols``) via :math:`\ttt{randblas\_require}`. +/// +/// This wrapper is the carrier for symmetric sparse matrices into the +/// :math:`\ttt{spsymm}`-family kernels (project-plans/randblas-symm-plan.md +/// Case C and beyond). It is intentionally separate from the underlying +/// :math:`\ttt{SparseMatrix}` concept so that calling +/// :math:`\ttt{spmm}` / :math:`\ttt{spgemm}` with a ``Symmetric`` +/// argument fails to compile rather than silently treating the matrix as +/// general --- a type-system guard against accidental loss of symmetry. +/// +/// The wrapper is non-owning: it holds a reference to :math:`A`, not a copy. +/// The caller must keep :math:`A` alive for the lifetime of the wrapper. +/// @endverbatim +template +struct Symmetric { + const SpMat& A; + const blas::Uplo uplo; + + using scalar_t = typename SpMat::scalar_t; + using index_t = typename SpMat::index_t; + + Symmetric(const SpMat& A_in, blas::Uplo uplo_in) : A(A_in), uplo(uplo_in) { + randblas_require(A_in.n_rows == A_in.n_cols); + } +}; + + +// ============================================================================= +/// Construct a ``Symmetric`` wrapper around ``A`` with the named triangle. +/// Syntactic sugar for the constructor; lets callers write +/// ``as_symmetric(A, blas::Uplo::Upper)`` and rely on template argument +/// deduction. +template +inline Symmetric as_symmetric(const SpMat& A, blas::Uplo uplo) { + return Symmetric(A, uplo); +} + +} // end namespace RandBLAS::sparse_data + + +namespace RandBLAS { + using RandBLAS::sparse_data::Symmetric; + using RandBLAS::sparse_data::as_symmetric; +} diff --git a/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 597b1e55..3834018e 100644 --- a/RandBLAS/util.hh +++ b/RandBLAS/util.hh @@ -61,6 +61,46 @@ void safe_scal(int64_t n, T a, T* x) { } } +// ============================================================================= +/// \fn lascl(blas::Layout layout, int64_t m, int64_t n, T alpha, T* A, int64_t lda) +/// @verbatim embed:rst:leading-slashes +/// In-place scale a dense :math:`m \times n` matrix: +/// +/// .. math:: +/// A \leftarrow \alpha \cdot A. +/// +/// Named after LAPACK's ``?lascl`` for the convention; the signature is +/// simpler than LAPACK's (no ``CFROM`` / ``CTO`` overflow protection, no +/// matrix-type flag --- general dense only). +/// +/// Fast paths: +/// +/// - :math:`\alpha = 1`: returns immediately. +/// - :math:`\alpha = 0`: ``std::fill``-based zero out. +/// - otherwise: per-column (ColMajor) or per-row (RowMajor) ``blas::scal``. +/// @endverbatim +template +inline void lascl(blas::Layout layout, int64_t m, int64_t n, T alpha, T* A, int64_t lda) { + if (alpha == T(1)) return; + if (alpha == T(0)) { + if (layout == blas::Layout::ColMajor) { + for (int64_t j = 0; j < n; ++j) + std::fill(A + j * lda, A + j * lda + m, T(0)); + } else { + for (int64_t i = 0; i < m; ++i) + std::fill(A + i * lda, A + i * lda + n, T(0)); + } + return; + } + if (layout == blas::Layout::ColMajor) { + for (int64_t j = 0; j < n; ++j) + blas::scal(m, alpha, A + j * lda, 1); + } else { + for (int64_t i = 0; i < m; ++i) + blas::scal(n, alpha, A + i * lda, 1); + } +} + template void omatcopy(int64_t m, int64_t n, const T* A, int64_t irs_a, int64_t ics_a, T* B, int64_t irs_b, int64_t ics_b) { // TODO: @@ -122,7 +162,11 @@ void flip_layout(blas::Layout layout_in, int64_t m, int64_t n, std::vector &A /// for all :math:`i,j \in \\{0,\ldots,n-1\\}.` An error is raised if any such check fails. /// This function returns immediately without performing any checks if :math:`\ttt{tol} < 0.` /// @endverbatim -/// sketch_symmetric calls this function with \math{\ttt{tol} = 0} by default. +/// (Historical note: pre-Phase-1 of project-plans/randblas-symm-plan.md, +/// sketch_symmetric called this function with \math{\ttt{tol} = 0} by default. +/// With the SYMM-based dispatch, only the triangle named by \math{\ttt{uplo}} +/// is read, so the symmetry check is no longer invoked from sketch_symmetric. +/// require_symmetric remains available as a standalone validator.) /// template void require_symmetric(blas::Layout layout, const T* A, int64_t n, int64_t lda, T tol) { diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index a953bb35..4a3b4938 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -103,5 +103,15 @@ target_link_libraries( spmm_performance PUBLIC RandBLAS blaspp lapackpp ) +add_executable( + spsymm_performance simple-kernel-benchmarks/spsymm_performance.cc +) +target_include_directories( + spsymm_performance PUBLIC ${Random123_DIR} +) +target_link_libraries( + spsymm_performance PUBLIC RandBLAS blaspp lapackpp +) + diff --git a/examples/simple-kernel-benchmarks/spmm_performance.cc b/examples/simple-kernel-benchmarks/spmm_performance.cc index 03bb7c43..411fb055 100644 --- a/examples/simple-kernel-benchmarks/spmm_performance.cc +++ b/examples/simple-kernel-benchmarks/spmm_performance.cc @@ -88,13 +88,7 @@ void handrolled_left_spmm_csr( const T *B, int64_t ldb, T beta, T *C, int64_t ldc ) { // Apply beta to C (same as dispatch) - if (layout == Layout::ColMajor) { - for (int64_t i = 0; i < n; ++i) - RandBLAS::util::safe_scal(d, beta, &C[i*ldc]); - } else { - for (int64_t i = 0; i < d; ++i) - RandBLAS::util::safe_scal(n, beta, &C[i*ldc]); - } + RandBLAS::util::lascl(layout, d, n, beta, C, ldc); if (alpha == (T)0) return; // Call the hand-rolled CSR kernel directly (same as dispatch fallback) diff --git a/examples/simple-kernel-benchmarks/spsymm_performance.cc b/examples/simple-kernel-benchmarks/spsymm_performance.cc new file mode 100644 index 00000000..bdebf44a --- /dev/null +++ b/examples/simple-kernel-benchmarks/spsymm_performance.cc @@ -0,0 +1,296 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// ============================================================================ +// SPSYMM / SKETCH_SYMMETRIC PERFORMANCE BENCHMARK +// ============================================================================ +// +// This benchmark answers two performance questions about the symm-kernels +// added in https://github.com/BallisticLA/RandBLAS/pull/163 : +// +// 1. Sparse: How much does the new RandBLAS::spsymm (one-triangle storage, +// MKL fast path via SPARSE_MATRIX_TYPE_SYMMETRIC) speed up over the +// "both-triangles + right_spmm" workaround that downstream code had to +// use before this PR? +// +// 2. Dense: How much does the rewritten RandBLAS::sketch_symmetric (now +// backed by blas::symm) speed up over its pre-PR behaviour, which +// silently forwarded to sketch_general -> lskge3 -> blas::gemm with no +// symmetry exploitation? +// +// NOTATION: +// A_symm - symmetric matrix of order n_A (dense or sparse, one triangle) +// B - dense matrix (n_A x d for side=Left; d x n_A for side=Right) +// Y - dense result matrix (same shape as B) +// density - fraction of upper-triangle entries of A_symm that are nonzero +// (the implied lower-triangle entries follow by symmetry) +// +// USAGE: +// ./spsymm_performance # default sweep +// ./spsymm_performance n_A d density [num_trials] # single config +// +// Defaults: num_trials=10, density=0.05. +// +// EXAMPLES: +// env OMP_NUM_THREADS=8 ./spsymm_performance +// env OMP_NUM_THREADS=8 ./spsymm_performance 2000 200 0.01 +// +// ============================================================================ + +#include + +#include +#include +#include +#include +#include +#include +#include + +using std::chrono::steady_clock; +using std::chrono::duration_cast; +using std::chrono::microseconds; +using blas::Layout; +using blas::Uplo; +using blas::Side; +using blas::Op; + +using T = double; +using sint_t = int64_t; +using SpMat = RandBLAS::sparse_data::CSRMatrix; + + +// ============================================================================ +// Helpers: random dense symmetric matrix + sparse counterparts. +// ============================================================================ + +void make_dense_symmetric(int64_t n, std::vector& A, uint64_t seed) { + A.assign(static_cast(n) * n, T(0)); + std::mt19937_64 rng(seed); + std::uniform_real_distribution uni(-1.0, 1.0); + for (int64_t j = 0; j < n; ++j) { + for (int64_t i = 0; i <= j; ++i) { + T v = uni(rng); + A[i + j * n] = v; + if (i != j) A[j + i * n] = v; // mirror + } + } +} + +// Sparsify the upper triangle of A (including diagonal) at the given density, +// then mirror into the lower triangle. +void sparsify_upper_then_mirror(int64_t n, std::vector& A, double density, uint64_t seed) { + std::mt19937_64 rng(seed); + std::uniform_real_distribution uni01(0.0, 1.0); + for (int64_t j = 0; j < n; ++j) { + for (int64_t i = 0; i <= j; ++i) { + if (uni01(rng) >= density) { + A[i + j * n] = T(0); + if (i != j) A[j + i * n] = T(0); + } + } + } +} + +// Build a one-triangle (Upper) CSR view by walking the upper triangle of A_dense. +SpMat build_csr_upper_only(int64_t n, const std::vector& A_dense) { + std::vector rowptr(n + 1, 0); + for (int64_t i = 0; i < n; ++i) { + for (int64_t j = i; j < n; ++j) + if (A_dense[i + j * n] != T(0)) ++rowptr[i + 1]; + } + for (int64_t i = 0; i < n; ++i) rowptr[i + 1] += rowptr[i]; + int64_t nnz = rowptr[n]; + + SpMat A_sparse(n, n); + A_sparse.reserve(nnz); + std::vector tmp_rp(rowptr); // running cursor per row + for (int64_t i = 0; i < n; ++i) { + for (int64_t j = i; j < n; ++j) { + T v = A_dense[i + j * n]; + if (v != T(0)) { + int64_t pos = tmp_rp[i]++; + A_sparse.colidxs[pos] = static_cast(j); + A_sparse.vals[pos] = v; + } + } + } + for (int64_t i = 0; i <= n; ++i) A_sparse.rowptr[i] = rowptr[i]; + return A_sparse; +} + +// Build a full (both-triangles-stored) CSR for the pre-PR workaround comparison. +SpMat build_csr_both_triangles(int64_t n, const std::vector& A_dense) { + std::vector rowptr(n + 1, 0); + for (int64_t i = 0; i < n; ++i) + for (int64_t j = 0; j < n; ++j) + if (A_dense[i + j * n] != T(0)) ++rowptr[i + 1]; + for (int64_t i = 0; i < n; ++i) rowptr[i + 1] += rowptr[i]; + int64_t nnz = rowptr[n]; + + SpMat A_sparse(n, n); + A_sparse.reserve(nnz); + std::vector tmp_rp(rowptr); + for (int64_t i = 0; i < n; ++i) { + for (int64_t j = 0; j < n; ++j) { + T v = A_dense[i + j * n]; + if (v != T(0)) { + int64_t pos = tmp_rp[i]++; + A_sparse.colidxs[pos] = static_cast(j); + A_sparse.vals[pos] = v; + } + } + } + for (int64_t i = 0; i <= n; ++i) A_sparse.rowptr[i] = rowptr[i]; + return A_sparse; +} + + +// ============================================================================ +// Timing helper: median + min over num_trials. +// ============================================================================ +template +std::pair run_trials(Func&& func, int num_trials) { + std::vector times; + times.reserve(num_trials); + for (int t = 0; t < num_trials; ++t) { + auto start = steady_clock::now(); + func(); + auto end = steady_clock::now(); + times.push_back(duration_cast(end - start).count()); + } + std::sort(times.begin(), times.end()); + return {times[0], times[num_trials / 2]}; +} + +void print_row(const std::string& name, long min_us, long med_us, long baseline) { + double ratio = (baseline > 0) ? (double)min_us / baseline : 1.0; + std::cout << " " << std::setw(38) << std::left << name + << std::setw(10) << std::right << min_us + << std::setw(10) << med_us + << std::setw(10) << std::fixed << std::setprecision(2) << ratio << "x\n"; +} + + +// ============================================================================ +// run_config: one (n_A, d, density) point. Compares the symm-aware path +// against the workaround / pre-PR baselines. +// ============================================================================ +void run_config(int64_t n_A, int64_t d, double density, int num_trials) { + uint64_t seed = 12345; + Layout layout = Layout::ColMajor; + T alpha = T(1.0), beta = T(0.0); + + std::cout << "--- A is " << n_A << "x" << n_A + << " symmetric, B is " << n_A << "x" << d + << ", density=" << std::setprecision(4) << density + << " (median + min over " << num_trials << " trials) ---\n"; + + std::vector A_full(static_cast(n_A) * n_A); + make_dense_symmetric(n_A, A_full, seed); + sparsify_upper_then_mirror(n_A, A_full, density, seed + 1); + + SpMat A_csr_upper = build_csr_upper_only(n_A, A_full); + SpMat A_csr_full = build_csr_both_triangles(n_A, A_full); + + int64_t actual_nnz_upper = A_csr_upper.rowptr[n_A]; + int64_t actual_nnz_full = A_csr_full.rowptr[n_A]; + std::cout << " upper-triangle nnz = " << actual_nnz_upper + << " (full = " << actual_nnz_full << ")\n"; + + std::vector B(static_cast(n_A) * d); + { + std::mt19937_64 rng(seed + 2); + std::uniform_real_distribution uni(-1.0, 1.0); + for (auto& x : B) x = uni(rng); + } + std::vector Y(static_cast(n_A) * d, T(0)); + + std::cout << "\n SPARSE (side=Left, Y = A*B):\n"; + std::cout << " " << std::setw(38) << std::left << "kernel" + << std::setw(10) << std::right << "min(us)" + << std::setw(10) << "med(us)" + << std::setw(10) << "ratio\n"; + + // 1. RandBLAS::spsymm on one-triangle CSR (PR #163's new path). + auto [t1_min, t1_med] = run_trials([&] { + RandBLAS::spsymm(layout, Uplo::Upper, n_A, d, + alpha, A_csr_upper, B.data(), n_A, + beta, Y.data(), n_A); + }, num_trials); + + // 2. RandBLAS::spmm with full storage (the pre-PR workaround: + // both triangles materialized into a general-sparse CSR). + auto [t2_min, t2_med] = run_trials([&] { + RandBLAS::spmm(layout, Op::NoTrans, Op::NoTrans, + n_A, d, n_A, + alpha, A_csr_full, B.data(), n_A, + beta, Y.data(), n_A); + }, num_trials); + + // 3. Dense blas::symm on the fully populated dense A (the "ideal" + // BLAS-symm baseline -- doesn't exploit sparsity, but is the + // fastest possible dense SYMM). + auto [t3_min, t3_med] = run_trials([&] { + blas::symm(layout, Side::Left, Uplo::Upper, n_A, d, + alpha, A_full.data(), n_A, B.data(), n_A, + beta, Y.data(), n_A); + }, num_trials); + + print_row("RandBLAS::spsymm (one tri + MKL)", t1_min, t1_med, t1_min); + print_row("RandBLAS::spmm (full + MKL)", t2_min, t2_med, t1_min); + print_row("blas::symm (dense ref)", t3_min, t3_med, t1_min); + + // Dense-symmetric sketching comparison: new SYMM-backed vs the + // pre-PR GEMM-forwarding behaviour. + std::cout << "\n DENSE (sketch_symmetric vs sketch_general on dense-symm A):\n"; + std::cout << " " << std::setw(38) << std::left << "kernel" + << std::setw(10) << std::right << "min(us)" + << std::setw(10) << "med(us)" + << std::setw(10) << "ratio\n"; + + RandBLAS::DenseDist DS(n_A, d, RandBLAS::ScalarDist::Uniform); + RandBLAS::DenseSkOp S(DS, static_cast(seed + 3)); + RandBLAS::fill_dense(S); + + // 1. New: RandBLAS::sketch_symmetric (SYMM-backed via blas::symm). + auto [s1_min, s1_med] = run_trials([&] { + RandBLAS::sketch_symmetric(layout, Uplo::Upper, n_A, d, + alpha, A_full.data(), n_A, S, 0, 0, + beta, Y.data(), n_A); + }, num_trials); + + // 2. Old: equivalent call via sketch_general (the pre-PR behaviour + // when sketch_symmetric silently forwarded to GEMM). + auto [s2_min, s2_med] = run_trials([&] { + RandBLAS::sketch_general(layout, Op::NoTrans, Op::NoTrans, + n_A, d, n_A, + alpha, A_full.data(), n_A, S, 0, 0, + beta, Y.data(), n_A); + }, num_trials); + + print_row("sketch_symmetric (new SYMM path)", s1_min, s1_med, s1_min); + print_row("sketch_general (pre-PR GEMM path)", s2_min, s2_med, s1_min); + + std::cout << "\n"; +} + + +// ============================================================================ +// main: default 4-point sweep, or single config from argv. +// ============================================================================ +int main(int argc, char** argv) { + if (argc >= 4) { + int64_t n_A = std::stoll(argv[1]); + int64_t d = std::stoll(argv[2]); + double density = std::stod(argv[3]); + int num_trials = (argc >= 5) ? std::stoi(argv[4]) : 10; + run_config(n_A, d, density, num_trials); + } else { + int num_trials = 10; + std::cout << "Default sweep (n_A in {500, 1000, 2000}, d=200, density=0.05)\n\n"; + for (int64_t n_A : {500, 1000, 2000}) { + run_config(n_A, /*d=*/200, /*density=*/0.05, num_trials); + } + } + return 0; +} diff --git a/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc b/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc index 7c4f23af..49a56cbf 100644 --- a/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc +++ b/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc @@ -341,9 +341,10 @@ void sketch_to_tqrcp(SpMat &A, int64_t k, T* Q, int64_t ldq, T* Y, int64_t ldy, lapack::geqp3(k, n, Y, ldy, piv, tau), "GEQP3 : ") // ================================================================ - // Step 2: copy A(:, piv(0)-1), ..., A(:, piv(k)-1) into dense Q + // Step 2: copy A(:, piv(0)-1), ..., A(:, piv(k)-1) into dense Q. + // Zero the m-by-k destination block in one pass before scattering. + RandBLAS::util::lascl(blas::Layout::ColMajor, m, k, 0.0, Q, ldq); for (int64_t j = 0; j < k; ++j) { - RandBLAS::util::safe_scal(m, 0.0, Q + j*ldq); for (int64_t ell = A.colptr[piv[j]-1]; ell < A.colptr[piv[j]]; ++ell) { int64_t i = A.rowidxs[ell]; Q[i + ldq*j] = A.vals[ell]; diff --git a/rtd/source/FAQ.rst b/rtd/source/FAQ.rst index ea9a0e1c..f4bed3bb 100644 --- a/rtd/source/FAQ.rst +++ b/rtd/source/FAQ.rst @@ -110,10 +110,16 @@ No support for negative values of "incx" or "incy" in sketch_vector. to extend our SPMV kernels to support that, then we'd happily accept such a contribution. (It shouldn't be hard! We just haven't gotten around to this.) -Symmetric matrices have to be stored as general matrices. - This stems partly from a desire for sketch_symmetric work equally well with DenseSkOp and SparseSkOp. - Another reason is that BLAS' SYMM function doesn't allow transposes, which is a key tool we use - in sketch_general to resolve layout discrepancies between the various arguments. +sketch_symmetric supports both DenseSkOp and SparseSkOp. + ``sketch_symmetric`` now takes a ``blas::Uplo`` and exploits symmetry of :math:`A`: + the ``DenseSkOp`` branch dispatches to ``blas::symm`` (via ``lsksy3`` / ``rsksy3``) + and the ``SparseSkOp`` branch dispatches to a hand-rolled kernel (``lsksys`` / + ``rsksys``) that emits two ``blas::axpy`` calls per stored nonzero of the SkOp, + reading only the triangle of :math:`A` named by ``uplo``. See + ``RandBLAS/sparse_data/DevNotes.md`` for the access-pattern detail. + +Layout-mismatched ``DenseSkOp`` in ``sketch_symmetric`` transpose-copies the operand. + When the ``DenseSkOp``'s storage layout differs from the caller's ``layout`` parameter, ``sketch_symmetric`` transpose-copies the operand into a tight buffer in the caller's layout and then calls ``blas::symm``. The copy is ``O(d * n)``; ``blas::symm`` has no on-the-fly transpose flag for the dense operand, so this is what it costs to keep the SYMM speedup over a ``blas::gemm`` fallback. The layout-matched case skips the copy. Language interoperability diff --git a/rtd/source/api_reference/sketch_dense.rst b/rtd/source/api_reference/sketch_dense.rst index 7dd0a571..1b784b23 100644 --- a/rtd/source/api_reference/sketch_dense.rst +++ b/rtd/source/api_reference/sketch_dense.rst @@ -67,18 +67,26 @@ Analogs to GEMM Analogs to SYMM --------------- +These overloads accept a ``blas::Uplo`` parameter naming the triangle of +:math:`\mtxA` that is structurally stored; the opposite triangle is implied +by symmetry and is **not** read. Both ``DenseSkOp`` and ``SparseSkOp`` are +supported: DenseSkOp dispatches to ``blas::symm`` via ``lsksy3`` / ``rsksy3``, +SparseSkOp dispatches to a hand-rolled two-axpy-per-nonzero kernel via +``lsksys`` / ``rsksys``. See ``RandBLAS/sparse_data/DevNotes.md`` for the +SparseSkOp access-pattern detail. + .. dropdown:: :math:`\mtxB = \alpha \cdot \mtxS \cdot \mtxA + \beta \cdot \mtxB` :animate: fade-in-slide-down :color: light - .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, T alpha, const SKOP &S, const T *A, int64_t lda, T beta, T *B, int64_t ldb, T sym_check_tol = 0) + .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, blas::Uplo uplo, T alpha, const SKOP &S, const T *A, int64_t lda, T beta, T *B, int64_t ldb) :project: RandBLAS .. dropdown:: :math:`\mtxB = \alpha \cdot \mtxA \cdot \mtxS + \beta \cdot \mtxB` :animate: fade-in-slide-down :color: light - .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, T alpha, const T *A, int64_t lda, const SKOP &S, T beta, T *B, int64_t ldb, T sym_check_tol = 0) + .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, blas::Uplo uplo, T alpha, const T *A, int64_t lda, const SKOP &S, T beta, T *B, int64_t ldb) :project: RandBLAS @@ -86,10 +94,10 @@ Analogs to SYMM :animate: fade-in-slide-down :color: light - .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, int64_t d, int64_t n, T alpha, const SKOP &S, int64_t ro_s, int64_t co_s, const T *A, int64_t lda, T beta, T *B, int64_t ldb, T sym_check_tol = 0) + .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, blas::Uplo uplo, int64_t d, int64_t n, T alpha, const SKOP &S, int64_t ro_s, int64_t co_s, const T *A, int64_t lda, T beta, T *B, int64_t ldb) :project: RandBLAS - .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, int64_t n, int64_t d, T alpha, const T *A, int64_t lda, const SKOP &S, int64_t ro_s, int64_t co_s, T beta, T *B, int64_t ldb, T sym_check_tol = 0) + .. doxygenfunction:: RandBLAS::sketch_symmetric(blas::Layout layout, blas::Uplo uplo, int64_t n, int64_t d, T alpha, const T *A, int64_t lda, const SKOP &S, int64_t ro_s, int64_t co_s, T beta, T *B, int64_t ldb) :project: RandBLAS diff --git a/rtd/source/api_reference/sketch_sparse.rst b/rtd/source/api_reference/sketch_sparse.rst index 98a14df9..51890b21 100644 --- a/rtd/source/api_reference/sketch_sparse.rst +++ b/rtd/source/api_reference/sketch_sparse.rst @@ -114,6 +114,47 @@ Deterministic operations This function requires Intel MKL and only supports single and double precision (``float`` and ``double``), in contrast to other RandBLAS kernels that work with any scalar type. +.. dropdown:: :math:`\mtxY = \alpha \cdot \mtxA \cdot \mtxB + \beta \cdot \mtxY,` with sparse symmetric :math:`\mtxA` (spsymm) + :animate: fade-in-slide-down + :color: light + + .. doxygenfunction:: RandBLAS::spsymm(blas::Layout layout, blas::Uplo uplo, int64_t m, int64_t n, T alpha, const SpMat &A, const T *B, int64_t ldb, T beta, T *Y, int64_t ldy) + :project: RandBLAS + + .. doxygenfunction:: RandBLAS::spsymm(blas::Layout layout, int64_t m, int64_t n, T alpha, const Symmetric &A_sym, const T *B, int64_t ldb, T beta, T *Y, int64_t ldy) + :project: RandBLAS + + Only the triangle named by ``uplo`` is read from :math:`\mtxA`. With Intel MKL, + the call dispatches to ``mkl_sparse_d_mm`` with a symmetric ``matrix_descr`` + for the fast path; otherwise it uses a hand-rolled per-format kernel + (CSR / CSC / COO). See ``RandBLAS/sparse_data/DevNotes.md`` for the dispatch + flow and the broader 4-case design that ``spsymm`` belongs to. + + .. note:: + + The ``Symmetric`` carrier (a non-owning wrapper holding a + ``const SpMat &`` and a ``blas::Uplo``) is the type-safety hook that + prevents a symmetric sparse matrix from being accidentally passed to + the general ``spmm`` / ``spgemm`` routines. + + The "dense-symm A times sparse SkOp" case (Case B) is also implemented: + it lives in ``sketch_symmetric`` (the SparseSkOp branch) and uses a + hand-rolled two-axpy-per-nonzero kernel that reads only the named + triangle of A. See ``RandBLAS/sparse_data/DevNotes.md`` for the + access pattern. + + The "sparse-symm A times sparse RHS, dense output" case (Case D) is + implemented via a two-``SparseMatrix``-arg ``spsymm`` overload that + densifies the sparse RHS into a temporary ``m``-by-``n`` buffer + (in the caller's layout) and then calls the Case-C dispatcher on + the densified buffer. This covers all 3 x 3 = 9 sparse-format + pairings for ``(A, B)``. ``mkl_sparse_sp2m`` rejects symmetric + descriptors and ``mkl_sparse_d_spmmd`` takes no descriptor, so + routing through MKL would not avoid the symmetric expansion. The + composition through Case C costs an ``O(m*n)`` temporary, which is + small for the typical RandNLA workload where ``B`` is a sketching + operator with ``nnz(B) << m*n``. + .. 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 cdf202a7..6d966522 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -38,12 +38,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..bce1d076 --- /dev/null +++ b/test/datastructures/test_symmetric_wrapper.cc @@ -0,0 +1,94 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#include "RandBLAS/sparse_data/coo_matrix.hh" +#include "RandBLAS/sparse_data/csr_matrix.hh" +#include "RandBLAS/sparse_data/csc_matrix.hh" +#include "RandBLAS/sparse_data/symmetric.hh" +#include "RandBLAS/exceptions.hh" + +#include +#include + +using namespace RandBLAS::sparse_data; + + +class TestSymmetricWrapper : public ::testing::Test {}; + + +// Construct from each of the three sparse formats and confirm field access + +// trait aliases. The wrapper does not require any actual data, since it only +// stores a reference + uplo. + +TEST_F(TestSymmetricWrapper, ConstructFromCOO) { + COOMatrix A(4, 4); + auto wrapped = as_symmetric(A, blas::Uplo::Upper); + EXPECT_EQ(&wrapped.A, &A); + EXPECT_EQ(wrapped.uplo, blas::Uplo::Upper); + EXPECT_EQ(wrapped.A.n_rows, 4); + EXPECT_EQ(wrapped.A.n_cols, 4); + static_assert(std::is_same_v, + "Symmetric>::scalar_t must be double"); +} + +TEST_F(TestSymmetricWrapper, ConstructFromCSR) { + CSRMatrix A(5, 5); + auto wrapped = as_symmetric(A, blas::Uplo::Lower); + EXPECT_EQ(&wrapped.A, &A); + EXPECT_EQ(wrapped.uplo, blas::Uplo::Lower); + EXPECT_EQ(wrapped.A.n_rows, 5); + static_assert(std::is_same_v, + "Symmetric>::scalar_t must be float"); +} + +TEST_F(TestSymmetricWrapper, ConstructFromCSC) { + CSCMatrix A(3, 3); + auto wrapped = as_symmetric(A, blas::Uplo::Upper); + EXPECT_EQ(&wrapped.A, &A); + EXPECT_EQ(wrapped.uplo, blas::Uplo::Upper); + EXPECT_EQ(wrapped.A.n_rows, 3); +} + +TEST_F(TestSymmetricWrapper, NamespacingExportsAtRandBLASScope) { + // The wrapper and the helper must be accessible via the top-level + // RandBLAS:: namespace (re-exported from sparse_data::). + COOMatrix A(2, 2); + RandBLAS::Symmetric> wrapped_typed(A, blas::Uplo::Upper); + EXPECT_EQ(&wrapped_typed.A, &A); + auto wrapped_sugar = RandBLAS::as_symmetric(A, blas::Uplo::Lower); + EXPECT_EQ(wrapped_sugar.uplo, blas::Uplo::Lower); +} + +TEST_F(TestSymmetricWrapper, NonSquareRejectedAtConstruction) { + COOMatrix A(3, 4); // not square + EXPECT_THROW({ + auto wrapped = as_symmetric(A, blas::Uplo::Upper); + (void) wrapped; + }, RandBLAS::Error); +} diff --git a/test/linops/test_sketch_symmetric.cc b/test/linops/test_sketch_symmetric.cc index 5ee3269f..e7697384 100644 --- a/test/linops/test_sketch_symmetric.cc +++ b/test/linops/test_sketch_symmetric.cc @@ -31,6 +31,7 @@ #include "RandBLAS/base.hh" #include "RandBLAS/random_gen.hh" #include "RandBLAS/dense_skops.hh" +#include "RandBLAS/sparse_skops.hh" #include "RandBLAS/util.hh" #include "RandBLAS/sksy.hh" @@ -39,6 +40,8 @@ using blas::Uplo; using RandBLAS::ScalarDist; using RandBLAS::DenseDist; using RandBLAS::DenseSkOp; +using RandBLAS::SparseDist; +using RandBLAS::SparseSkOp; using RandBLAS::RNGState; using RandBLAS::Axis; @@ -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: an earlier `test_error_on_asymmetric` test verified that + // sketch_symmetric threw on asymmetric input via require_symmetric. With + // the SYMM-based dispatch (project-plans/randblas-symm-plan.md, Phase 1), + // only the triangle named by `uplo` is read --- the opposite triangle + // contributing wrong values is undefined behavior at the caller, not a + // RandBLAS check. The test was removed accordingly. + + // ============================================================================= + // Case B exerciser: sparse SkOp x dense symmetric A. Reference is the + // densified-sparse-skop fed to blas::symm. layout is explicit (SparseSkOp + // has no S.layout the way DenseSkOp does --- the COO storage is layout- + // agnostic, and the test layout determines how both B and the dense + // reference of S are laid out). template - static void test_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; + static void test_sparse_skop( + Layout layout, + uint32_t seed_a, uint32_t seed_skop, Axis major_axis, + T alpha, int64_t d, int64_t n, int64_t lda, T beta, + blas::Side side_skop, + int64_t vec_nnz = 2, + Uplo uplo = Uplo::Upper + ) { + auto [rows_out, cols_out] = dims_of_sketch_symmetric_output(d, n, side_skop); + std::vector A(lda * lda, T(0)); + random_symmetric_mat(n, A.data(), lda, RNGState(seed_a)); + + SparseDist DS(rows_out, cols_out, vec_nnz, major_axis); + SparseSkOp S(DS, seed_skop); + RandBLAS::fill_sparse(S); + + // Densify S into a buffer matching the requested layout, for the SYMM + // reference. lds is the major-axis leading dim. + int64_t lds = (layout == Layout::ColMajor) ? rows_out : cols_out; + int64_t ldb = lds; + std::vector S_dense(static_cast(rows_out) * cols_out, T(0)); + auto Scoo = RandBLAS::coo_view_of_skop(S); + for (int64_t p = 0; p < Scoo.nnz; ++p) { + int64_t r = Scoo.rows[p], c = Scoo.cols[p]; + if (layout == Layout::ColMajor) { + S_dense[r + c * lds] = Scoo.vals[p]; + } else { + S_dense[r * lds + c] = Scoo.vals[p]; + } } - return; + + uint32_t seed_b = seed_a + 42; + std::vector B_actual(static_cast(rows_out) * cols_out); + DenseDist DB(rows_out, cols_out, ScalarDist::Uniform); + RandBLAS::fill_dense_unpacked(layout, DB, rows_out, cols_out, 0, 0, B_actual.data(), RNGState(seed_b)); + std::vector B_expect = B_actual; + + auto side_a = sketch_symmetric_side( + side_skop, layout, uplo, rows_out, cols_out, + alpha, A.data(), lda, S, 0, 0, beta, B_actual.data(), ldb + ); + blas::symm(layout, side_a, uplo, rows_out, cols_out, + alpha, A.data(), lda, S_dense.data(), lds, beta, B_expect.data(), ldb); + + // Same FMA-order-divergence tolerance as the spsymm tests: the sparse + // path emits two AXPYs per stored nonzero (different accumulation + // order than dense SYMM on a fully-stored matrix). + T atol = T(100) * std::numeric_limits::epsilon(); + T rtol = T(10) * std::numeric_limits::epsilon(); + auto msg = RandBLAS::testing::matrices_approx_equal( + layout, blas::Op::NoTrans, rows_out, cols_out, + B_actual.data(), ldb, B_expect.data(), ldb, + __PRETTY_FUNCTION__, __FILE__, __LINE__, atol, rtol + ); + if (!msg.empty()) FAIL() << msg; } }; @@ -391,9 +437,50 @@ TEST_F(TestSketchSymmetric, right_lift_opposing_layouts) { test_opposing_layouts(31, 33, Axis::Long, 0.5, 50, 10, 19, -1.0, blas::Side::Right); } -// MARK: 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 inner-loop axpy strides). + test_sparse_skop(Layout::ColMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, 0.0, blas::Side::Left, 2, Uplo::Lower); + test_sparse_skop(Layout::RowMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, 0.0, blas::Side::Left, 2, Uplo::Lower); + test_sparse_skop(Layout::ColMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, 0.0, blas::Side::Right, 2, Uplo::Lower); + test_sparse_skop(Layout::RowMajor, 0, 1, Axis::Short, 0.5, 3, 10, 10, 0.0, blas::Side::Right, 2, Uplo::Lower); } + +TEST_F(TestSketchSymmetric, sparse_skop_lift) { + // Embedding goes the wrong way (d > n): tests that the kernel handles + // non-square sketches in both directions. + test_sparse_skop(Layout::ColMajor, 0, 1, Axis::Short, 0.5, 13, 10, 10, 0.0, blas::Side::Left, 2, Uplo::Upper); + test_sparse_skop(Layout::ColMajor, 0, 1, Axis::Short, 0.5, 13, 10, 10, 0.0, blas::Side::Right, 2, Uplo::Upper); +} + diff --git a/test/linops/test_spsymm.cc b/test/linops/test_spsymm.cc new file mode 100644 index 00000000..5792df3e --- /dev/null +++ b/test/linops/test_spsymm.cc @@ -0,0 +1,330 @@ +// 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, leaving only the + // structurally-stored side + diagonal. + if (uplo == Uplo::Upper) { + for (int64_t j = 0; j < n; ++j) + for (int64_t i = j + 1; i < n; ++i) + A[i + j * lda] = T(0); + } else { + for (int64_t j = 0; j < n; ++j) + for (int64_t i = 0; i < j; ++i) + A[i + j * lda] = T(0); + } + } + + template + static void dense_to_sparse_format(Layout layout, T* dense_buf, T abs_tol, SpMat& sp) { + using sint_t = typename SpMat::index_t; + if constexpr (std::is_same_v>) { + dense_to_coo(layout, dense_buf, abs_tol, sp); + } else if constexpr (std::is_same_v>) { + dense_to_csr(layout, dense_buf, abs_tol, sp); + } else if constexpr (std::is_same_v>) { + dense_to_csc(layout, dense_buf, abs_tol, sp); + } else { + static_assert(sizeof(SpMat) == 0, "Unsupported sparse format."); + } + } + + template + static void run_case( + Layout layout, Side side, Uplo uplo, + int64_t n_A, int64_t d, + T alpha, T beta, + uint32_t seed_A, uint32_t seed_B, + bool route_via_wrapper = false + ) { + // For side=Left: Y = alpha*A*B + beta*Y, A is n_A x n_A, B and Y are n_A x d. + // For side=Right: Y = alpha*B*A + beta*Y, B and Y are d x n_A, A is n_A x n_A. + int64_t m_BY, n_BY; + if (side == Side::Left) { m_BY = n_A; n_BY = d; } + else { m_BY = d; n_BY = n_A; } + + // Build dense symmetric A. + int64_t lda = n_A; + std::vector A_full(lda * n_A, T(0)); + fill_sym_dense(n_A, A_full.data(), lda, seed_A); + + // Build sparse A: take A_full, zero out the non-named triangle, then convert. + std::vector A_triangle(A_full); + zero_other_triangle(n_A, A_triangle.data(), lda, uplo); + SpMat A_sparse(n_A, n_A); + dense_to_sparse_format(Layout::ColMajor, A_triangle.data(), T(0), A_sparse); + + // Build random B and an initial Y (for beta != 0 to be non-trivial). + int64_t ldb = (layout == Layout::ColMajor) ? m_BY : n_BY; + int64_t ldy = ldb; + std::vector B(m_BY * n_BY); + DenseDist DB(m_BY, n_BY, ScalarDist::Uniform); + RandBLAS::fill_dense_unpacked(layout, DB, m_BY, n_BY, 0, 0, B.data(), RNGState(seed_B)); + + std::vector Y_actual(m_BY * n_BY); + RandBLAS::fill_dense_unpacked(layout, DB, m_BY, n_BY, 0, 0, Y_actual.data(), RNGState(seed_B + 7)); + std::vector Y_expect = Y_actual; + + // Reference using dense blas::symm on the fully-populated A_full + // (both triangles match because A_full is symmetrized; choice of uplo + // for the reference doesn't matter, but we pass `uplo` for consistency). + blas::symm(layout, side, uplo, m_BY, n_BY, + alpha, A_full.data(), lda, B.data(), ldb, + beta, Y_expect.data(), ldy); + + // Under test: spsymm on the one-triangle sparse A. By default we + // call the low-level dispatcher directly; if `route_via_wrapper` is + // set, we go through the public RandBLAS::spsymm(Symmetric) + // overload to exercise the wrapper-routing path. side=Left only + // for the wrapper path since the public wrapper defaults side=Left. + if (route_via_wrapper) { + randblas_require(side == Side::Left); + auto A_sym = RandBLAS::as_symmetric(A_sparse, uplo); + RandBLAS::spsymm(layout, m_BY, n_BY, + alpha, A_sym, B.data(), ldb, + beta, Y_actual.data(), ldy); + } else { + RandBLAS::sparse_data::spsymm(layout, side, uplo, m_BY, n_BY, + alpha, A_sparse, B.data(), ldb, + beta, Y_actual.data(), ldy); + } + + // Tolerance: the dense reference (blas::symm) and the sparse path + // accumulate in different orders, so we get a few ULPs of accumulation + // divergence. Use 100*eps absolute tolerance to absorb that. + T atol = T(100) * std::numeric_limits::epsilon(); + T rtol = T(10) * std::numeric_limits::epsilon(); + auto msg = RandBLAS::testing::matrices_approx_equal( + layout, blas::Op::NoTrans, m_BY, n_BY, + Y_actual.data(), ldy, Y_expect.data(), ldy, + __PRETTY_FUNCTION__, __FILE__, __LINE__, atol, rtol + ); + if (!msg.empty()) FAIL() << msg; + } + + template + static void sweep_layout_uplo(Side side, int64_t n_A, int64_t d, double alpha, double beta) { + using T = typename SpMat::scalar_t; + for (auto layout : {Layout::ColMajor, Layout::RowMajor}) { + for (auto uplo : {Uplo::Upper, Uplo::Lower}) { + SCOPED_TRACE(testing::Message() << + "layout=" << (layout == Layout::ColMajor ? "Col" : "Row") + << " uplo=" << (uplo == Uplo::Upper ? "U" : "L")); + run_case(layout, side, uplo, n_A, d, T(alpha), T(beta), 0, 1); + } + } + } + + // Case D: sparse-symmetric A times sparse B -> dense Y. Reference is + // dense blas::symm on a fully-populated A and a densified B. + template + static void run_case_d( + Layout layout, Side side, Uplo uplo, + int64_t n_A, int64_t d, + T alpha, T beta, + uint32_t seed_A, uint32_t seed_B, + double density_B = 0.3 + ) { + int64_t m_BY, n_BY; + if (side == Side::Left) { m_BY = n_A; n_BY = d; } + else { m_BY = d; n_BY = n_A; } + + // Build dense symm A (full-storage reference). + int64_t lda = n_A; + std::vector A_full(lda * n_A, T(0)); + fill_sym_dense(n_A, A_full.data(), lda, seed_A); + // Build sparse A from a one-triangle copy. + std::vector A_tri(A_full); + zero_other_triangle(n_A, A_tri.data(), lda, uplo); + SpMatA A_sparse(n_A, n_A); + dense_to_sparse_format(Layout::ColMajor, A_tri.data(), T(0), A_sparse); + + // Random sparse B as a ColMajor dense buffer first, then convert to SpMatB. + std::vector B_dense(m_BY * n_BY, T(0)); + { + std::mt19937_64 rng(static_cast(seed_B)); + std::uniform_real_distribution uni01(0.0, 1.0); + std::uniform_real_distribution univ(-1.0, 1.0); + for (int64_t j = 0; j < n_BY; ++j) { + for (int64_t i = 0; i < m_BY; ++i) { + if (uni01(rng) < density_B) { + B_dense[i + j * m_BY] = static_cast(univ(rng)); + } + } + } + } + SpMatB B_sparse(m_BY, n_BY); + dense_to_sparse_format(Layout::ColMajor, B_dense.data(), T(0), B_sparse); + + int64_t ldb = (layout == Layout::ColMajor) ? m_BY : n_BY; + int64_t ldy = ldb; + // The dense reference call to blas::symm uses the requested layout. + std::vector B_dense_layout(m_BY * n_BY); + if (layout == Layout::ColMajor) { + std::copy(B_dense.begin(), B_dense.end(), B_dense_layout.begin()); + } else { + for (int64_t i = 0; i < m_BY; ++i) + for (int64_t j = 0; j < n_BY; ++j) + B_dense_layout[i * ldb + j] = B_dense[i + j * m_BY]; + } + + std::vector Y_actual(m_BY * n_BY); + DenseDist DY(m_BY, n_BY, ScalarDist::Uniform); + RandBLAS::fill_dense_unpacked(layout, DY, m_BY, n_BY, 0, 0, Y_actual.data(), RNGState(seed_B + 13)); + std::vector Y_expect = Y_actual; + + // Reference: dense blas::symm on full-storage A and dense B. + blas::symm(layout, side, uplo, m_BY, n_BY, + alpha, A_full.data(), lda, B_dense_layout.data(), ldb, + beta, Y_expect.data(), ldy); + + // Under test: sparse-symm A times sparse B via Case D. + RandBLAS::sparse_data::spsymm(layout, side, uplo, m_BY, n_BY, + alpha, A_sparse, B_sparse, + beta, Y_actual.data(), ldy); + + T atol = T(100) * std::numeric_limits::epsilon(); + T rtol = T(10) * std::numeric_limits::epsilon(); + auto msg = RandBLAS::testing::matrices_approx_equal( + layout, blas::Op::NoTrans, m_BY, n_BY, + Y_actual.data(), ldy, Y_expect.data(), ldy, + __PRETTY_FUNCTION__, __FILE__, __LINE__, atol, rtol + ); + if (!msg.empty()) FAIL() << msg; + } +}; + + +// 24-cell coverage: {COO, CSR, CSC} x {ColMajor, RowMajor} x {Upper, Lower} x {Left, Right}. +// Each TEST_F sweeps the (layout, uplo) plane internally for one (format, side) combination. + +TEST_F(TestSpsymm, CSR_Left) { sweep_layout_uplo>(Side::Left, 10, 4, 1.5, -0.5); } +TEST_F(TestSpsymm, CSR_Right) { sweep_layout_uplo>(Side::Right, 10, 4, 1.5, -0.5); } +TEST_F(TestSpsymm, CSC_Left) { sweep_layout_uplo>(Side::Left, 10, 4, 1.5, -0.5); } +TEST_F(TestSpsymm, CSC_Right) { sweep_layout_uplo>(Side::Right, 10, 4, 1.5, -0.5); } +TEST_F(TestSpsymm, COO_Left) { sweep_layout_uplo>(Side::Left, 10, 4, 1.5, -0.5); } +TEST_F(TestSpsymm, COO_Right) { sweep_layout_uplo>(Side::Right, 10, 4, 1.5, -0.5); } + +// Float coverage for one representative format +TEST_F(TestSpsymm, CSR_Left_Float) { sweep_layout_uplo>(Side::Left, 10, 4, 1.5, -0.5); } + +// beta=0 (init-from-zero) edge case +TEST_F(TestSpsymm, CSR_Left_Beta0) { + for (auto layout : {Layout::ColMajor, Layout::RowMajor}) + for (auto uplo : {Uplo::Upper, Uplo::Lower}) + run_case>(layout, Side::Left, uplo, 10, 4, 1.0, 0.0, 0, 2); +} + +// alpha=0 (only beta-scaling) edge case +TEST_F(TestSpsymm, CSR_Left_Alpha0) { + run_case>(Layout::ColMajor, Side::Left, Uplo::Upper, 10, 4, 0.0, 0.5, 0, 3); +} + + +// Format-pair sweep: 3 A-formats x 3 B-formats x both sides + an uplo and +// edge-case sample. MKL handles all 9 pairs via make_mkl_handle. +TEST_F(TestSpsymm, CaseD_CSR_CSR_Left) { + run_case_d, CSRMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 1.5, -0.5, 0, 1); + run_case_d, CSRMatrix>(Layout::RowMajor, Side::Left, Uplo::Lower, 8, 3, 1.5, -0.5, 0, 1); +} +TEST_F(TestSpsymm, CaseD_CSC_CSC_Left) { + run_case_d, CSCMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 1.5, -0.5, 0, 1); + run_case_d, CSCMatrix>(Layout::RowMajor, Side::Left, Uplo::Lower, 8, 3, 1.5, -0.5, 0, 1); +} +TEST_F(TestSpsymm, CaseD_COO_COO_Left) { + run_case_d, COOMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 1.5, -0.5, 0, 1); + run_case_d, COOMatrix>(Layout::RowMajor, Side::Left, Uplo::Lower, 8, 3, 1.5, -0.5, 0, 1); +} +TEST_F(TestSpsymm, CaseD_Mixed_Format_Left) { + run_case_d, CSCMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 1.0, 0.0, 0, 1); + run_case_d, CSRMatrix>(Layout::ColMajor, Side::Left, Uplo::Lower, 8, 3, 1.0, 0.0, 0, 1); + run_case_d, CSRMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 1.0, 0.0, 0, 1); +} +TEST_F(TestSpsymm, CaseD_CSR_CSR_Right) { + run_case_d, CSRMatrix>(Layout::ColMajor, Side::Right, Uplo::Upper, 8, 3, 1.5, -0.5, 0, 1); + run_case_d, CSRMatrix>(Layout::RowMajor, Side::Right, Uplo::Lower, 8, 3, 1.5, -0.5, 0, 1); +} +TEST_F(TestSpsymm, CaseD_Float) { + run_case_d, CSRMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 1.5f, -0.5f, 0, 1); +} +TEST_F(TestSpsymm, CaseD_AlphaZero_BetaScale) { + // alpha=0 path: just beta-scales Y, doesn't even touch A or B. + run_case_d, CSRMatrix>(Layout::ColMajor, Side::Left, Uplo::Upper, 8, 3, 0.0, 0.5, 0, 1); +} + +// Routes through the public RandBLAS::spsymm(Symmetric) wrapper +// overload instead of the lower-level RandBLAS::sparse_data::spsymm. +// All other setup (dense reference, comparison tolerance) is identical +// to run_case; we set route_via_wrapper=true to flip the dispatch. +TEST_F(TestSpsymm, SymmetricWrapper) { + 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 + ); +}