Allow to choose Pardiso vendor

Project.toml

@@ -96,7 +96,7 @@ MPI = "0.20"
 Markdown = "1.10"
 Metal = "1"
 MultiFloats = "1"
-Pardiso = "0.5"
+Pardiso = "0.5.7"
 Pkg = "1"
 PrecompileTools = "1.2"
 Preferences = "1.4"

ext/LinearSolvePardisoExt.jl

@@ -22,26 +22,44 @@
         reltol,
         verbose::Bool,
         assumptions::LinearSolve.OperatorAssumptions)
-    @unpack nprocs, solver_type, matrix_type, cache_analysis, iparm, dparm = alg
+    @unpack nprocs, solver_type, matrix_type, cache_analysis, iparm, dparm, vendor = alg
     A = convert(AbstractMatrix, A)
 
+    if isnothing(vendor)
+        if Pardiso.panua_is_available()
+            vendor = :Panua
+        else
+            vendor = :MKL
+        end
+    end
+
     transposed_iparm = 1
-    solver = if Pardiso.PARDISO_LOADED[]
-        solver = Pardiso.PardisoSolver()
-        Pardiso.pardisoinit(solver)
-        solver_type !== nothing && Pardiso.set_solver!(solver, solver_type)
+    solver = if vendor == :MKL
+        solver = if Pardiso.mkl_is_available()
+            solver = Pardiso.MKLPardisoSolver()
+            Pardiso.pardisoinit(solver)
+            nprocs !== nothing && Pardiso.set_nprocs!(solver, nprocs)
 
-        solver
-    else
-        solver = Pardiso.MKLPardisoSolver()
-        Pardiso.pardisoinit(solver)
-        nprocs !== nothing && Pardiso.set_nprocs!(solver, nprocs)
+            # for mkl 1 means conjugated and 2 means transposed.
+            # https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-0/pardiso-iparm-parameter.html#IPARM37
+            transposed_iparm = 2
 
-        # for mkl 1 means conjugated and 2 means transposed.
-        # https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-0/pardiso-iparm-parameter.html#IPARM37
-        transposed_iparm = 2
+            solver
+        else
+            error("MKL Pardiso is not available. On MacOSX, possibly, try Panua Pardiso.")
+        end
+    elseif vendor == :Panua
+        solver = if Pardiso.panua_is_available()
+            solver = Pardiso.PardisoSolver()
+            Pardiso.pardisoinit(solver)
+            solver_type !== nothing && Pardiso.set_solver!(solver, solver_type)
 
-        solver
+            solver
+        else
+            error("Panua Pardiso is not available.")
+        end
+    else
+        error("Pardiso vendor must be either `:MKL` or `:Panua`")
     end
 
     if matrix_type !== nothing
@@ -113,7 +131,6 @@
 function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::PardisoJL; kwargs...)
     @unpack A, b, u = cache
     A = convert(AbstractMatrix, A)
-
     if cache.isfresh
         phase = alg.cache_analysis ? Pardiso.NUM_FACT : Pardiso.ANALYSIS_NUM_FACT
         Pardiso.set_phase!(cache.cacheval, phase)

src/LinearSolve.jl

@@ -253,6 +253,7 @@ export SimpleGMRES
 export HYPREAlgorithm
 export CudaOffloadFactorization
 export MKLPardisoFactorize, MKLPardisoIterate
+export PanuaPardisoFactorize, PanuaPardisoIterate
 export PardisoJL
 export MKLLUFactorization
 export AppleAccelerateLUFactorization

src/extension_algs.jl

@@ -108,7 +108,7 @@
 given the input types, and these keyword arguments are only for overriding the
 default handling process. This should not be required by most users.
 """
-MKLPardisoFactorize(; kwargs...) = PardisoJL(; solver_type = 0, kwargs...)
+MKLPardisoFactorize(; kwargs...) = PardisoJL(; vendor = :MKL, solver_type = 0, kwargs...)
 
 """
 ```julia
@@ -136,19 +136,18 @@
 given the input types, and these keyword arguments are only for overriding the
 default handling process. This should not be required by most users.
 """
-MKLPardisoIterate(; kwargs...) = PardisoJL(; solver_type = 1, kwargs...)
+MKLPardisoIterate(; kwargs...) = PardisoJL(; vendor = :MKL, solver_type = 1, kwargs...)
 
 """
 ```julia
