From 29b7e067b12528473ed1e75077d8d55b4a11e04a Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 14 Jul 2026 09:27:46 -0400 Subject: [PATCH 1/3] rewrite blas_mul! for CPU so less methods have to be compiled --- src/linalg.jl | 65 ++++++++++++++++++++++++--------------------------- 1 file changed, 31 insertions(+), 34 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index cb56c51..c5befd1 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -111,50 +111,47 @@ function blas_mul!( C::StridedView{T, 2}, A::StridedView{T, 2}, B::StridedView{T, 2}, α::Number, β::Number ) where {T <: LinearAlgebra.BlasFloat} + m, n = size(C) + m == size(A, 1) && n == size(B, 2) || throw(DimensionMismatch()) nthreads = use_threaded_mul() ? get_num_threads() : 1 - return _threaded_blas_mul!(C, A, B, α, β, nthreads) + A2, CA = getblasmatrix(A) + B2, CB = getblasmatrix(B) + return _threaded_blas_mul!(C, A2, CA, B2, CB, convert(T, α), convert(T, β), nthreads) end +_split_row(A, m) = (A[1:m, :], A[(m + 1):end, :]) +_split_col(A, m) = (A[:, 1:m], A[:, (m + 1):end]) + function _threaded_blas_mul!( - C::StridedView{T, 2}, A::StridedView{T, 2}, B::StridedView{T, 2}, + C::StridedView{T, 2}, + A::StridedView{T, 2}, CA, + B::StridedView{T, 2}, CB, α::Number, β::Number, nthreads ) where {T <: LinearAlgebra.BlasFloat} m, n = size(C) - m == size(A, 1) && n == size(B, 2) || throw(DimensionMismatch()) - return if nthreads == 1 || m * n < 1024 - A2, CA = getblasmatrix(A) - B2, CB = getblasmatrix(B) - LinearAlgebra.BLAS.gemm!(CA, CB, convert(T, α), A2, B2, convert(T, β), C) + if nthreads == 1 || m * n < 1024 + LinearAlgebra.BLAS.gemm!(CA, CB, α, A, B, β, C) + return C + end + + if m > n + m2 = round(Int, m / 16) * 8 + Chead, Ctail = _split_row(C, m2) + Ahead, Atail = CA == 'N' ? _split_row(A, m2) : _split_col(A, m2) + Bhead = Btail = B else - if m > n - m2 = round(Int, m / 16) * 8 - nthreads2 = nthreads >> 1 - t = Threads.@spawn _threaded_blas_mul!( - C[1:($m2), :], A[1:($m2), :], B, α, β, - $nthreads2 - ) - _threaded_blas_mul!( - C[(m2 + 1):m, :], A[(m2 + 1):m, :], B, α, β, - nthreads - nthreads2 - ) - wait(t) - return C - else - n2 = round(Int, n / 16) * 8 - nthreads2 = nthreads >> 1 - t = Threads.@spawn _threaded_blas_mul!( - C[:, 1:($n2)], A, B[:, 1:($n2)], α, β, - $nthreads2 - ) - _threaded_blas_mul!( - C[:, (n2 + 1):n], A, B[:, (n2 + 1):n], α, β, - nthreads - nthreads2 - ) - wait(t) - return C - end + n2 = round(Int, n / 16) * 8 + Chead, Ctail = _split_col(C, n2) + Ahead = Atail = A + Bhead, Btail = CB == 'N' ? _split_col(B, n2) : _split_row(B, n2) end + nthreads1 = nthreads >> 1 + nthreads2 = nthreads - nthreads1 + t = Threads.@spawn _threaded_blas_mul!(Chead, Ahead, CA, Bhead, CB, α, β, nthreads1) + _threaded_blas_mul!(Ctail, Atail, CA, Btail, CB, α, β, nthreads2) + wait(t) + return C end # This implementation is faster than LinearAlgebra.generic_matmatmul From 6a35dec13c52088bc2e65bddb462dc47a22e43a6 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 14 Jul 2026 10:30:52 -0400 Subject: [PATCH 2/3] move generated function to be precompilation-compatible --- src/mapreduce.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/mapreduce.jl b/src/mapreduce.jl index 5de8e6c..7a48290 100644 --- a/src/mapreduce.jl +++ b/src/mapreduce.jl @@ -253,17 +253,6 @@ function _mapreduce_threaded!( return nothing end -@generated function _mapreduce_kernel!( - (f), (op), - (initop), dims::NTuple{N, Int}, - blocks::NTuple{N, Int}, - arrays::NTuple{M, StridedView}, - strides::NTuple{M, NTuple{N, Int}}, - offsets::NTuple{M, Int} - ) where {N, M} - return _mapreduce_kernel_expr(f, op, initop, N, M) -end - function _mapreduce_kernel_expr(f, op, initop, N::Int, M::Int) # many variables @@ -503,6 +492,17 @@ function _mapreduce_kernel_expr(f, op, initop, N::Int, M::Int) return ex end +@generated function _mapreduce_kernel!( + (f), (op), + (initop), dims::NTuple{N, Int}, + blocks::NTuple{N, Int}, + arrays::NTuple{M, StridedView}, + strides::NTuple{M, NTuple{N, Int}}, + offsets::NTuple{M, Int} + ) where {N, M} + return _mapreduce_kernel_expr(f, op, initop, N, M) +end + function indexorder(strides::NTuple{N, Int}) where {N} # returns order such that strides[i] is the order[i]th smallest element of strides, not # counting zero strides zero strides have order 1 From 6027a6e73d35a6d6103396199f1b9df1b181781e Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 14 Jul 2026 10:31:29 -0400 Subject: [PATCH 3/3] add precompilaton --- Project.toml | 2 ++ src/Strided.jl | 2 ++ src/precompile.jl | 23 +++++++++++++++++++++++ 3 files changed, 27 insertions(+) create mode 100644 src/precompile.jl diff --git a/Project.toml b/Project.toml index da9069b..b76a0b4 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ authors = ["Lukas Devos ", "Maarten Van Damme