diff --git a/.github/actions/setup-randblas-deps/action.yml b/.github/actions/setup-randblas-deps/action.yml index bbbdd294..18182d7d 100644 --- a/.github/actions/setup-randblas-deps/action.yml +++ b/.github/actions/setup-randblas-deps/action.yml @@ -39,6 +39,14 @@ runs: ;; mkl) sudo apt-get install -qq -y libmkl-dev + # Debian packages MKL headers under /usr/include/mkl and sets no + # MKLROOT. blaspp's MKL detection (config/mkl_version.cc) needs the + # headers on the include path or it never defines BLAS_HAVE_MKL, and + # RandBLAS's find_path(mkl_spblas.h) needs MKLROOT to point at the + # header dir. Discover the dir rather than hardcoding it. + mkl_inc="$(dirname "$(find /usr/include -name mkl_spblas.h | head -n1)")" + echo "MKLROOT=${mkl_inc}" >> "$GITHUB_ENV" + echo "CPATH=${mkl_inc}${CPATH:+:$CPATH}" >> "$GITHUB_ENV" ;; auto) sudo apt-get install -qq -y libopenblas-openmp-dev @@ -94,7 +102,7 @@ runs: uses: actions/cache@v4 with: path: ${{ github.workspace }}/../blaspp-install - key: blaspp-${{ runner.os }}-${{ runner.arch }}-${{ inputs.cc }}-${{ inputs.blas-backend }}-${{ steps.blaspp-rev.outputs.sha }} + key: blaspp-mklfix2-${{ runner.os }}-${{ runner.arch }}-${{ inputs.cc }}-${{ inputs.blas-backend }}-${{ steps.blaspp-rev.outputs.sha }} - name: build blaspp if: steps.blaspp-cache.outputs.cache-hit != 'true' @@ -110,9 +118,15 @@ runs: mkdir blaspp-build cd blaspp-build blas_arg="" + # ILP64 only for MKL: it makes MKL_INT 64-bit, which RandBLAS's MKL + # spgemm path (test_spgemm, compiled only under RandBLAS_HAS_MKL) + # requires to match its int64_t sparse indices. Other backends stay + # LP64 -- blaspp's public API is int64_t regardless, so RandBLAS is + # unaffected, and this avoids needing ILP64 OpenBLAS/Accelerate libs. + blas_int_arg="" case "${{ inputs.blas-backend }}" in openblas) blas_arg="-Dblas=openblas" ;; - mkl) blas_arg="-Dblas=mkl" ;; + mkl) blas_arg="-Dblas=mkl"; blas_int_arg="-Dblas_int=int64" ;; accelerate) blas_arg="-Dblas=accelerate" ;; auto) blas_arg="" ;; esac @@ -122,6 +136,7 @@ runs: -DCMAKE_BUILD_TYPE=Release \ -DCMAKE_INSTALL_PREFIX="$(pwd)/../blaspp-install" \ -Dbuild_tests=OFF \ + ${blas_int_arg} \ ${blas_arg} \ ../blaspp if [[ "$(uname)" == "Darwin" ]]; then @@ -137,7 +152,7 @@ runs: uses: actions/cache@v4 with: path: ${{ github.workspace }}/../lapackpp-install - key: lapackpp-${{ runner.os }}-${{ runner.arch }}-${{ inputs.cc }}-${{ inputs.blas-backend }}-${{ steps.blaspp-rev.outputs.sha }}-${{ steps.lapackpp-rev.outputs.sha }} + key: lapackpp-mklfix2-${{ runner.os }}-${{ runner.arch }}-${{ inputs.cc }}-${{ inputs.blas-backend }}-${{ steps.blaspp-rev.outputs.sha }}-${{ steps.lapackpp-rev.outputs.sha }} - name: build lapackpp if: inputs.install-lapackpp == 'true' && steps.lapackpp-cache.outputs.cache-hit != 'true' @@ -152,6 +167,11 @@ runs: git clone --depth 1 https://github.com/icl-utk-edu/lapackpp.git mkdir lapackpp-build cd lapackpp-build + # Match the int width of the blaspp we link against (ILP64 only for MKL). + blas_int_arg="" + if [[ "${{ inputs.blas-backend }}" == "mkl" ]]; then + blas_int_arg="-Dblas_int=int64" + fi cmake \ -DCMAKE_C_COMPILER="${CC}" \ -DCMAKE_CXX_COMPILER="${CXX}" \ @@ -159,6 +179,7 @@ runs: -DCMAKE_INSTALL_PREFIX="$(pwd)/../lapackpp-install" \ -Dblaspp_DIR="$(pwd)/../blaspp-install/lib/cmake/blaspp" \ -Dbuild_tests=OFF \ + ${blas_int_arg} \ ../lapackpp if [[ "$(uname)" == "Darwin" ]]; then jobs="$(sysctl -n hw.logicalcpu)" diff --git a/CMake/MKL_sparse.cmake b/CMake/MKL_sparse.cmake index 7b1abd9e..b4cf22f8 100644 --- a/CMake/MKL_sparse.cmake +++ b/CMake/MKL_sparse.cmake @@ -33,9 +33,15 @@ endif() if (_mkl_found_via_blaspp) # BLAS++ was built with MKL. Find the sparse BLAS header. + # mkl_spblas.h lives in different places depending on the MKL packaging: + # - oneAPI: $MKLROOT/include + # - Debian apt: /usr/include/mkl (and MKLROOT is usually unset) + # Search MKLROOT and the standard system prefixes, trying both the bare + # include dir and an mkl/ subdir. find_path(MKL_SPARSE_INCLUDE_DIR mkl_spblas.h HINTS ENV MKLROOT $ENV{MKLROOT}/include - PATH_SUFFIXES include) + PATHS /usr/include /usr/local/include + PATH_SUFFIXES include mkl include/mkl latest/include) if (MKL_SPARSE_INCLUDE_DIR) set(tmp TRUE) message(STATUS " MKL sparse BLAS header found: ${MKL_SPARSE_INCLUDE_DIR}/mkl_spblas.h") diff --git a/RandBLAS/sparse_data/csr_spmm_impl.hh b/RandBLAS/sparse_data/csr_spmm_impl.hh index cf0338fb..0705bd4c 100644 --- a/RandBLAS/sparse_data/csr_spmm_impl.hh +++ b/RandBLAS/sparse_data/csr_spmm_impl.hh @@ -48,6 +48,40 @@ using RandBLAS::SignedInteger; #define SignedInteger typename #endif +// Core CSR SpMV accumulation: Av[i] += alpha * sum_ell vals[ell] * v[colidxs[ell]]. +// Templated on UnitStride so the hot ColMajor case (incv == incAv == 1) drops the +// runtime index multiplies; `if constexpr` keeps a single source for both variants. +// The simd reduction permits FP reassociation, which breaks the serial dependence +// on the accumulator (no effect when compiled without OpenMP -- the pragma is +// ignored and the loop simply runs unvectorized, as with the rest of this file). +template +static void apply_csr_to_vector_ik_impl( + T alpha, + const T *vals, const sint_t *rowptr, const sint_t *colidxs, + const T *v, int64_t incv, + int64_t len_Av, T *Av, int64_t incAv +) { + UNUSED(incv); UNUSED(incAv); + // ^ silence compiler complaints if UnitStride == true. + for (int64_t i = 0; i < len_Av; ++i) { + T Av_i_diff = 0.0; + #pragma omp simd reduction(+:Av_i_diff) + for (int64_t ell = rowptr[i]; ell < rowptr[i+1]; ++ell) { + int64_t j = colidxs[ell]; + if constexpr (UnitStride) { + Av_i_diff += vals[ell] * v[j]; + } else { + Av_i_diff += vals[ell] * v[j*incv]; + } + } + if constexpr (UnitStride) { + Av[i] += alpha * Av_i_diff; + } else { + Av[i*incAv] += alpha * Av_i_diff; + } + } +} + template static void apply_csr_to_vector_ik( T alpha, @@ -62,14 +96,14 @@ static void apply_csr_to_vector_ik( T *Av, // Av += alpha * A * v. int64_t incAv // stride between elements of Av ) { - for (int64_t i = 0; i < len_Av; ++i) { - T Av_i_diff = 0.0; - for (int64_t ell = rowptr[i]; ell < rowptr[i+1]; ++ell) { - int j = colidxs[ell]; - T Aij = vals[ell]; - Av_i_diff += Aij * v[j*incv]; - } - Av[i*incAv] += alpha * Av_i_diff; + // Unit stride (incv == incAv == 1) is the ColMajor case; specialize on it so + // the inner reduction carries no index multiplies. Other strides fall through. + if (incv == 1 && incAv == 1) { + apply_csr_to_vector_ik_impl( + alpha, vals, rowptr, colidxs, v, incv, len_Av, Av, incAv); + } else { + apply_csr_to_vector_ik_impl( + alpha, vals, rowptr, colidxs, v, incv, len_Av, Av, incAv); } } diff --git a/RandBLAS/sparse_data/mkl_spmm_impl.hh b/RandBLAS/sparse_data/mkl_spmm_impl.hh index 6d2af501..ba65b2f3 100644 --- a/RandBLAS/sparse_data/mkl_spmm_impl.hh +++ b/RandBLAS/sparse_data/mkl_spmm_impl.hh @@ -287,11 +287,6 @@ bool mkl_left_spmm( constexpr bool is_coo = std::is_same_v>; constexpr bool is_csc = std::is_same_v>; - // MKL's mkl_sparse_d_mm does not support CSC format (returns NOT_SUPPORTED). - // Fall back to hand-rolled kernels for CSC matrices. - if constexpr (is_csc) - return false; - // For COO matrices, RandBLAS allows A.n_rows >= d and A.n_cols >= m // (operating on a submatrix). MKL operates on the full matrix, so we // can only use MKL when dimensions match exactly and offsets are zero. @@ -308,21 +303,40 @@ bool mkl_left_spmm( if (opB != blas::Op::NoTrans) return false; - auto h = make_mkl_handle(A); + // Build the MKL sparse handle and the operation to apply to it. + // CSR / COO: the handle wraps A directly; the operation is opA. + // CSC: mkl_sparse_?_mm does not accept CSC, but a CSC matrix's arrays ARE + // a CSR matrix for A^T (transpose_as_csr shares colptr/rowidxs/vals). So + // build the CSR-of-A^T handle and flip the operation (NoTrans <-> Trans), + // which leaves op(A)*B unchanged. MKL supports SPARSE_OPERATION_TRANSPOSE + // on CSR natively, so no explicit transpose or copy is materialized. + MKLSparseHandle h; + sparse_operation_t mkl_op; + if constexpr (is_csc) { + auto At = A.transpose(); // CSRMatrix view of A^T (shares A's arrays) + h = make_mkl_handle_csr(At); + mkl_op = (opA == blas::Op::NoTrans) + ? SPARSE_OPERATION_TRANSPOSE + : SPARSE_OPERATION_NON_TRANSPOSE; + } else { + h = make_mkl_handle(A); + mkl_op = to_mkl_op(opA); + } + 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, + mkl_op, alpha, h.handle, descr, to_mkl_layout(layout), B, (MKL_INT)n, (MKL_INT)ldb, beta, C, (MKL_INT)ldc ); } else if constexpr (std::is_same_v) { status = mkl_sparse_s_mm( - to_mkl_op(opA), alpha, h.handle, descr, + mkl_op, alpha, h.handle, descr, to_mkl_layout(layout), B, (MKL_INT)n, (MKL_INT)ldb, beta, C, (MKL_INT)ldc diff --git a/RandBLAS/sparse_data/spmm_dispatch.hh b/RandBLAS/sparse_data/spmm_dispatch.hh index c2edf112..2c6344fe 100644 --- a/RandBLAS/sparse_data/spmm_dispatch.hh +++ b/RandBLAS/sparse_data/spmm_dispatch.hh @@ -69,21 +69,13 @@ void left_spmm( ) { using blas::Layout; using blas::Op; - // handle applying a transposed sparse matrix. + // Applying a transposed sparse matrix reduces to the NoTrans case on a + // zero-copy transpose view (CSR<->CSC, COO<->COO). MKL, when present, + // engages on the recursive NoTrans call: mkl_left_spmm handles all three + // formats -- including CSC, which it consumes as a CSR-of-transpose view -- + // so the transposed CSR no longer needs to be pre-routed to MKL here to + // avoid a CSC fallback. if (opA == Op::Trans) { - // Try MKL first: it supports SPARSE_OPERATION_TRANSPOSE on CSR natively. - // Without this, the transpose trick below converts CSR→CSC, which MKL - // can't handle, causing unnecessary fallback RandBLAS-defined kernels. - #if defined(RandBLAS_HAS_MKL) - if constexpr (sizeof(typename SpMat::index_t) == sizeof(MKL_INT)) { - bool handled = RandBLAS::sparse_data::mkl::mkl_left_spmm( - layout, opA, opB, d, n, m, alpha, - A, ro_a, co_a, B, ldb, beta, C, ldc - ); - if (handled) - return; - } - #endif auto At = A.transpose(); left_spmm(layout, Op::NoTrans, opB, d, n, m, alpha, At, co_a, ro_a, B, ldb, beta, C, ldc); return; @@ -141,7 +133,8 @@ void left_spmm( #if defined(RandBLAS_HAS_MKL) if constexpr (sizeof(typename SpMat::index_t) == sizeof(MKL_INT)) { // mkl_left_spmm returns false if it can't handle this case - // (e.g., COO with submatrix offsets, CSC format). + // (e.g., COO with submatrix offsets, or opB == Trans). CSC is handled + // via a CSR-of-transpose view inside mkl_left_spmm. // Beta is already applied to C above, so pass beta=1 to MKL // so it adds alpha*A*B to the existing (pre-scaled) C. bool handled = RandBLAS::sparse_data::mkl::mkl_left_spmm( diff --git a/RandBLAS/testing/lapack_like.hh b/RandBLAS/testing/lapack_like.hh index 724effb2..dd07ff98 100644 --- a/RandBLAS/testing/lapack_like.hh +++ b/RandBLAS/testing/lapack_like.hh @@ -61,22 +61,24 @@ void apply_reflectors(int64_t m, int64_t n, int64_t k, const T* vecs, int64_t ld template int potrf_upper_sequential(int64_t n, T* A, int64_t lda) { - // Cache access is much better if the matrix is lower triangular. - // Could implement as lower triangular and then call transpose_square. - for (int64_t j = 0; j < n; ++j) { - if (A[j + j * lda] <= 0) { - std::cout << "Cholesky failed at index " << (j+1) << " of " << n << "."; - return j+1; - } - A[j + j * lda] = std::sqrt(A[j + j * lda]); - for (int64_t i = j + 1; i < n; ++i) { - A[j + i * lda] /= A[j + j * lda]; - } - for (int64_t k = j + 1; k < n; ++k) { - for (int64_t i = k; i < n; ++i) { - A[k + i * lda] -= A[j + i * lda] * A[j + k * lda]; + int64_t j, i, k; + for (j = 0; j < n; ++j) { + T* A_j = A + j*lda; // column j: rows 0..n-1, unit stride + for (i = 0; i < j; ++i) { + const T* A_i = A + i*lda; // column i: already finalized + T s = A_j[i]; + for (int64_t k = 0; k < i; ++k) { + s -= A_i[k] * A_j[k]; // dot of two column-prefixes — both stride 1 } + A_j[i] = s / A_i[i]; // divide by U(i,i), the diagonal of column i } + T d = A_j[j]; + for (k = 0; k < j; ++k) + d -= A_j[k] * A_j[k]; // walks down column j, stride 1 + if (d <= 0) { + return j + 1; // leading minor of order j+1 not SPD + } + A_j[j] = std::sqrt(d); } return 0; } @@ -107,7 +109,9 @@ int potrf_upper(int64_t n, T* A, int64_t lda, int64_t b = 64) { blas::syrk(layout, uplo, blas::Op::Trans, n - curr_b, curr_b, (T) -1.0, A12, lda, (T) 1.0, A22, lda ); - potrf_upper(n-curr_b, A22, lda, b); + if (int info = potrf_upper(n-curr_b, A22, lda, b)) { + return info + curr_b; + } } return 0; } diff --git a/examples/simple-kernel-benchmarks/spmm_performance.cc b/examples/simple-kernel-benchmarks/spmm_performance.cc index ecdcf8e6..84c74e92 100644 --- a/examples/simple-kernel-benchmarks/spmm_performance.cc +++ b/examples/simple-kernel-benchmarks/spmm_performance.cc @@ -4,39 +4,52 @@ // SPARSE-DENSE MATRIX MULTIPLICATION PERFORMANCE BENCHMARK // ============================================================================ // -// This benchmark answers four questions about RandBLAS SpMM performance: +// Benchmarks the RandBLAS left_spmm kernels (C = A*B, with A sparse) across the +// three sparse formats (CSR, CSC, COO) and BOTH dense-matrix layouts +// (ColMajor, RowMajor). Driving left_spmm with both layouts is enough to +// exercise every hand-rolled SpMM kernel, because the dispatch +// (spmm_dispatch.hh) selects the kernel purely from (layout_opB, layout_C): // -// 1. Which sparse format is best for left SpMM (B = S*A)? -// -> CSR, CSC, and COO are compared; CSR and COO (via MKL) tend to win. +// ColMajor dense -> apply_csr_jik_p11 / apply_csc_jki_p11 +// RowMajor dense -> apply_csr_ikb_p1b_rowmajor / apply_csc_kib_1p1_rowmajor +// (COO delegates to the CSC kernels in both layouts, via apply_coo_via_csc) // -// 2. Which sparse format is best for right SpMM (B = A*S)? -// -> CSC wins, because the dispatch transposes CSC to CSR and uses MKL -// with SPARSE_OPERATION_NON_TRANSPOSE — the fastest MKL path. +// A fourth combination -- op(B)=Trans (computing C = A * B^T, the dense operand +// fed transposed) -- also routes to apply_csr_jik_p11 / apply_csc_jki_p11, but +// with a STRIDED dense access (gather when C is ColMajor, strided store when C +// is RowMajor) rather than the unit-stride NoTrans ColMajor pattern. MKL +// declines op(B)=Trans, so it is hand-rolled even on MKL builds. // -// 3. How much does MKL accelerate SpMM over hand-rolled kernels? -// -> For left SpMM, the benchmark runs CSR through both MKL and the -// hand-rolled kernel on the same data, giving a direct speedup ratio. +// These are exactly the kernels that left_spmm and right_spmm TOGETHER reach. +// right_spmm reduces to left_spmm by transposing the sparse operand (CSR<->CSC +// view) and flipping the layout, so the old "right CSR (ColMajor)" path runs +// the very same kernel as "left RowMajor CSC" does here (and old "right CSC" == +// "left RowMajor CSR"). Benchmarking left_spmm directly in both layouts hits +// the same kernels while measuring each one without the right_spmm wrapper. // -// 4. Are left and right SpMM comparable in performance? -// -> The "SQUARE PROBLEMS" section uses d = m = n so both directions -// have identical FLOP counts and output sizes. A summary line after -// each config reports the best-of-each comparison. +// This benchmark answers: +// 1. Which sparse format is fastest, for each dense layout? +// 2. ColMajor vs RowMajor: how much does the dense layout matter? RowMajor +// routes to the BLAS-axpy-based kernels; ColMajor routes to the scalar +// gather (CSR) / scatter (CSC) kernels. +// 3. (MKL builds only) How much does MKL accelerate over the hand-rolled +// kernels? The ColMajor table runs CSR through both left_spmm (which uses +// MKL when available) and the hand-rolled kernel called directly. This is +// the one path that does NOT go through left_spmm -- on an MKL build the +// dispatch always picks MKL, so the hand-rolled kernel is otherwise +// unreachable. It compiles out entirely when MKL is disabled. +// 4. op(B)=Trans (C = A * B^T): how much does the strided dense access cost +// relative to the unit-stride NoTrans path? Reported in two extra tables. // -// NOTE: This benchmarks the SpMM kernel directly (left_spmm / right_spmm), -// not sketch_general. Any overhead from SKOP generation or sketch_general's -// internal dispatch is not captured here. In practice SpMM dominates the -// cost, so these results are a good proxy for sketching performance. +// NOTE: benchmarks the SpMM kernel directly, not sketch_general. SpMM dominates +// sketching cost, so these results are a good proxy for sketching performance. // // NOTATION: -// S - Sparse matrix (m x n) -// A - Dense input matrix -// B - Dense result matrix -// -// Left SpMM: B(m x d) = S(m x n) * A(n x d) -// Right SpMM: B(d x n) = A(d x m) * S(m x n) -// -// A "densify + GEMM" reference converts S to dense and calls BLAS GEMM, -// reporting the densify and GEMM times separately. +// A - sparse (m x n), B - dense (n x d), C - dense result (m x d) +// left_spmm: C(m x d) = A(m x n) * B(n x d) +// The same logical B (and hence result C) is stored in BOTH ColMajor and +// RowMajor, so the two layouts compute the identical product and are +// directly comparable. // // USAGE: // ./spmm_performance # default sweep @@ -44,8 +57,8 @@ // // Default sweep (no arguments): // Two sections, 10 trials each: -// 1. SQUARE (d = m = n): sizes 100..2000, fair left-vs-right comparison -// 2. RECTANGULAR: square S with small d, tall S, wide S +// 1. SQUARE (d = m = n): sizes 100..2000 +// 2. RECTANGULAR: square A with small d, tall A, wide A // // Single config (4+ arguments): // One (m, n, d, density) configuration, default 20 trials. @@ -79,8 +92,10 @@ using blas::Layout; using blas::Op; // Calls the hand-rolled CSR left-multiply kernel directly, bypassing MKL. -// Used to answer question 3 (MKL vs hand-rolled speedup). Replicates the -// beta-scaling and kernel selection logic from spmm_dispatch.hh. +// Used to answer question 3 (MKL vs hand-rolled speedup) on MKL builds, where +// left_spmm would otherwise always pick MKL. Replicates the beta-scaling and +// kernel selection logic from spmm_dispatch.hh. Here A is sparse (CSR), B is the +// dense input, C is the dense result -- matching the left_spmm convention. template void handrolled_left_spmm_csr( blas::Layout layout, int64_t d, int64_t n, int64_t m, @@ -169,156 +184,170 @@ void print_densify_row(const std::string& name, long total, long dens, long gemm << " (densify " << dens << " + GEMM " << gemm << ")\n"; } -// Run one complete benchmark configuration: generate matrices, time all -// format x direction combinations, verify correctness, and print results. +void print_table_header(const std::string& title) { + std::cout << " " << title << "\n"; + std::cout << " " << std::setw(24) << std::left << "Kernel" + << std::setw(10) << std::right << "Min (us)" + << std::setw(10) << "Med (us)" + << std::setw(11) << "vs best" << "\n"; + std::cout << " " << std::string(55, '-') << "\n"; +} + +// Run one complete benchmark configuration: generate matrices, time every +// format x layout combination through left_spmm, verify correctness, print. void run_config(int64_t m, int64_t n, int64_t d, double density, int num_trials) { using T = double; + namespace rbsd = RandBLAS::sparse_data; uint64_t seed = 12345; // Header std::string shape; if (m == n && d == m) shape = "all square"; - else if (m == n) shape = "square S"; + else if (m == n) shape = "square A"; else if (m > n) shape = "tall"; else shape = "wide"; - std::cout << "--- S(" << m << "x" << n << "), d=" << d + std::cout << "--- A(" << m << "x" << n << "), d=" << d << ", density=" << std::setprecision(4) << density << " (" << shape << ") ---\n"; // Generate one random COO matrix and convert to CSR/CSC so all three // formats represent the identical mathematical matrix. - auto [S_coo, next_state] = RandBLAS::testing::random_coo(m, n, density, RandBLAS::RNGState<>(seed)); - auto S_csr = S_coo.as_owning_csr(); - auto S_csc = S_coo.as_owning_csc(); - - std::cout << " nnz=" << S_coo.nnz << ", trials=" << num_trials << "\n\n"; - - // Generate dense matrices - auto state = next_state; - std::vector A_left(n * d); - RandBLAS::DenseDist D_left(n, d); - state = RandBLAS::fill_dense(D_left, A_left.data(), state); - - std::vector A_right(d * m); - RandBLAS::DenseDist D_right(d, m); - state = RandBLAS::fill_dense(D_right, A_right.data(), state); - - // Result and workspace buffers - std::vector B_left(m * d); - std::vector B_right(d * n); - std::vector S_dense(m * n); - - // ---- Left SpMM: B = S * A ---- - auto [min_l_csr, med_l_csr] = run_trials([&]() { - std::fill(B_left.begin(), B_left.end(), 0.0); - RandBLAS::sparse_data::left_spmm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, - m, d, n, 1.0, S_csr, 0, 0, A_left.data(), n, 0.0, B_left.data(), m); - }, num_trials); - - // Hand-rolled CSR: same data as above, but bypasses MKL (question 3). - long min_l_csr_hr = 0, med_l_csr_hr = 0; + auto [A_coo, next_state] = RandBLAS::testing::random_coo(m, n, density, RandBLAS::RNGState<>(seed)); + auto A_csr = A_coo.as_owning_csr(); + auto A_csc = A_coo.as_owning_csc(); + + std::cout << " nnz=" << A_coo.nnz << ", trials=" << num_trials << "\n\n"; + + // One logical dense B (n x d), stored in both layouts so the ColMajor and + // RowMajor runs compute the identical product C = A*B. + // ColMajor: B(i,j) = B_cm[i + j*n] (ldB = n) + // RowMajor: B(i,j) = B_rm[i*d + j] (ldB = d) + std::vector B_cm(n * d); + RandBLAS::DenseDist DB(n, d); + RandBLAS::fill_dense(DB, B_cm.data(), next_state); + std::vector B_rm(n * d); + for (int64_t i = 0; i < n; ++i) + for (int64_t j = 0; j < d; ++j) + B_rm[i*d + j] = B_cm[i + j*n]; + + // Result buffers (logical C is m x d in both layouts) and densify workspace. + std::vector C_cm(m * d), C_rm(m * d); + std::vector A_dense(m * n); + + // Time left_spmm for one (format, layout): C = A * B, with B and C stored + // in `layout`. ldB / ldC follow from the layout (see notation above). + auto time_spmm = [&](const auto& A, Layout layout, std::vector& C) { + const T* Bptr = (layout == Layout::ColMajor) ? B_cm.data() : B_rm.data(); + int64_t ldB = (layout == Layout::ColMajor) ? n : d; + int64_t ldC = (layout == Layout::ColMajor) ? m : d; + return run_trials([&]() { + std::fill(C.begin(), C.end(), 0.0); + rbsd::left_spmm(layout, Op::NoTrans, Op::NoTrans, m, d, n, + 1.0, A, 0, 0, Bptr, ldB, 0.0, C.data(), ldC); + }, num_trials); + }; + + // ---- ColMajor left_spmm (hits jik / jki kernels) ---- + auto [min_cm_csr, med_cm_csr] = time_spmm(A_csr, Layout::ColMajor, C_cm); + auto [min_cm_csc, med_cm_csc] = time_spmm(A_csc, Layout::ColMajor, C_cm); + auto [min_cm_coo, med_cm_coo] = time_spmm(A_coo, Layout::ColMajor, C_cm); + + // Hand-rolled CSR ColMajor: bypasses MKL for the MKL-vs-hand-rolled + // comparison. This is the only path that does not go through left_spmm. #if defined(RandBLAS_HAS_MKL) + long min_cm_csr_hr = 0, med_cm_csr_hr = 0; { auto [mn, md] = run_trials([&]() { - std::fill(B_left.begin(), B_left.end(), 0.0); + std::fill(C_cm.begin(), C_cm.end(), 0.0); handrolled_left_spmm_csr(Layout::ColMajor, m, d, n, 1.0, - S_csr, A_left.data(), n, 0.0, B_left.data(), m); + A_csr, B_cm.data(), n, 0.0, C_cm.data(), m); }, num_trials); - min_l_csr_hr = mn; - med_l_csr_hr = md; + min_cm_csr_hr = mn; + med_cm_csr_hr = md; } #endif - auto [min_l_csc, med_l_csc] = run_trials([&]() { - std::fill(B_left.begin(), B_left.end(), 0.0); - RandBLAS::sparse_data::left_spmm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, - m, d, n, 1.0, S_csc, 0, 0, A_left.data(), n, 0.0, B_left.data(), m); - }, num_trials); - - auto [min_l_coo, med_l_coo] = run_trials([&]() { - std::fill(B_left.begin(), B_left.end(), 0.0); - RandBLAS::sparse_data::left_spmm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, - m, d, n, 1.0, S_coo, 0, 0, A_left.data(), n, 0.0, B_left.data(), m); - }, num_trials); - - // ---- Right SpMM: B = A * S ---- - auto [min_r_csr, med_r_csr] = run_trials([&]() { - std::fill(B_right.begin(), B_right.end(), 0.0); - RandBLAS::sparse_data::right_spmm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, - d, n, m, 1.0, A_right.data(), d, S_csr, 0, 0, 0.0, B_right.data(), d); - }, num_trials); - - auto [min_r_csc, med_r_csc] = run_trials([&]() { - std::fill(B_right.begin(), B_right.end(), 0.0); - RandBLAS::sparse_data::right_spmm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, - d, n, m, 1.0, A_right.data(), d, S_csc, 0, 0, 0.0, B_right.data(), d); - }, num_trials); - - auto [min_r_coo, med_r_coo] = run_trials([&]() { - std::fill(B_right.begin(), B_right.end(), 0.0); - RandBLAS::sparse_data::right_spmm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, - d, n, m, 1.0, A_right.data(), d, S_coo, 0, 0, 0.0, B_right.data(), d); - }, num_trials); - - // ---- Reference: densify + GEMM ---- - auto [min_dens, min_gemm_left] = run_split_trials( - [&]() { RandBLAS::sparse_data::csr::csr_to_dense(S_csr, Layout::ColMajor, S_dense.data()); }, + // ---- RowMajor left_spmm (hits ikb / kib kernels; the paths the old + // right_spmm section reached, via the CSR<->CSC transpose) ---- + auto [min_rm_csr, med_rm_csr] = time_spmm(A_csr, Layout::RowMajor, C_rm); + auto [min_rm_csc, med_rm_csc] = time_spmm(A_csc, Layout::RowMajor, C_rm); + auto [min_rm_coo, med_rm_coo] = time_spmm(A_coo, Layout::RowMajor, C_rm); + + // ---- Transposed dense (opB=Trans): C = A * B^T, dense operand fed + // transposed so the SAME product A*B is computed. This routes to the + // jik/jki kernels with a STRIDED dense access (incv != 1 when C is + // ColMajor, incAv != 1 when C is RowMajor) -- distinct from the + // unit-stride NoTrans ColMajor path, and the only sparse-NoTrans path + // MKL declines (so it stays hand-rolled on MKL builds). + // B^T (d x n) needs no new storage: the transpose of the ColMajor n-x-d + // buffer IS the RowMajor buffer, so ColMajor-B^T == B_rm, RowMajor == B_cm. + auto time_spmm_t = [&](const auto& A, Layout layout, std::vector& C) { + const T* BTptr = (layout == Layout::ColMajor) ? B_rm.data() : B_cm.data(); + int64_t ldBT = (layout == Layout::ColMajor) ? d : n; + int64_t ldC = (layout == Layout::ColMajor) ? m : d; + return run_trials([&]() { + std::fill(C.begin(), C.end(), 0.0); + rbsd::left_spmm(layout, Op::NoTrans, Op::Trans, m, d, n, + 1.0, A, 0, 0, BTptr, ldBT, 0.0, C.data(), ldC); + }, num_trials); + }; + + auto [min_tcm_csr, med_tcm_csr] = time_spmm_t(A_csr, Layout::ColMajor, C_cm); + auto [min_tcm_csc, med_tcm_csc] = time_spmm_t(A_csc, Layout::ColMajor, C_cm); + auto [min_tcm_coo, med_tcm_coo] = time_spmm_t(A_coo, Layout::ColMajor, C_cm); + auto [min_trm_csr, med_trm_csr] = time_spmm_t(A_csr, Layout::RowMajor, C_rm); + auto [min_trm_csc, med_trm_csc] = time_spmm_t(A_csc, Layout::RowMajor, C_rm); + auto [min_trm_coo, med_trm_coo] = time_spmm_t(A_coo, Layout::RowMajor, C_rm); + + // ---- Reference: densify + GEMM, per layout ---- + // The ColMajor reference (ref_cm) doubles as the correctness oracle. + std::vector ref_cm(m * d); + auto [min_dens_cm, min_gemm_cm] = run_split_trials( + [&]() { rbsd::csr::csr_to_dense(A_csr, Layout::ColMajor, A_dense.data()); }, [&]() { - std::fill(B_left.begin(), B_left.end(), 0.0); + std::fill(ref_cm.begin(), ref_cm.end(), 0.0); blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, - m, d, n, 1.0, S_dense.data(), m, A_left.data(), n, 0.0, B_left.data(), m); + m, d, n, 1.0, A_dense.data(), m, B_cm.data(), n, 0.0, ref_cm.data(), m); }, num_trials); - long min_ref_left = min_dens + min_gemm_left; + long min_ref_cm = min_dens_cm + min_gemm_cm; - auto [min_dens2, min_gemm_right] = run_split_trials( - [&]() { RandBLAS::sparse_data::csr::csr_to_dense(S_csr, Layout::ColMajor, S_dense.data()); }, + auto [min_dens_rm, min_gemm_rm] = run_split_trials( + [&]() { rbsd::csr::csr_to_dense(A_csr, Layout::RowMajor, A_dense.data()); }, [&]() { - std::fill(B_right.begin(), B_right.end(), 0.0); - blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, - d, n, m, 1.0, A_right.data(), d, S_dense.data(), m, 0.0, B_right.data(), d); + std::fill(C_rm.begin(), C_rm.end(), 0.0); + blas::gemm(Layout::RowMajor, Op::NoTrans, Op::NoTrans, + m, d, n, 1.0, A_dense.data(), n, B_rm.data(), d, 0.0, C_rm.data(), d); }, num_trials); - long min_ref_right = min_dens + min_gemm_right; + long min_ref_rm = min_dens_rm + min_gemm_rm; - // ---- Correctness: compare each SpMM path against densify + GEMM ---- + // ---- Correctness: re-run each (format, layout) once, compare to ref_cm ---- bool all_pass = true; int num_checks = 0; - - auto check = [&](const std::string& label, auto& sparse, auto& dense_in, auto& result, - int64_t r, int64_t c, int64_t k, bool is_left) { - // Densify this sparse matrix and compute reference via GEMM - std::vector S_dens(m * n); - if constexpr (std::is_same_v, - RandBLAS::sparse_data::CSRMatrix>) { - RandBLAS::sparse_data::csr::csr_to_dense(sparse, Layout::ColMajor, S_dens.data()); - } else if constexpr (std::is_same_v, - RandBLAS::sparse_data::CSCMatrix>) { - RandBLAS::sparse_data::csc::csc_to_dense(sparse, Layout::ColMajor, S_dens.data()); - } else { - RandBLAS::sparse_data::coo::coo_to_dense(sparse, Layout::ColMajor, S_dens.data()); - } - - std::vector ref(r * c, 0.0); - if (is_left) { - blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, - r, c, k, 1.0, S_dens.data(), m, dense_in.data(), n, 0.0, ref.data(), m); + std::vector result(m * d); + + // All NoTrans and Trans paths compute the same logical product A*B, so they + // are all checked against the single ColMajor oracle ref_cm. For opB=Trans + // the dense operand is B^T (== the other layout's buffer; see timing above). + auto check = [&](const std::string& label, const auto& A, Layout layout, Op opB) { + const T* Bptr; + int64_t ldB; + if (opB == Op::NoTrans) { + Bptr = (layout == Layout::ColMajor) ? B_cm.data() : B_rm.data(); + ldB = (layout == Layout::ColMajor) ? n : d; } else { - blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, - r, c, k, 1.0, dense_in.data(), d, S_dens.data(), m, 0.0, ref.data(), d); + Bptr = (layout == Layout::ColMajor) ? B_rm.data() : B_cm.data(); + ldB = (layout == Layout::ColMajor) ? d : n; } - - // Compute via SpMM and compare + int64_t ldC = (layout == Layout::ColMajor) ? m : d; std::fill(result.begin(), result.end(), 0.0); - if (is_left) { - RandBLAS::sparse_data::left_spmm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, - r, c, k, 1.0, sparse, 0, 0, dense_in.data(), n, 0.0, result.data(), m); - } else { - RandBLAS::sparse_data::right_spmm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, - r, c, k, 1.0, dense_in.data(), d, sparse, 0, 0, 0.0, result.data(), d); - } - + rbsd::left_spmm(layout, Op::NoTrans, opB, m, d, n, + 1.0, A, 0, 0, Bptr, ldB, 0.0, result.data(), ldC); double maxdiff = 0; - for (int64_t i = 0; i < r * c; ++i) - maxdiff = std::max(maxdiff, std::abs(result[i] - ref[i])); + for (int64_t i = 0; i < m; ++i) + for (int64_t j = 0; j < d; ++j) { + int64_t idx = (layout == Layout::ColMajor) ? (i + j*m) : (i*d + j); + maxdiff = std::max(maxdiff, std::abs(result[idx] - ref_cm[i + j*m])); + } num_checks++; if (maxdiff > 1e-10) { std::cout << " FAIL: " << label << " max|diff|=" << std::scientific << maxdiff << "\n"; @@ -326,32 +355,31 @@ void run_config(int64_t m, int64_t n, int64_t d, double density, int num_trials) } }; - check("left+CSR", S_csr, A_left, B_left, m, d, n, true); - check("left+CSC", S_csc, A_left, B_left, m, d, n, true); - check("left+COO", S_coo, A_left, B_left, m, d, n, true); - check("right+CSR", S_csr, A_right, B_right, d, n, m, false); - check("right+CSC", S_csc, A_right, B_right, d, n, m, false); - check("right+COO", S_coo, A_right, B_right, d, n, m, false); - - // Also verify hand-rolled CSR correctness + check("ColMajor+CSR", A_csr, Layout::ColMajor, Op::NoTrans); + check("ColMajor+CSC", A_csc, Layout::ColMajor, Op::NoTrans); + check("ColMajor+COO", A_coo, Layout::ColMajor, Op::NoTrans); + check("RowMajor+CSR", A_csr, Layout::RowMajor, Op::NoTrans); + check("RowMajor+CSC", A_csc, Layout::RowMajor, Op::NoTrans); + check("RowMajor+COO", A_coo, Layout::RowMajor, Op::NoTrans); + check("ColMajor+CSR(B^T)", A_csr, Layout::ColMajor, Op::Trans); + check("ColMajor+CSC(B^T)", A_csc, Layout::ColMajor, Op::Trans); + check("ColMajor+COO(B^T)", A_coo, Layout::ColMajor, Op::Trans); + check("RowMajor+CSR(B^T)", A_csr, Layout::RowMajor, Op::Trans); + check("RowMajor+CSC(B^T)", A_csc, Layout::RowMajor, Op::Trans); + check("RowMajor+COO(B^T)", A_coo, Layout::RowMajor, Op::Trans); + + // Also verify the hand-rolled CSR path. #if defined(RandBLAS_HAS_MKL) { - std::vector S_dens(m * n); - RandBLAS::sparse_data::csr::csr_to_dense(S_csr, Layout::ColMajor, S_dens.data()); - std::vector ref(m * d, 0.0); - blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, - m, d, n, 1.0, S_dens.data(), m, A_left.data(), n, 0.0, ref.data(), m); - - std::fill(B_left.begin(), B_left.end(), 0.0); + std::fill(result.begin(), result.end(), 0.0); handrolled_left_spmm_csr(Layout::ColMajor, m, d, n, 1.0, - S_csr, A_left.data(), n, 0.0, B_left.data(), m); - + A_csr, B_cm.data(), n, 0.0, result.data(), m); double maxdiff = 0; for (int64_t i = 0; i < m * d; ++i) - maxdiff = std::max(maxdiff, std::abs(B_left[i] - ref[i])); + maxdiff = std::max(maxdiff, std::abs(result[i] - ref_cm[i])); num_checks++; if (maxdiff > 1e-10) { - std::cout << " FAIL: left+CSR(hand-rolled) max|diff|=" << std::scientific << maxdiff << "\n"; + std::cout << " FAIL: ColMajor+CSR(hand-rolled) max|diff|=" << std::scientific << maxdiff << "\n"; all_pass = false; } } @@ -364,62 +392,61 @@ void run_config(int64_t m, int64_t n, int64_t d, double density, int num_trials) // ---- Results tables ---- std::cout << std::fixed << std::setprecision(2); - // Left SpMM table + // ColMajor table #if defined(RandBLAS_HAS_MKL) - long bl = std::min({min_l_csr, min_l_csr_hr, min_l_csc, min_l_coo}); + long bcm = std::min({min_cm_csr, min_cm_csr_hr, min_cm_csc, min_cm_coo}); #else - long bl = std::min({min_l_csr, min_l_csc, min_l_coo}); + long bcm = std::min({min_cm_csr, min_cm_csc, min_cm_coo}); #endif - - std::cout << " LEFT SPMM: B(" << m << "x" << d << ") = S * A\n"; - std::cout << " " << std::setw(24) << std::left << "Kernel" - << std::setw(10) << std::right << "Min (us)" - << std::setw(10) << "Med (us)" - << std::setw(11) << "vs best" << "\n"; - std::cout << " " << std::string(55, '-') << "\n"; - + print_table_header("LEFT SPMM, ColMajor dense: C(" + std::to_string(m) + "x" + + std::to_string(d) + ") = A * B (jik/jki kernels)"); #if defined(RandBLAS_HAS_MKL) - print_row("CSR (MKL)", min_l_csr, med_l_csr, bl); - print_row("CSR (hand-rolled)", min_l_csr_hr, med_l_csr_hr, bl); - print_row("CSC (hand-rolled)", min_l_csc, med_l_csc, bl); - print_row("COO (MKL)", min_l_coo, med_l_coo, bl); + print_row("CSR (MKL)", min_cm_csr, med_cm_csr, bcm); + print_row("CSR (hand-rolled)", min_cm_csr_hr, med_cm_csr_hr, bcm); #else - print_row("CSR (hand-rolled)", min_l_csr, med_l_csr, bl); - print_row("CSC (hand-rolled)", min_l_csc, med_l_csc, bl); - print_row("COO (hand-rolled)", min_l_coo, med_l_coo, bl); + print_row("CSR", min_cm_csr, med_cm_csr, bcm); #endif - print_densify_row("densify+GEMM", min_ref_left, min_dens, min_gemm_left, bl); + print_row("CSC", min_cm_csc, med_cm_csc, bcm); + print_row("COO", min_cm_coo, med_cm_coo, bcm); + print_densify_row("densify+GEMM", min_ref_cm, min_dens_cm, min_gemm_cm, bcm); std::cout << "\n"; - // Right SpMM table - long br = std::min({min_r_csr, min_r_csc, min_r_coo}); + // RowMajor table + long brm = std::min({min_rm_csr, min_rm_csc, min_rm_coo}); + print_table_header("LEFT SPMM, RowMajor dense: C(" + std::to_string(m) + "x" + + std::to_string(d) + ") = A * B (ikb/kib kernels)"); + print_row("CSR", min_rm_csr, med_rm_csr, brm); + print_row("CSC", min_rm_csc, med_rm_csc, brm); + print_row("COO", min_rm_coo, med_rm_coo, brm); + print_densify_row("densify+GEMM", min_ref_rm, min_dens_rm, min_gemm_rm, brm); + std::cout << "\n"; - std::cout << " RIGHT SPMM: B(" << d << "x" << n << ") = A * S\n"; - std::cout << " " << std::setw(24) << std::left << "Kernel" - << std::setw(10) << std::right << "Min (us)" - << std::setw(10) << "Med (us)" - << std::setw(11) << "vs best" << "\n"; - std::cout << " " << std::string(55, '-') << "\n"; + // Transposed-dense tables (opB=Trans -> jik/jki with strided dense access). + // MKL declines opB=Trans, so these stay hand-rolled even on MKL builds. + long btc = std::min({min_tcm_csr, min_tcm_csc, min_tcm_coo}); + print_table_header("LEFT SPMM, transposed dense, ColMajor C: C = A * B^T (opB=Trans -> jik/jki, strided gather)"); + print_row("CSR", min_tcm_csr, med_tcm_csr, btc); + print_row("CSC", min_tcm_csc, med_tcm_csc, btc); + print_row("COO", min_tcm_coo, med_tcm_coo, btc); + std::cout << "\n"; - #if defined(RandBLAS_HAS_MKL) - print_row("CSR (MKL)", min_r_csr, med_r_csr, br); - print_row("CSC (MKL via CSR)", min_r_csc, med_r_csc, br); - print_row("COO (MKL)", min_r_coo, med_r_coo, br); - #else - print_row("CSR (hand-rolled)", min_r_csr, med_r_csr, br); - print_row("CSC (hand-rolled)", min_r_csc, med_r_csc, br); - print_row("COO (hand-rolled)", min_r_coo, med_r_coo, br); - #endif - print_densify_row("densify+GEMM", min_ref_right, min_dens, min_gemm_right, br); + long btr = std::min({min_trm_csr, min_trm_csc, min_trm_coo}); + print_table_header("LEFT SPMM, transposed dense, RowMajor C: C = A * B^T (opB=Trans -> jik/jki, strided store)"); + print_row("CSR", min_trm_csr, med_trm_csr, btr); + print_row("CSC", min_trm_csc, med_trm_csc, btr); + print_row("COO", min_trm_coo, med_trm_coo, btr); std::cout << "\n"; - // Summary: best left vs best right - std::cout << " SUMMARY: best left " << bl << " us | best right " << br << " us | "; - if (bl <= br) - std::cout << "left " << std::fixed << std::setprecision(2) << (double)br / bl << "x faster\n"; + // Summary: best NoTrans (ColMajor vs RowMajor) and the transposed-dense cost. + long bnt = std::min(bcm, brm); + long bt = std::min(btc, btr); + std::cout << " SUMMARY: best NoTrans ColMajor " << bcm << " us | best NoTrans RowMajor " << brm << " us | "; + if (brm <= bcm) + std::cout << "RowMajor " << std::fixed << std::setprecision(2) << (double)bcm / brm << "x faster\n"; else - std::cout << "right " << std::fixed << std::setprecision(2) << (double)bl / br << "x faster\n"; - std::cout << "\n"; + std::cout << "ColMajor " << std::fixed << std::setprecision(2) << (double)brm / bcm << "x faster\n"; + std::cout << " best transposed-dense (opB=Trans) " << bt << " us (" + << std::fixed << std::setprecision(2) << (double)bt / bnt << "x vs best NoTrans)\n\n"; } int main(int argc, char** argv) { @@ -433,9 +460,10 @@ int main(int argc, char** argv) { std::cout << "MKL support: DISABLED\n"; #endif std::cout << "\n"; - std::cout << " S is m-by-n (sparse), A and B are dense.\n"; - std::cout << " Left SpMM: B(m x d) = S(m x n) * A(n x d)\n"; - std::cout << " Right SpMM: B(d x n) = A(d x m) * S(m x n)\n\n"; + std::cout << " A is m-by-n (sparse), B and C are dense.\n"; + std::cout << " left_spmm: C(m x d) = A(m x n) * B(n x d)\n"; + std::cout << " All kernels are reached through left_spmm; the ColMajor and\n"; + std::cout << " RowMajor dense layouts select the two kernel families.\n\n"; if (argc >= 5) { // Single config mode diff --git a/test/test_exceptions.cc b/test/test_exceptions.cc index 28aa8f13..3fa04b2a 100644 --- a/test/test_exceptions.cc +++ b/test/test_exceptions.cc @@ -28,7 +28,7 @@ TEST_F(TestExceptions, randblas_require_expr_arg) { randblas_require(flag > 1); } catch (RandBLAS::Error &e) { std::string message{e.what()}; - flag = (int) message.find("flag > 1") != std::string::npos; + flag = message.find("flag > 1") != std::string::npos; } ASSERT_TRUE(flag); }