Merge pull request #506 from j-fu/pardiso-vendor

Allow to choose Pardiso vendor

Author: ChrisRackauckas

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 @@ function LinearSolve.init_cacheval(alg::PardisoJL,
         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 @@ end
 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 @@ 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.
 """
-MKLPardisoFactorize(; kwargs...) = PardisoJL(; solver_type = 0, kwargs...)
+MKLPardisoFactorize(; kwargs...) = PardisoJL(; vendor = :MKL, solver_type = 0, kwargs...)
 
 """
 ```julia
@@ -136,19 +136,18 @@ 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.
 """
-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 @@ 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.
 """
+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 @@ struct PardisoJL{T1, T2} <: LinearSolve.SciMLLinearSolveAlgorithm
             @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