Kernel speed improvements

.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,13 +167,19 @@ 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}" \
           -DCMAKE_BUILD_TYPE=Release \
           -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)"

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")

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 <bool UnitStride, typename T, SignedInteger sint_t>
+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 <typename T, SignedInteger sint_t = int64_t>
 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<true>(
+            alpha, vals, rowptr, colidxs, v, incv, len_Av, Av, incAv);
+    } else {
+        apply_csr_to_vector_ik_impl<false>(
+            alpha, vals, rowptr, colidxs, v, incv, len_Av, Av, incAv);
     }
 }
 

RandBLAS/sparse_data/mkl_spmm_impl.hh

@@ -287,11 +287,6 @@ bool mkl_left_spmm(
     constexpr bool is_coo = std::is_same_v<SpMat, COOMatrix<T, sint_t>>;
     constexpr bool is_csc = std::is_same_v<SpMat, CSCMatrix<T, sint_t>>;
 
-    // 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<T, double>) {
         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<T, float>) {
         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

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(

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 <typename T>
 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;
 }