-PardisoJL(; nprocs::Union{Int, Nothing} = nothing,
-    solver_type = nothing,
+PanuaPardisoFactorize(; nprocs::Union{Int, Nothing} = nothing,
     matrix_type = nothing,
     cache_analysis = false,
     iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
     dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing)
 ```
 
-A generic method using MKL Pardiso. Specifying `solver_type` is required.
+A sparse factorization method using Panua Pardiso.
 
 !!! note
 
@@ -165,20 +164,75 @@
 given the input types, and these keyword arguments are only for overriding the
 default handling process. This should not be required by most users.
 """
+PanuaPardisoFactorize(; kwargs...) = PardisoJL(;
+    vendor = :Panua, solver_type = 0, kwargs...)
+
+"""
+```julia
+PanuaPardisoIterate(; nprocs::Union{Int, Nothing} = nothing,
+    matrix_type = nothing,
+    iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
+    dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing)
+```
+
+A mixed factorization+iterative method using Panua Pardiso.
+
+!!! note
+
+    Using this solver requires adding the package Pardiso.jl, i.e. `using Pardiso`
+
+## Keyword Arguments
+
+For the definition of the keyword arguments, see the Pardiso.jl documentation.
+All values default to `nothing` and the solver internally determines the values
+given the input types, and these keyword arguments are only for overriding the
+default handling process. This should not be required by most users.
+"""
+PanuaPardisoIterate(; kwargs...) = PardisoJL(; vendor = :Panua, solver_type = 1, kwargs...)
+
+"""
+```julia
+PardisoJL(; nprocs::Union{Int, Nothing} = nothing,
+    solver_type = nothing,
+    matrix_type = nothing,
+    iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
+    dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
+    vendor::Union{Symbol, Nothing} = nothing
+)
+```
+
+A generic method using  Pardiso. Specifying `solver_type` is required.
+
+!!! note
+
+    Using this solver requires adding the package Pardiso.jl, i.e. `using Pardiso`
+
+## Keyword Arguments
+
+The `vendor` keyword allows to choose between Panua pardiso  (former pardiso-project.org; `vendor=:Panua`)
+and  MKL Pardiso (`vendor=:MKL`). If `vendor==nothing`, Panua pardiso is preferred over MKL Pardiso.
+
+For the definition of the other keyword arguments, see the Pardiso.jl documentation.
+All values default to `nothing` and the solver internally determines the values
+given the input types, and these keyword arguments are only for overriding the
+default handling process. This should not be required by most users.
+"""
 struct PardisoJL{T1, T2} <: LinearSolve.SciMLLinearSolveAlgorithm
     nprocs::Union{Int, Nothing}
     solver_type::T1
     matrix_type::T2
     cache_analysis::Bool
     iparm::Union{Vector{Tuple{Int, Int}}, Nothing}
     dparm::Union{Vector{Tuple{Int, Int}}, Nothing}
+    vendor::Union{Symbol, Nothing}
 
     function PardisoJL(; nprocs::Union{Int, Nothing} = nothing,
             solver_type = nothing,
             matrix_type = nothing,
             cache_analysis = false,
             iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
-            dparm::Union{Vector{Tuple{Int, Float64}}, Nothing} = nothing)
+            dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
+            vendor::Union{Symbol, Nothing} = nothing)
         ext = Base.get_extension(@__MODULE__, :LinearSolvePardisoExt)
         if ext === nothing
             error("PardisoJL requires that Pardiso is loaded, i.e. `using Pardiso`")
@@ -188,7 +242,7 @@
             @assert T1 <: Union{Int, Nothing, ext.Pardiso.Solver}
             @assert T2 <: Union{Int, Nothing, ext.Pardiso.MatrixType}
             return new{T1, T2}(
-                nprocs, solver_type, matrix_type, cache_analysis, iparm, dparm)
+                nprocs, solver_type, matrix_type, cache_analysis, iparm, dparm, vendor)
         end
     end
 end