From 5a97dbb36d93527cf076a0182c995d3bc1de6aa3 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 20 Aug 2025 09:04:47 -0400 Subject: [PATCH 1/7] Add OpenBLASLUFactorization implementation - Implement OpenBLASLUFactorization as a direct wrapper over OpenBLAS_jll - Add getrf! and getrs! functions for LU factorization and solving - Support Float32, Float64, ComplexF32, and ComplexF64 types - Include proper module structure and exports - Add OpenBLAS_jll as a dependency - Tests confirm functionality matches existing LUFactorization --- Project.toml | 1 + src/LinearSolve.jl | 4 + src/openblas.jl | 252 +++++++++++++++++++++++++++++++++++++++++++++ test_openblas.jl | 59 +++++++++++ 4 files changed, 316 insertions(+) create mode 100644 src/openblas.jl create mode 100644 test_openblas.jl diff --git a/Project.toml b/Project.toml index ad009a58c..482834e32 100644 --- a/Project.toml +++ b/Project.toml @@ -17,6 +17,7 @@ Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MKL_jll = "856f044c-d86e-5d09-b602-aeab76dc8ba7" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" +OpenBLAS_jll = "4536629a-c528-5b80-bd46-f80d51c5b363" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Preferences = "21216c6a-2e73-6563-6e65-726566657250" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index e8bcc5b3e..098222689 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -59,6 +59,8 @@ else const usemkl = false end +using OpenBLAS_jll + @reexport using SciMLBase @@ -345,6 +347,7 @@ include("extension_algs.jl") include("factorization.jl") include("appleaccelerate.jl") include("mkl.jl") +include("openblas.jl") include("simplelu.jl") include("simplegmres.jl") include("iterative_wrappers.jl") @@ -461,6 +464,7 @@ export MKLPardisoFactorize, MKLPardisoIterate export PanuaPardisoFactorize, PanuaPardisoIterate export PardisoJL export MKLLUFactorization +export OpenBLASLUFactorization export AppleAccelerateLUFactorization export MetalLUFactorization diff --git a/src/openblas.jl b/src/openblas.jl new file mode 100644 index 000000000..6381507a9 --- /dev/null +++ b/src/openblas.jl @@ -0,0 +1,252 @@ +""" +```julia +OpenBLASLUFactorization() +``` + +A wrapper over OpenBLAS. Direct calls to OpenBLAS in a way that pre-allocates workspace +to avoid allocations and does not require libblastrampoline. +""" +struct OpenBLASLUFactorization <: AbstractFactorization end + +module OpenBLASLU + +using LinearAlgebra +using LinearAlgebra: BlasInt, LU, require_one_based_indexing, checksquare +using LinearAlgebra.LAPACK: chkfinite, chkstride1, @blasfunc, chkargsok, chktrans, chklapackerror +using OpenBLAS_jll + +function getrf!(A::AbstractMatrix{<:ComplexF64}; + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{BlasInt}(), + check = false) + require_one_based_indexing(A) + check && chkfinite(A) + chkstride1(A) + m, n = size(A) + lda = max(1, stride(A, 2)) + if isempty(ipiv) + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))) + end + ccall((@blasfunc(zgetrf_), OpenBLAS_jll.libopenblas), Cvoid, + (Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF64}, + Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), + m, n, A, lda, ipiv, info) + chkargsok(info[]) + A, ipiv, info[], info #Error code is stored in LU factorization type +end + +function getrf!(A::AbstractMatrix{<:ComplexF32}; + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{BlasInt}(), + check = false) + require_one_based_indexing(A) + check && chkfinite(A) + chkstride1(A) + m, n = size(A) + lda = max(1, stride(A, 2)) + if isempty(ipiv) + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))) + end + ccall((@blasfunc(cgetrf_), OpenBLAS_jll.libopenblas), Cvoid, + (Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF32}, + Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), + m, n, A, lda, ipiv, info) + chkargsok(info[]) + A, ipiv, info[], info #Error code is stored in LU factorization type +end + +function getrf!(A::AbstractMatrix{<:Float64}; + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{BlasInt}(), + check = false) + require_one_based_indexing(A) + check && chkfinite(A) + chkstride1(A) + m, n = size(A) + lda = max(1, stride(A, 2)) + if isempty(ipiv) + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))) + end + ccall((@blasfunc(dgetrf_), OpenBLAS_jll.libopenblas), Cvoid, + (Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, + Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), + m, n, A, lda, ipiv, info) + chkargsok(info[]) + A, ipiv, info[], info #Error code is stored in LU factorization type +end + +function getrf!(A::AbstractMatrix{<:Float32}; + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{BlasInt}(), + check = false) + require_one_based_indexing(A) + check && chkfinite(A) + chkstride1(A) + m, n = size(A) + lda = max(1, stride(A, 2)) + if isempty(ipiv) + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))) + end + ccall((@blasfunc(sgetrf_), OpenBLAS_jll.libopenblas), Cvoid, + (Ref{BlasInt}, Ref{BlasInt}, Ptr{Float32}, + Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), + m, n, A, lda, ipiv, info) + chkargsok(info[]) + A, ipiv, info[], info #Error code is stored in LU factorization type +end + +function getrs!(trans::AbstractChar, + A::AbstractMatrix{<:ComplexF64}, + ipiv::AbstractVector{BlasInt}, + B::AbstractVecOrMat{<:ComplexF64}; + info = Ref{BlasInt}()) + require_one_based_indexing(A, ipiv, B) + LinearAlgebra.LAPACK.chktrans(trans) + chkstride1(A, B, ipiv) + n = LinearAlgebra.checksquare(A) + if n != size(B, 1) + throw(DimensionMismatch("B has leading dimension $(size(B,1)), but needs $n")) + end + if n != length(ipiv) + throw(DimensionMismatch("ipiv has length $(length(ipiv)), but needs to be $n")) + end + nrhs = size(B, 2) + ccall((@blasfunc(zgetrs_), OpenBLAS_jll.libopenblas), Cvoid, + (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF64}, Ref{BlasInt}, + Ptr{BlasInt}, Ptr{ComplexF64}, Ref{BlasInt}, Ptr{BlasInt}, Clong), + trans, n, size(B, 2), A, max(1, stride(A, 2)), ipiv, B, max(1, stride(B, 2)), info, + 1) + LinearAlgebra.LAPACK.chklapackerror(BlasInt(info[])) + B +end + +function getrs!(trans::AbstractChar, + A::AbstractMatrix{<:ComplexF32}, + ipiv::AbstractVector{BlasInt}, + B::AbstractVecOrMat{<:ComplexF32}; + info = Ref{BlasInt}()) + require_one_based_indexing(A, ipiv, B) + LinearAlgebra.LAPACK.chktrans(trans) + chkstride1(A, B, ipiv) + n = LinearAlgebra.checksquare(A) + if n != size(B, 1) + throw(DimensionMismatch("B has leading dimension $(size(B,1)), but needs $n")) + end + if n != length(ipiv) + throw(DimensionMismatch("ipiv has length $(length(ipiv)), but needs to be $n")) + end + nrhs = size(B, 2) + ccall((@blasfunc(cgetrs_), OpenBLAS_jll.libopenblas), Cvoid, + (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF32}, Ref{BlasInt}, + Ptr{BlasInt}, Ptr{ComplexF32}, Ref{BlasInt}, Ptr{BlasInt}, Clong), + trans, n, size(B, 2), A, max(1, stride(A, 2)), ipiv, B, max(1, stride(B, 2)), info, + 1) + LinearAlgebra.LAPACK.chklapackerror(BlasInt(info[])) + B +end + +function getrs!(trans::AbstractChar, + A::AbstractMatrix{<:Float64}, + ipiv::AbstractVector{BlasInt}, + B::AbstractVecOrMat{<:Float64}; + info = Ref{BlasInt}()) + require_one_based_indexing(A, ipiv, B) + LinearAlgebra.LAPACK.chktrans(trans) + chkstride1(A, B, ipiv) + n = LinearAlgebra.checksquare(A) + if n != size(B, 1) + throw(DimensionMismatch("B has leading dimension $(size(B,1)), but needs $n")) + end + if n != length(ipiv) + throw(DimensionMismatch("ipiv has length $(length(ipiv)), but needs to be $n")) + end + nrhs = size(B, 2) + ccall((@blasfunc(dgetrs_), OpenBLAS_jll.libopenblas), Cvoid, + (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, + Ptr{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{BlasInt}, Clong), + trans, n, size(B, 2), A, max(1, stride(A, 2)), ipiv, B, max(1, stride(B, 2)), info, + 1) + LinearAlgebra.LAPACK.chklapackerror(BlasInt(info[])) + B +end + +function getrs!(trans::AbstractChar, + A::AbstractMatrix{<:Float32}, + ipiv::AbstractVector{BlasInt}, + B::AbstractVecOrMat{<:Float32}; + info = Ref{BlasInt}()) + require_one_based_indexing(A, ipiv, B) + LinearAlgebra.LAPACK.chktrans(trans) + chkstride1(A, B, ipiv) + n = LinearAlgebra.checksquare(A) + if n != size(B, 1) + throw(DimensionMismatch("B has leading dimension $(size(B,1)), but needs $n")) + end + if n != length(ipiv) + throw(DimensionMismatch("ipiv has length $(length(ipiv)), but needs to be $n")) + end + nrhs = size(B, 2) + ccall((@blasfunc(sgetrs_), OpenBLAS_jll.libopenblas), Cvoid, + (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{Float32}, Ref{BlasInt}, + Ptr{BlasInt}, Ptr{Float32}, Ref{BlasInt}, Ptr{BlasInt}, Clong), + trans, n, size(B, 2), A, max(1, stride(A, 2)), ipiv, B, max(1, stride(B, 2)), info, + 1) + LinearAlgebra.LAPACK.chklapackerror(BlasInt(info[])) + B +end + +end # module OpenBLASLU + +default_alias_A(::OpenBLASLUFactorization, ::Any, ::Any) = false +default_alias_b(::OpenBLASLUFactorization, ::Any, ::Any) = false + +const PREALLOCATED_OPENBLAS_LU = begin + A = rand(0, 0) + luinst = ArrayInterface.lu_instance(A), Ref{BlasInt}() +end + +function LinearSolve.init_cacheval(alg::OpenBLASLUFactorization, A, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::LinearVerbosity, + assumptions::OperatorAssumptions) + PREALLOCATED_OPENBLAS_LU +end + +function LinearSolve.init_cacheval(alg::OpenBLASLUFactorization, + A::AbstractMatrix{<:Union{Float32, ComplexF32, ComplexF64}}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::LinearVerbosity, + assumptions::OperatorAssumptions) + A = rand(eltype(A), 0, 0) + ArrayInterface.lu_instance(A), Ref{BlasInt}() +end + +function SciMLBase.solve!(cache::LinearCache, alg::OpenBLASLUFactorization; + kwargs...) + A = cache.A + A = convert(AbstractMatrix, A) + if cache.isfresh + cacheval = @get_cacheval(cache, :OpenBLASLUFactorization) + res = OpenBLASLU.getrf!(A; ipiv = cacheval[1].ipiv, info = cacheval[2]) + fact = LU(res[1:3]...), res[4] + cache.cacheval = fact + + if !LinearAlgebra.issuccess(fact[1]) + return SciMLBase.build_linear_solution( + alg, cache.u, nothing, cache; retcode = ReturnCode.Failure) + end + cache.isfresh = false + end + + A, info = @get_cacheval(cache, :OpenBLASLUFactorization) + require_one_based_indexing(cache.u, cache.b) + m, n = size(A, 1), size(A, 2) + if m > n + Bc = copy(cache.b) + OpenBLASLU.getrs!('N', A.factors, A.ipiv, Bc; info) + copyto!(cache.u, 1, Bc, 1, n) + else + copyto!(cache.u, cache.b) + OpenBLASLU.getrs!('N', A.factors, A.ipiv, cache.u; info) + end + + SciMLBase.build_linear_solution(alg, cache.u, nothing, cache; retcode = ReturnCode.Success) +end \ No newline at end of file diff --git a/test_openblas.jl b/test_openblas.jl new file mode 100644 index 000000000..9d31a35c5 --- /dev/null +++ b/test_openblas.jl @@ -0,0 +1,59 @@ +using LinearAlgebra +using LinearSolve +using Test + +@testset "OpenBLASLUFactorization Tests" begin + # Test with Float64 + @testset "Float64" begin + A = rand(10, 10) + b = rand(10) + prob = LinearProblem(A, b) + + sol_openblas = solve(prob, OpenBLASLUFactorization()) + sol_default = solve(prob, LUFactorization()) + + @test norm(A * sol_openblas.u - b) < 1e-10 + @test norm(sol_openblas.u - sol_default.u) < 1e-10 + end + + # Test with Float32 + @testset "Float32" begin + A = rand(Float32, 10, 10) + b = rand(Float32, 10) + prob = LinearProblem(A, b) + + sol_openblas = solve(prob, OpenBLASLUFactorization()) + sol_default = solve(prob, LUFactorization()) + + @test norm(A * sol_openblas.u - b) < 1e-5 + @test norm(sol_openblas.u - sol_default.u) < 1e-5 + end + + # Test with ComplexF64 + @testset "ComplexF64" begin + A = rand(ComplexF64, 10, 10) + b = rand(ComplexF64, 10) + prob = LinearProblem(A, b) + + sol_openblas = solve(prob, OpenBLASLUFactorization()) + sol_default = solve(prob, LUFactorization()) + + @test norm(A * sol_openblas.u - b) < 1e-10 + @test norm(sol_openblas.u - sol_default.u) < 1e-10 + end + + # Test with ComplexF32 + @testset "ComplexF32" begin + A = rand(ComplexF32, 10, 10) + b = rand(ComplexF32, 10) + prob = LinearProblem(A, b) + + sol_openblas = solve(prob, OpenBLASLUFactorization()) + sol_default = solve(prob, LUFactorization()) + + @test norm(A * sol_openblas.u - b) < 1e-5 + @test norm(sol_openblas.u - sol_default.u) < 1e-5 + end +end + +println("All tests passed!") \ No newline at end of file From 4b7ef0a66f517e146f5ec7726a7bc2c29729795e Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 20 Aug 2025 09:41:24 -0400 Subject: [PATCH 2/7] Integrate OpenBLASLUFactorization into existing test suite - Add OpenBLASLUFactorization to basictests.jl alongside other factorizations - Include in resolve.jl subtype iteration tests - Add to preferences.jl algorithm availability tests - Remove separate test file in favor of integrated testing --- test/basictests.jl | 3 +++ test/preferences.jl | 6 +++++ test/resolve.jl | 3 ++- test_openblas.jl | 59 --------------------------------------------- 4 files changed, 11 insertions(+), 60 deletions(-) delete mode 100644 test_openblas.jl diff --git a/test/basictests.jl b/test/basictests.jl index 8aa7547a9..f360d87a1 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -287,6 +287,9 @@ end push!(test_algs, MKLLUFactorization()) end + # Always test OpenBLAS since it's a direct dependency + push!(test_algs, OpenBLASLUFactorization()) + # Test BLIS if extension is available if Base.get_extension(LinearSolve, :LinearSolveBLISExt) !== nothing push!(test_algs, BLISLUFactorization()) diff --git a/test/preferences.jl b/test/preferences.jl index bd6cf4516..0182657b4 100644 --- a/test/preferences.jl +++ b/test/preferences.jl @@ -190,6 +190,12 @@ using Preferences println("✅ MKLLUFactorization confirmed working") end + # Test OpenBLAS (always available as a dependency) + sol_openblas = solve(prob, OpenBLASLUFactorization()) + @test sol_openblas.retcode == ReturnCode.Success + @test norm(A * sol_openblas.u - b) < 1e-8 + println("✅ OpenBLASLUFactorization confirmed working") + # Test Apple Accelerate if available if LinearSolve.appleaccelerate_isavailable() sol_apple = solve(prob, AppleAccelerateLUFactorization()) diff --git a/test/resolve.jl b/test/resolve.jl index 3fe0a2fea..0f8bdbea0 100644 --- a/test/resolve.jl +++ b/test/resolve.jl @@ -19,7 +19,8 @@ for alg in vcat(InteractiveUtils.subtypes(AbstractDenseFactorization), ]) && (!(alg == AppleAccelerateLUFactorization) || LinearSolve.appleaccelerate_isavailable()) && - (!(alg == MKLLUFactorization) || LinearSolve.usemkl) + (!(alg == MKLLUFactorization) || LinearSolve.usemkl) && + (!(alg == OpenBLASLUFactorization) || true) # OpenBLAS is always available as a dependency A = [1.0 2.0; 3.0 4.0] alg in [KLUFactorization, UMFPACKFactorization, SparspakFactorization] && (A = sparse(A)) diff --git a/test_openblas.jl b/test_openblas.jl deleted file mode 100644 index 9d31a35c5..000000000 --- a/test_openblas.jl +++ /dev/null @@ -1,59 +0,0 @@ -using LinearAlgebra -using LinearSolve -using Test - -@testset "OpenBLASLUFactorization Tests" begin - # Test with Float64 - @testset "Float64" begin - A = rand(10, 10) - b = rand(10) - prob = LinearProblem(A, b) - - sol_openblas = solve(prob, OpenBLASLUFactorization()) - sol_default = solve(prob, LUFactorization()) - - @test norm(A * sol_openblas.u - b) < 1e-10 - @test norm(sol_openblas.u - sol_default.u) < 1e-10 - end - - # Test with Float32 - @testset "Float32" begin - A = rand(Float32, 10, 10) - b = rand(Float32, 10) - prob = LinearProblem(A, b) - - sol_openblas = solve(prob, OpenBLASLUFactorization()) - sol_default = solve(prob, LUFactorization()) - - @test norm(A * sol_openblas.u - b) < 1e-5 - @test norm(sol_openblas.u - sol_default.u) < 1e-5 - end - - # Test with ComplexF64 - @testset "ComplexF64" begin - A = rand(ComplexF64, 10, 10) - b = rand(ComplexF64, 10) - prob = LinearProblem(A, b) - - sol_openblas = solve(prob, OpenBLASLUFactorization()) - sol_default = solve(prob, LUFactorization()) - - @test norm(A * sol_openblas.u - b) < 1e-10 - @test norm(sol_openblas.u - sol_default.u) < 1e-10 - end - - # Test with ComplexF32 - @testset "ComplexF32" begin - A = rand(ComplexF32, 10, 10) - b = rand(ComplexF32, 10) - prob = LinearProblem(A, b) - - sol_openblas = solve(prob, OpenBLASLUFactorization()) - sol_default = solve(prob, LUFactorization()) - - @test norm(A * sol_openblas.u - b) < 1e-5 - @test norm(sol_openblas.u - sol_default.u) < 1e-5 - end -end - -println("All tests passed!") \ No newline at end of file From 676c27f56fb7e856d6cc43b5b000c638a5509459 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 20 Aug 2025 12:57:33 -0400 Subject: [PATCH 3/7] Add comprehensive documentation for OpenBLASLUFactorization - Add detailed docstring with performance characteristics and usage examples - Include OpenBLASLUFactorization in solver documentation alongside MKL and AppleAccelerate - Update algorithm selection guide to mention OpenBLAS as an option for large dense matrices - Document when to use OpenBLAS vs other BLAS implementations --- docs/src/basics/algorithm_selection.md | 5 +++-- docs/src/solvers/solvers.md | 14 ++++++++++--- src/openblas.jl | 29 ++++++++++++++++++++++++-- 3 files changed, 41 insertions(+), 7 deletions(-) diff --git a/docs/src/basics/algorithm_selection.md b/docs/src/basics/algorithm_selection.md index f7ae4feda..37058c4ef 100644 --- a/docs/src/basics/algorithm_selection.md +++ b/docs/src/basics/algorithm_selection.md @@ -71,9 +71,10 @@ sol = solve(LinearProblem(A_small, rand(50)), SimpleLUFactorization()) A_medium = rand(200, 200) sol = solve(LinearProblem(A_medium, rand(200)), RFLUFactorization()) -# Large matrices (> 500×500): MKLLUFactorization or AppleAccelerate +# Large matrices (> 500×500): MKLLUFactorization, OpenBLASLUFactorization, or AppleAccelerate A_large = rand(1000, 1000) sol = solve(LinearProblem(A_large, rand(1000)), MKLLUFactorization()) +# Alternative: OpenBLASLUFactorization() for direct OpenBLAS calls ``` ### Sparse Matrices @@ -141,7 +142,7 @@ Is A symmetric positive definite? → CholeskyFactorization Is A symmetric indefinite? → BunchKaufmanFactorization Is A sparse? → UMFPACKFactorization or KLUFactorization Is A small dense? → RFLUFactorization or SimpleLUFactorization -Is A large dense? → MKLLUFactorization or AppleAccelerateLUFactorization +Is A large dense? → MKLLUFactorization, OpenBLASLUFactorization, or AppleAccelerateLUFactorization Is A GPU array? → QRFactorization or LUFactorization Is A an operator/function? → KrylovJL_GMRES Is the system overdetermined? → QRFactorization or KrylovJL_LSMR diff --git a/docs/src/solvers/solvers.md b/docs/src/solvers/solvers.md index 1290a8fa2..e10a4d7a9 100644 --- a/docs/src/solvers/solvers.md +++ b/docs/src/solvers/solvers.md @@ -17,9 +17,11 @@ the best choices, with SVD being the slowest but most precise. For efficiency, `RFLUFactorization` is the fastest for dense LU-factorizations until around 150x150 matrices, though this can be dependent on the exact details of the hardware. After this point, `MKLLUFactorization` is usually faster on most hardware. Note that on Mac computers -that `AppleAccelerateLUFactorization` is generally always the fastest. `LUFactorization` will -use your base system BLAS which can be fast or slow depending on the hardware configuration. -`SimpleLUFactorization` will be fast only on very small matrices but can cut down on compile times. +that `AppleAccelerateLUFactorization` is generally always the fastest. `OpenBLASLUFactorization` +provides direct OpenBLAS calls without going through libblastrampoline and can be faster than +`LUFactorization` in some configurations. `LUFactorization` will use your base system BLAS which +can be fast or slow depending on the hardware configuration. `SimpleLUFactorization` will be fast +only on very small matrices but can cut down on compile times. For very large dense factorizations, offloading to the GPU can be preferred. Metal.jl can be used on Mac hardware to offload, and has a cutoff point of being faster at around size 20,000 x 20,000 @@ -207,6 +209,12 @@ KrylovJL MKLLUFactorization ``` +### OpenBLAS + +```@docs +OpenBLASLUFactorization +``` + ### AppleAccelerate.jl !!! note diff --git a/src/openblas.jl b/src/openblas.jl index 6381507a9..fe5b2f3e4 100644 --- a/src/openblas.jl +++ b/src/openblas.jl @@ -3,8 +3,33 @@ OpenBLASLUFactorization() ``` -A wrapper over OpenBLAS. Direct calls to OpenBLAS in a way that pre-allocates workspace -to avoid allocations and does not require libblastrampoline. +A direct wrapper over OpenBLAS's LU factorization (`getrf!` and `getrs!`). +This solver makes direct calls to OpenBLAS_jll without going through Julia's +libblastrampoline, which can provide performance benefits in certain configurations. + +## Performance Characteristics + +- Pre-allocates workspace to avoid allocations during solving +- Makes direct `ccall`s to OpenBLAS routines +- Can be faster than `LUFactorization` when OpenBLAS is well-optimized for the hardware +- Supports `Float32`, `Float64`, `ComplexF32`, and `ComplexF64` element types + +## When to Use + +- When you want to ensure OpenBLAS is used regardless of the system BLAS configuration +- When benchmarking shows better performance than `LUFactorization` on your specific hardware +- When you need consistent behavior across different systems (always uses OpenBLAS) + +## Example + +```julia +using LinearSolve, LinearAlgebra + +A = rand(100, 100) +b = rand(100) +prob = LinearProblem(A, b) +sol = solve(prob, OpenBLASLUFactorization()) +``` """ struct OpenBLASLUFactorization <: AbstractFactorization end From d4da31b4f93d9da8a851626333795752f95d802e Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 20 Aug 2025 15:33:22 -0400 Subject: [PATCH 4/7] Fix stale imports and add compat entry for OpenBLAS_jll - Remove redundant import of LinearAlgebra items that are already available via using LinearAlgebra - Qualify all uses of BlasInt, LU, require_one_based_indexing, and checksquare with LinearAlgebra prefix - Add OpenBLAS_jll = "0.3" to compat section in Project.toml - Apply JuliaFormatter with SciMLStyle to ensure consistent formatting --- Project.toml | 1 + src/openblas.jl | 137 +++++++++++++++++++++++++----------------------- 2 files changed, 73 insertions(+), 65 deletions(-) diff --git a/Project.toml b/Project.toml index 482834e32..c9759eeb9 100644 --- a/Project.toml +++ b/Project.toml @@ -113,6 +113,7 @@ MPI = "0.20" Markdown = "1.10" Metal = "1.4" MultiFloats = "2.3" +OpenBLAS_jll = "0.3" Pardiso = "1" Pkg = "1.10" PrecompileTools = "1.2" diff --git a/src/openblas.jl b/src/openblas.jl index fe5b2f3e4..f5e7fbc63 100644 --- a/src/openblas.jl +++ b/src/openblas.jl @@ -3,22 +3,22 @@ OpenBLASLUFactorization() ``` -A direct wrapper over OpenBLAS's LU factorization (`getrf!` and `getrs!`). -This solver makes direct calls to OpenBLAS_jll without going through Julia's +A direct wrapper over OpenBLAS's LU factorization (`getrf!` and `getrs!`). +This solver makes direct calls to OpenBLAS_jll without going through Julia's libblastrampoline, which can provide performance benefits in certain configurations. ## Performance Characteristics -- Pre-allocates workspace to avoid allocations during solving -- Makes direct `ccall`s to OpenBLAS routines -- Can be faster than `LUFactorization` when OpenBLAS is well-optimized for the hardware -- Supports `Float32`, `Float64`, `ComplexF32`, and `ComplexF64` element types + - Pre-allocates workspace to avoid allocations during solving + - Makes direct `ccall`s to OpenBLAS routines + - Can be faster than `LUFactorization` when OpenBLAS is well-optimized for the hardware + - Supports `Float32`, `Float64`, `ComplexF32`, and `ComplexF64` element types ## When to Use -- When you want to ensure OpenBLAS is used regardless of the system BLAS configuration -- When benchmarking shows better performance than `LUFactorization` on your specific hardware -- When you need consistent behavior across different systems (always uses OpenBLAS) + - When you want to ensure OpenBLAS is used regardless of the system BLAS configuration + - When benchmarking shows better performance than `LUFactorization` on your specific hardware + - When you need consistent behavior across different systems (always uses OpenBLAS) ## Example @@ -36,85 +36,85 @@ struct OpenBLASLUFactorization <: AbstractFactorization end module OpenBLASLU using LinearAlgebra -using LinearAlgebra: BlasInt, LU, require_one_based_indexing, checksquare -using LinearAlgebra.LAPACK: chkfinite, chkstride1, @blasfunc, chkargsok, chktrans, chklapackerror +using LinearAlgebra.LAPACK: chkfinite, chkstride1, @blasfunc, chkargsok, chktrans, + chklapackerror using OpenBLAS_jll function getrf!(A::AbstractMatrix{<:ComplexF64}; - ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), - info = Ref{BlasInt}(), + ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{LinearAlgebra.BlasInt}(), check = false) - require_one_based_indexing(A) + LinearAlgebra.require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) m, n = size(A) lda = max(1, stride(A, 2)) if isempty(ipiv) - ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))) + ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))) end ccall((@blasfunc(zgetrf_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF64}, - Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), + (Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ptr{ComplexF64}, + Ref{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}), m, n, A, lda, ipiv, info) chkargsok(info[]) A, ipiv, info[], info #Error code is stored in LU factorization type end function getrf!(A::AbstractMatrix{<:ComplexF32}; - ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), - info = Ref{BlasInt}(), + ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{LinearAlgebra.BlasInt}(), check = false) - require_one_based_indexing(A) + LinearAlgebra.require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) m, n = size(A) lda = max(1, stride(A, 2)) if isempty(ipiv) - ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))) + ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))) end ccall((@blasfunc(cgetrf_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF32}, - Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), + (Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ptr{ComplexF32}, + Ref{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}), m, n, A, lda, ipiv, info) chkargsok(info[]) A, ipiv, info[], info #Error code is stored in LU factorization type end function getrf!(A::AbstractMatrix{<:Float64}; - ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), - info = Ref{BlasInt}(), + ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{LinearAlgebra.BlasInt}(), check = false) - require_one_based_indexing(A) + LinearAlgebra.require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) m, n = size(A) lda = max(1, stride(A, 2)) if isempty(ipiv) - ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))) + ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))) end ccall((@blasfunc(dgetrf_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, - Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), + (Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ptr{Float64}, + Ref{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}), m, n, A, lda, ipiv, info) chkargsok(info[]) A, ipiv, info[], info #Error code is stored in LU factorization type end function getrf!(A::AbstractMatrix{<:Float32}; - ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), - info = Ref{BlasInt}(), + ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{LinearAlgebra.BlasInt}(), check = false) - require_one_based_indexing(A) + LinearAlgebra.require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) m, n = size(A) lda = max(1, stride(A, 2)) if isempty(ipiv) - ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))) + ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))) end ccall((@blasfunc(sgetrf_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{BlasInt}, Ref{BlasInt}, Ptr{Float32}, - Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), + (Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ptr{Float32}, + Ref{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}), m, n, A, lda, ipiv, info) chkargsok(info[]) A, ipiv, info[], info #Error code is stored in LU factorization type @@ -122,10 +122,10 @@ end function getrs!(trans::AbstractChar, A::AbstractMatrix{<:ComplexF64}, - ipiv::AbstractVector{BlasInt}, + ipiv::AbstractVector{LinearAlgebra.BlasInt}, B::AbstractVecOrMat{<:ComplexF64}; - info = Ref{BlasInt}()) - require_one_based_indexing(A, ipiv, B) + info = Ref{LinearAlgebra.BlasInt}()) + LinearAlgebra.require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) n = LinearAlgebra.checksquare(A) @@ -137,20 +137,22 @@ function getrs!(trans::AbstractChar, end nrhs = size(B, 2) ccall((@blasfunc(zgetrs_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF64}, Ref{BlasInt}, - Ptr{BlasInt}, Ptr{ComplexF64}, Ref{BlasInt}, Ptr{BlasInt}, Clong), + (Ref{UInt8}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, + Ptr{ComplexF64}, Ref{LinearAlgebra.BlasInt}, + Ptr{LinearAlgebra.BlasInt}, Ptr{ComplexF64}, Ref{LinearAlgebra.BlasInt}, + Ptr{LinearAlgebra.BlasInt}, Clong), trans, n, size(B, 2), A, max(1, stride(A, 2)), ipiv, B, max(1, stride(B, 2)), info, 1) - LinearAlgebra.LAPACK.chklapackerror(BlasInt(info[])) + LinearAlgebra.LAPACK.chklapackerror(LinearAlgebra.BlasInt(info[])) B end function getrs!(trans::AbstractChar, A::AbstractMatrix{<:ComplexF32}, - ipiv::AbstractVector{BlasInt}, + ipiv::AbstractVector{LinearAlgebra.BlasInt}, B::AbstractVecOrMat{<:ComplexF32}; - info = Ref{BlasInt}()) - require_one_based_indexing(A, ipiv, B) + info = Ref{LinearAlgebra.BlasInt}()) + LinearAlgebra.require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) n = LinearAlgebra.checksquare(A) @@ -162,20 +164,22 @@ function getrs!(trans::AbstractChar, end nrhs = size(B, 2) ccall((@blasfunc(cgetrs_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF32}, Ref{BlasInt}, - Ptr{BlasInt}, Ptr{ComplexF32}, Ref{BlasInt}, Ptr{BlasInt}, Clong), + (Ref{UInt8}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, + Ptr{ComplexF32}, Ref{LinearAlgebra.BlasInt}, + Ptr{LinearAlgebra.BlasInt}, Ptr{ComplexF32}, Ref{LinearAlgebra.BlasInt}, + Ptr{LinearAlgebra.BlasInt}, Clong), trans, n, size(B, 2), A, max(1, stride(A, 2)), ipiv, B, max(1, stride(B, 2)), info, 1) - LinearAlgebra.LAPACK.chklapackerror(BlasInt(info[])) + LinearAlgebra.LAPACK.chklapackerror(LinearAlgebra.BlasInt(info[])) B end function getrs!(trans::AbstractChar, A::AbstractMatrix{<:Float64}, - ipiv::AbstractVector{BlasInt}, + ipiv::AbstractVector{LinearAlgebra.BlasInt}, B::AbstractVecOrMat{<:Float64}; - info = Ref{BlasInt}()) - require_one_based_indexing(A, ipiv, B) + info = Ref{LinearAlgebra.BlasInt}()) + LinearAlgebra.require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) n = LinearAlgebra.checksquare(A) @@ -187,20 +191,21 @@ function getrs!(trans::AbstractChar, end nrhs = size(B, 2) ccall((@blasfunc(dgetrs_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, - Ptr{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{BlasInt}, Clong), + (Ref{UInt8}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, + Ptr{Float64}, Ref{LinearAlgebra.BlasInt}, + Ptr{LinearAlgebra.BlasInt}, Ptr{Float64}, Ref{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}, Clong), trans, n, size(B, 2), A, max(1, stride(A, 2)), ipiv, B, max(1, stride(B, 2)), info, 1) - LinearAlgebra.LAPACK.chklapackerror(BlasInt(info[])) + LinearAlgebra.LAPACK.chklapackerror(LinearAlgebra.BlasInt(info[])) B end function getrs!(trans::AbstractChar, A::AbstractMatrix{<:Float32}, - ipiv::AbstractVector{BlasInt}, + ipiv::AbstractVector{LinearAlgebra.BlasInt}, B::AbstractVecOrMat{<:Float32}; - info = Ref{BlasInt}()) - require_one_based_indexing(A, ipiv, B) + info = Ref{LinearAlgebra.BlasInt}()) + LinearAlgebra.require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) n = LinearAlgebra.checksquare(A) @@ -212,11 +217,12 @@ function getrs!(trans::AbstractChar, end nrhs = size(B, 2) ccall((@blasfunc(sgetrs_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{Float32}, Ref{BlasInt}, - Ptr{BlasInt}, Ptr{Float32}, Ref{BlasInt}, Ptr{BlasInt}, Clong), + (Ref{UInt8}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, + Ptr{Float32}, Ref{LinearAlgebra.BlasInt}, + Ptr{LinearAlgebra.BlasInt}, Ptr{Float32}, Ref{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}, Clong), trans, n, size(B, 2), A, max(1, stride(A, 2)), ipiv, B, max(1, stride(B, 2)), info, 1) - LinearAlgebra.LAPACK.chklapackerror(BlasInt(info[])) + LinearAlgebra.LAPACK.chklapackerror(LinearAlgebra.BlasInt(info[])) B end @@ -227,7 +233,7 @@ default_alias_b(::OpenBLASLUFactorization, ::Any, ::Any) = false const PREALLOCATED_OPENBLAS_LU = begin A = rand(0, 0) - luinst = ArrayInterface.lu_instance(A), Ref{BlasInt}() + luinst = ArrayInterface.lu_instance(A), Ref{LinearAlgebra.BlasInt}() end function LinearSolve.init_cacheval(alg::OpenBLASLUFactorization, A, b, u, Pl, Pr, @@ -241,7 +247,7 @@ function LinearSolve.init_cacheval(alg::OpenBLASLUFactorization, maxiters::Int, abstol, reltol, verbose::LinearVerbosity, assumptions::OperatorAssumptions) A = rand(eltype(A), 0, 0) - ArrayInterface.lu_instance(A), Ref{BlasInt}() + ArrayInterface.lu_instance(A), Ref{LinearAlgebra.BlasInt}() end function SciMLBase.solve!(cache::LinearCache, alg::OpenBLASLUFactorization; @@ -251,7 +257,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::OpenBLASLUFactorization; if cache.isfresh cacheval = @get_cacheval(cache, :OpenBLASLUFactorization) res = OpenBLASLU.getrf!(A; ipiv = cacheval[1].ipiv, info = cacheval[2]) - fact = LU(res[1:3]...), res[4] + fact = LinearAlgebra.LU(res[1:3]...), res[4] cache.cacheval = fact if !LinearAlgebra.issuccess(fact[1]) @@ -262,7 +268,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::OpenBLASLUFactorization; end A, info = @get_cacheval(cache, :OpenBLASLUFactorization) - require_one_based_indexing(cache.u, cache.b) + LinearAlgebra.require_one_based_indexing(cache.u, cache.b) m, n = size(A, 1), size(A, 2) if m > n Bc = copy(cache.b) @@ -273,5 +279,6 @@ function SciMLBase.solve!(cache::LinearCache, alg::OpenBLASLUFactorization; OpenBLASLU.getrs!('N', A.factors, A.ipiv, cache.u; info) end - SciMLBase.build_linear_solution(alg, cache.u, nothing, cache; retcode = ReturnCode.Success) -end \ No newline at end of file + SciMLBase.build_linear_solution( + alg, cache.u, nothing, cache; retcode = ReturnCode.Success) +end From a89490a365e46646c649181236893e0f9cbf9138 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 20 Aug 2025 16:28:13 -0400 Subject: [PATCH 5/7] Refactor OpenBLAS implementation to match MKL pattern - Remove separate module structure, follow MKL's pattern exactly - Conditionally check for OpenBLAS_jll availability using is_available() - Keep OpenBLAS_jll as a dependency (required even for stdlib packages) - Simplify implementation without nested @static checks - Tests conditionally run based on LinearSolve.useopenblas flag --- src/LinearSolve.jl | 2 + src/openblas.jl | 141 ++++++++++++++++++++------------------------ test/basictests.jl | 6 +- test/preferences.jl | 12 ++-- test/resolve.jl | 2 +- 5 files changed, 78 insertions(+), 85 deletions(-) diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index 25b19bbf0..10ae83894 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -59,7 +59,9 @@ else const usemkl = false end +# OpenBLAS_jll is a standard library, always available using OpenBLAS_jll +const useopenblas = OpenBLAS_jll.is_available() @reexport using SciMLBase diff --git a/src/openblas.jl b/src/openblas.jl index f5e7fbc63..ea6cb4210 100644 --- a/src/openblas.jl +++ b/src/openblas.jl @@ -33,99 +33,94 @@ sol = solve(prob, OpenBLASLUFactorization()) """ struct OpenBLASLUFactorization <: AbstractFactorization end -module OpenBLASLU +# OpenBLAS methods - OpenBLAS_jll is always available as a standard library -using LinearAlgebra -using LinearAlgebra.LAPACK: chkfinite, chkstride1, @blasfunc, chkargsok, chktrans, - chklapackerror -using OpenBLAS_jll - -function getrf!(A::AbstractMatrix{<:ComplexF64}; - ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))), - info = Ref{LinearAlgebra.BlasInt}(), +function openblas_getrf!(A::AbstractMatrix{<:ComplexF64}; + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{BlasInt}(), check = false) - LinearAlgebra.require_one_based_indexing(A) + require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) m, n = size(A) lda = max(1, stride(A, 2)) if isempty(ipiv) - ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))) + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))) end ccall((@blasfunc(zgetrf_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ptr{ComplexF64}, - Ref{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}), + (Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF64}, + Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), m, n, A, lda, ipiv, info) chkargsok(info[]) A, ipiv, info[], info #Error code is stored in LU factorization type end -function getrf!(A::AbstractMatrix{<:ComplexF32}; - ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))), - info = Ref{LinearAlgebra.BlasInt}(), +function openblas_getrf!(A::AbstractMatrix{<:ComplexF32}; + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{BlasInt}(), check = false) - LinearAlgebra.require_one_based_indexing(A) + require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) m, n = size(A) lda = max(1, stride(A, 2)) if isempty(ipiv) - ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))) + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))) end ccall((@blasfunc(cgetrf_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ptr{ComplexF32}, - Ref{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}), + (Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF32}, + Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), m, n, A, lda, ipiv, info) chkargsok(info[]) A, ipiv, info[], info #Error code is stored in LU factorization type end -function getrf!(A::AbstractMatrix{<:Float64}; - ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))), - info = Ref{LinearAlgebra.BlasInt}(), +function openblas_getrf!(A::AbstractMatrix{<:Float64}; + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{BlasInt}(), check = false) - LinearAlgebra.require_one_based_indexing(A) + require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) m, n = size(A) lda = max(1, stride(A, 2)) if isempty(ipiv) - ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))) + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))) end ccall((@blasfunc(dgetrf_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ptr{Float64}, - Ref{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}), + (Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, + Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), m, n, A, lda, ipiv, info) chkargsok(info[]) A, ipiv, info[], info #Error code is stored in LU factorization type end -function getrf!(A::AbstractMatrix{<:Float32}; - ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))), - info = Ref{LinearAlgebra.BlasInt}(), +function openblas_getrf!(A::AbstractMatrix{<:Float32}; + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{BlasInt}(), check = false) - LinearAlgebra.require_one_based_indexing(A) + require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) m, n = size(A) lda = max(1, stride(A, 2)) if isempty(ipiv) - ipiv = similar(A, LinearAlgebra.BlasInt, min(size(A, 1), size(A, 2))) + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))) end ccall((@blasfunc(sgetrf_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ptr{Float32}, - Ref{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}), + (Ref{BlasInt}, Ref{BlasInt}, Ptr{Float32}, + Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), m, n, A, lda, ipiv, info) chkargsok(info[]) A, ipiv, info[], info #Error code is stored in LU factorization type end -function getrs!(trans::AbstractChar, +function openblas_getrs!(trans::AbstractChar, A::AbstractMatrix{<:ComplexF64}, - ipiv::AbstractVector{LinearAlgebra.BlasInt}, + ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{<:ComplexF64}; - info = Ref{LinearAlgebra.BlasInt}()) - LinearAlgebra.require_one_based_indexing(A, ipiv, B) + info = Ref{BlasInt}()) + require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) n = LinearAlgebra.checksquare(A) @@ -137,22 +132,20 @@ function getrs!(trans::AbstractChar, end nrhs = size(B, 2) ccall((@blasfunc(zgetrs_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{UInt8}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, - Ptr{ComplexF64}, Ref{LinearAlgebra.BlasInt}, - Ptr{LinearAlgebra.BlasInt}, Ptr{ComplexF64}, Ref{LinearAlgebra.BlasInt}, - Ptr{LinearAlgebra.BlasInt}, Clong), + (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF64}, Ref{BlasInt}, + Ptr{BlasInt}, Ptr{ComplexF64}, Ref{BlasInt}, Ptr{BlasInt}, Clong), trans, n, size(B, 2), A, max(1, stride(A, 2)), ipiv, B, max(1, stride(B, 2)), info, 1) - LinearAlgebra.LAPACK.chklapackerror(LinearAlgebra.BlasInt(info[])) + LinearAlgebra.LAPACK.chklapackerror(BlasInt(info[])) B end -function getrs!(trans::AbstractChar, +function openblas_getrs!(trans::AbstractChar, A::AbstractMatrix{<:ComplexF32}, - ipiv::AbstractVector{LinearAlgebra.BlasInt}, + ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{<:ComplexF32}; - info = Ref{LinearAlgebra.BlasInt}()) - LinearAlgebra.require_one_based_indexing(A, ipiv, B) + info = Ref{BlasInt}()) + require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) n = LinearAlgebra.checksquare(A) @@ -164,22 +157,20 @@ function getrs!(trans::AbstractChar, end nrhs = size(B, 2) ccall((@blasfunc(cgetrs_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{UInt8}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, - Ptr{ComplexF32}, Ref{LinearAlgebra.BlasInt}, - Ptr{LinearAlgebra.BlasInt}, Ptr{ComplexF32}, Ref{LinearAlgebra.BlasInt}, - Ptr{LinearAlgebra.BlasInt}, Clong), + (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF32}, Ref{BlasInt}, + Ptr{BlasInt}, Ptr{ComplexF32}, Ref{BlasInt}, Ptr{BlasInt}, Clong), trans, n, size(B, 2), A, max(1, stride(A, 2)), ipiv, B, max(1, stride(B, 2)), info, 1) - LinearAlgebra.LAPACK.chklapackerror(LinearAlgebra.BlasInt(info[])) + LinearAlgebra.LAPACK.chklapackerror(BlasInt(info[])) B end -function getrs!(trans::AbstractChar, +function openblas_getrs!(trans::AbstractChar, A::AbstractMatrix{<:Float64}, - ipiv::AbstractVector{LinearAlgebra.BlasInt}, + ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{<:Float64}; - info = Ref{LinearAlgebra.BlasInt}()) - LinearAlgebra.require_one_based_indexing(A, ipiv, B) + info = Ref{BlasInt}()) + require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) n = LinearAlgebra.checksquare(A) @@ -191,21 +182,20 @@ function getrs!(trans::AbstractChar, end nrhs = size(B, 2) ccall((@blasfunc(dgetrs_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{UInt8}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, - Ptr{Float64}, Ref{LinearAlgebra.BlasInt}, - Ptr{LinearAlgebra.BlasInt}, Ptr{Float64}, Ref{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}, Clong), + (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, + Ptr{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{BlasInt}, Clong), trans, n, size(B, 2), A, max(1, stride(A, 2)), ipiv, B, max(1, stride(B, 2)), info, 1) - LinearAlgebra.LAPACK.chklapackerror(LinearAlgebra.BlasInt(info[])) + LinearAlgebra.LAPACK.chklapackerror(BlasInt(info[])) B end -function getrs!(trans::AbstractChar, +function openblas_getrs!(trans::AbstractChar, A::AbstractMatrix{<:Float32}, - ipiv::AbstractVector{LinearAlgebra.BlasInt}, + ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{<:Float32}; - info = Ref{LinearAlgebra.BlasInt}()) - LinearAlgebra.require_one_based_indexing(A, ipiv, B) + info = Ref{BlasInt}()) + require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) n = LinearAlgebra.checksquare(A) @@ -217,23 +207,20 @@ function getrs!(trans::AbstractChar, end nrhs = size(B, 2) ccall((@blasfunc(sgetrs_), OpenBLAS_jll.libopenblas), Cvoid, - (Ref{UInt8}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, - Ptr{Float32}, Ref{LinearAlgebra.BlasInt}, - Ptr{LinearAlgebra.BlasInt}, Ptr{Float32}, Ref{LinearAlgebra.BlasInt}, Ptr{LinearAlgebra.BlasInt}, Clong), + (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{Float32}, Ref{BlasInt}, + Ptr{BlasInt}, Ptr{Float32}, Ref{BlasInt}, Ptr{BlasInt}, Clong), trans, n, size(B, 2), A, max(1, stride(A, 2)), ipiv, B, max(1, stride(B, 2)), info, 1) - LinearAlgebra.LAPACK.chklapackerror(LinearAlgebra.BlasInt(info[])) + LinearAlgebra.LAPACK.chklapackerror(BlasInt(info[])) B end -end # module OpenBLASLU - default_alias_A(::OpenBLASLUFactorization, ::Any, ::Any) = false default_alias_b(::OpenBLASLUFactorization, ::Any, ::Any) = false const PREALLOCATED_OPENBLAS_LU = begin A = rand(0, 0) - luinst = ArrayInterface.lu_instance(A), Ref{LinearAlgebra.BlasInt}() + luinst = ArrayInterface.lu_instance(A), Ref{BlasInt}() end function LinearSolve.init_cacheval(alg::OpenBLASLUFactorization, A, b, u, Pl, Pr, @@ -247,7 +234,7 @@ function LinearSolve.init_cacheval(alg::OpenBLASLUFactorization, maxiters::Int, abstol, reltol, verbose::LinearVerbosity, assumptions::OperatorAssumptions) A = rand(eltype(A), 0, 0) - ArrayInterface.lu_instance(A), Ref{LinearAlgebra.BlasInt}() + ArrayInterface.lu_instance(A), Ref{BlasInt}() end function SciMLBase.solve!(cache::LinearCache, alg::OpenBLASLUFactorization; @@ -256,8 +243,8 @@ function SciMLBase.solve!(cache::LinearCache, alg::OpenBLASLUFactorization; A = convert(AbstractMatrix, A) if cache.isfresh cacheval = @get_cacheval(cache, :OpenBLASLUFactorization) - res = OpenBLASLU.getrf!(A; ipiv = cacheval[1].ipiv, info = cacheval[2]) - fact = LinearAlgebra.LU(res[1:3]...), res[4] + res = openblas_getrf!(A; ipiv = cacheval[1].ipiv, info = cacheval[2]) + fact = LU(res[1:3]...), res[4] cache.cacheval = fact if !LinearAlgebra.issuccess(fact[1]) @@ -268,15 +255,15 @@ function SciMLBase.solve!(cache::LinearCache, alg::OpenBLASLUFactorization; end A, info = @get_cacheval(cache, :OpenBLASLUFactorization) - LinearAlgebra.require_one_based_indexing(cache.u, cache.b) + require_one_based_indexing(cache.u, cache.b) m, n = size(A, 1), size(A, 2) if m > n Bc = copy(cache.b) - OpenBLASLU.getrs!('N', A.factors, A.ipiv, Bc; info) + openblas_getrs!('N', A.factors, A.ipiv, Bc; info) copyto!(cache.u, 1, Bc, 1, n) else copyto!(cache.u, cache.b) - OpenBLASLU.getrs!('N', A.factors, A.ipiv, cache.u; info) + openblas_getrs!('N', A.factors, A.ipiv, cache.u; info) end SciMLBase.build_linear_solution( diff --git a/test/basictests.jl b/test/basictests.jl index f360d87a1..f7f115245 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -287,8 +287,10 @@ end push!(test_algs, MKLLUFactorization()) end - # Always test OpenBLAS since it's a direct dependency - push!(test_algs, OpenBLASLUFactorization()) + # Test OpenBLAS if available + if LinearSolve.useopenblas + push!(test_algs, OpenBLASLUFactorization()) + end # Test BLIS if extension is available if Base.get_extension(LinearSolve, :LinearSolveBLISExt) !== nothing diff --git a/test/preferences.jl b/test/preferences.jl index 0182657b4..d3118610e 100644 --- a/test/preferences.jl +++ b/test/preferences.jl @@ -190,11 +190,13 @@ using Preferences println("✅ MKLLUFactorization confirmed working") end - # Test OpenBLAS (always available as a dependency) - sol_openblas = solve(prob, OpenBLASLUFactorization()) - @test sol_openblas.retcode == ReturnCode.Success - @test norm(A * sol_openblas.u - b) < 1e-8 - println("✅ OpenBLASLUFactorization confirmed working") + # Test OpenBLAS if available + if LinearSolve.useopenblas + sol_openblas = solve(prob, OpenBLASLUFactorization()) + @test sol_openblas.retcode == ReturnCode.Success + @test norm(A * sol_openblas.u - b) < 1e-8 + println("✅ OpenBLASLUFactorization confirmed working") + end # Test Apple Accelerate if available if LinearSolve.appleaccelerate_isavailable() diff --git a/test/resolve.jl b/test/resolve.jl index 0f8bdbea0..932df173a 100644 --- a/test/resolve.jl +++ b/test/resolve.jl @@ -20,7 +20,7 @@ for alg in vcat(InteractiveUtils.subtypes(AbstractDenseFactorization), (!(alg == AppleAccelerateLUFactorization) || LinearSolve.appleaccelerate_isavailable()) && (!(alg == MKLLUFactorization) || LinearSolve.usemkl) && - (!(alg == OpenBLASLUFactorization) || true) # OpenBLAS is always available as a dependency + (!(alg == OpenBLASLUFactorization) || LinearSolve.useopenblas) A = [1.0 2.0; 3.0 4.0] alg in [KLUFactorization, UMFPACKFactorization, SparspakFactorization] && (A = sparse(A)) From 1bb9e5a94fbe810a4325705d9f5e3fad9919851f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 20 Aug 2025 17:21:32 -0400 Subject: [PATCH 6/7] Update src/LinearSolve.jl --- src/LinearSolve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index 10ae83894..4be5cb595 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -60,7 +60,7 @@ else end # OpenBLAS_jll is a standard library, always available -using OpenBLAS_jll +using OpenBLAS_jll: OpenBLAS_jll const useopenblas = OpenBLAS_jll.is_available() From aecd1afa010702f5f8da074a86f40d60ab168843 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 20 Aug 2025 17:27:49 -0400 Subject: [PATCH 7/7] Add preference for loading OpenBLAS, matching MKL pattern - Added LoadOpenBLAS_JLL preference with default value true - Users can now disable OpenBLAS via preferences - Matches the existing MKL preference pattern --- src/LinearSolve.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index 4be5cb595..fa20b7bd1 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -59,9 +59,13 @@ else const usemkl = false end -# OpenBLAS_jll is a standard library, always available -using OpenBLAS_jll: OpenBLAS_jll -const useopenblas = OpenBLAS_jll.is_available() +# OpenBLAS_jll is a standard library, but allow users to disable it via preferences +if Preferences.@load_preference("LoadOpenBLAS_JLL", true) + using OpenBLAS_jll: OpenBLAS_jll + const useopenblas = OpenBLAS_jll.is_available() +else + const useopenblas = false +end @reexport using SciMLBase