From e864de3d25c1dbb0198c25fce4f75b301672f8a5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Mon, 3 Jun 2024 10:50:54 +0200 Subject: [PATCH 1/3] Allow code based control of usage of (Panua) Pardiso and MKL Pardiso. --- Project.toml | 2 +- ext/LinearSolvePardisoExt.jl | 38 ++++++++++++++------ src/LinearSolve.jl | 1 + src/extension_algs.jl | 69 ++++++++++++++++++++++++++++++++---- test/pardiso/pardiso.jl | 17 ++++++--- 5 files changed, 105 insertions(+), 22 deletions(-) diff --git a/Project.toml b/Project.toml index 834a6b246..5056f3ff1 100644 --- a/Project.toml +++ b/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" diff --git a/ext/LinearSolvePardisoExt.jl b/ext/LinearSolvePardisoExt.jl index 60312ce66..c454bd2f5 100644 --- a/ext/LinearSolvePardisoExt.jl +++ b/ext/LinearSolvePardisoExt.jl @@ -22,21 +22,39 @@ function LinearSolve.init_cacheval(alg::PardisoJL, reltol, verbose::Bool, assumptions::LinearSolve.OperatorAssumptions) - @unpack nprocs, solver_type, matrix_type, iparm, dparm = alg + @unpack nprocs, solver_type, matrix_type, iparm, dparm, vendor = alg A = convert(AbstractMatrix, A) - solver = if Pardiso.PARDISO_LOADED[] - solver = Pardiso.PardisoSolver() - solver_type !== nothing && Pardiso.set_solver!(solver, solver_type) + if isnothing(vendor) + if Pardiso.panua_is_available() + vendor=:Panua + else + vendor=:MKL + end + end + + solver = if vendor == :MKL + solver = if Pardiso.mkl_is_available() + solver = Pardiso.MKLPardisoSolver() + nprocs !== nothing && Pardiso.set_nprocs!(solver, nprocs) - solver + 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() + solver_type !== nothing && Pardiso.set_solver!(solver, solver_type) + + solver + else + error("Panua Pardiso is not available.") + end else - solver = Pardiso.MKLPardisoSolver() - nprocs !== nothing && Pardiso.set_nprocs!(solver, nprocs) - - solver + error("Pardiso vendor must be either `:MKL` or `:Panua`") end - + Pardiso.pardisoinit(solver) # default initialization if matrix_type !== nothing diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index 5df86c925..fcec925d8 100644 --- a/src/LinearSolve.jl +++ b/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 diff --git a/src/extension_algs.jl b/src/extension_algs.jl index 4c2d4fecb..b2349b74f 100644 --- a/src/extension_algs.jl +++ b/src/extension_algs.jl @@ -103,7 +103,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 @@ -126,18 +126,41 @@ 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, + iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing, + dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing) +``` + +A sparse factorization 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. +""" +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 generic method using MKL Pardiso. Specifying `solver_type` is required. +A mixed factorization+iterative method using Panua Pardiso. !!! note @@ -150,18 +173,50 @@ 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 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, iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing, - dparm::Union{Vector{Tuple{Int, Int}}, 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`") @@ -170,7 +225,7 @@ struct PardisoJL{T1, T2} <: LinearSolve.SciMLLinearSolveAlgorithm T2 = typeof(matrix_type) @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, iparm, dparm) + return new{T1, T2}(nprocs, solver_type, matrix_type, iparm, dparm, vendor) end end end diff --git a/test/pardiso/pardiso.jl b/test/pardiso/pardiso.jl index 6753be21b..58a04af46 100644 --- a/test/pardiso/pardiso.jl +++ b/test/pardiso/pardiso.jl @@ -14,17 +14,26 @@ e = ones(n) e2 = ones(n - 1) A2 = spdiagm(-1 => im * e2, 0 => lambda * e, 1 => -im * e2) b2 = rand(n) + im * zeros(n) -cache_kwargs = (; verbose = true, abstol = 1e-8, reltol = 1e-8, maxiter = 30) +cache_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiter = 30) prob2 = LinearProblem(A2, b2) -for alg in (PardisoJL(), MKLPardisoFactorize(), MKLPardisoIterate()) +algs=[] +# if Pardiso.mkl_is_available() +# algs=vcat(algs,[PardisoJL(), MKLPardisoFactorize(), MKLPardisoIterate()]) +# end + +if Pardiso.panua_is_available() + algs=vcat(algs,[PardisoJL(), PanuaPardisoFactorize(), PanuaPardisoIterate()]) +end +@info algs +for alg in algs u = solve(prob1, alg; cache_kwargs...).u @test A1 * u ≈ b1 u = solve(prob2, alg; cache_kwargs...).u - @test eltype(u) <: Complex - @test_broken A2 * u ≈ b2 +# @test eltype(u) <: Complex +# @test A2 * u ≈ b2 end Random.seed!(10) From 835f224b7c00c3928a4c349a5496d5c6741fd601 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Mon, 3 Jun 2024 11:08:31 +0200 Subject: [PATCH 2/3] Fix iparm value for transposed matrix for MKL Pardiso --- ext/LinearSolvePardisoExt.jl | 14 ++++++++------ test/pardiso/pardiso.jl | 18 ++++++++++-------- 2 files changed, 18 insertions(+), 14 deletions(-) diff --git a/ext/LinearSolvePardisoExt.jl b/ext/LinearSolvePardisoExt.jl index c454bd2f5..5bf3094cc 100644 --- a/ext/LinearSolvePardisoExt.jl +++ b/ext/LinearSolvePardisoExt.jl @@ -32,11 +32,13 @@ function LinearSolve.init_cacheval(alg::PardisoJL, vendor=:MKL end end - + + transposed_iparm = 1 solver = if vendor == :MKL solver = if Pardiso.mkl_is_available() solver = Pardiso.MKLPardisoSolver() nprocs !== nothing && Pardiso.set_nprocs!(solver, nprocs) + transposed_iparm = 2 # see https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-0/pardiso-iparm-parameter.html#IPARM37 solver else @@ -84,7 +86,7 @@ function LinearSolve.init_cacheval(alg::PardisoJL, end # Make sure to say it's transposed because its CSC not CSR - Pardiso.set_iparm!(solver, 12, 1) + Pardiso.set_iparm!(solver, 12, transposed_iparm) #= Note: It is recommended to use IPARM(11)=1 (scaling) and IPARM(13)=1 (matchings) for @@ -107,12 +109,12 @@ function LinearSolve.init_cacheval(alg::PardisoJL, # applies these exact factors L and U for the next steps in a # preconditioned Krylov-Subspace iteration. If the iteration does not # converge, the solver will automatically switch back to the numerical factorization. - Pardiso.set_iparm!(solver, 3, round(Int, abs(log10(reltol)), RoundDown) * 10 + 1) + Pardiso.set_iparm!(solver, 4, round(Int, abs(log10(reltol)), RoundDown) * 10 + 1) end Pardiso.pardiso(solver, - u, - SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)), + u, + SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)), b) return solver @@ -121,7 +123,6 @@ end function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::PardisoJL; kwargs...) @unpack A, b, u = cache A = convert(AbstractMatrix, A) - if cache.isfresh Pardiso.set_phase!(cache.cacheval, Pardiso.NUM_FACT) Pardiso.pardiso(cache.cacheval, A, eltype(A)[]) @@ -130,6 +131,7 @@ function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::PardisoJL; kwargs Pardiso.set_phase!(cache.cacheval, Pardiso.SOLVE_ITERATIVE_REFINE) Pardiso.pardiso(cache.cacheval, u, A, b) + return SciMLBase.build_linear_solution(alg, cache.u, nothing, cache) end diff --git a/test/pardiso/pardiso.jl b/test/pardiso/pardiso.jl index 58a04af46..957238180 100644 --- a/test/pardiso/pardiso.jl +++ b/test/pardiso/pardiso.jl @@ -18,22 +18,24 @@ cache_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiter = 30) prob2 = LinearProblem(A2, b2) -algs=[] -# if Pardiso.mkl_is_available() -# algs=vcat(algs,[PardisoJL(), MKLPardisoFactorize(), MKLPardisoIterate()]) -# end +algs=[PardisoJL()] + +if Pardiso.mkl_is_available() + algs=vcat(algs,[MKLPardisoFactorize(), MKLPardisoIterate()]) +end if Pardiso.panua_is_available() - algs=vcat(algs,[PardisoJL(), PanuaPardisoFactorize(), PanuaPardisoIterate()]) + algs=vcat(algs,[PanuaPardisoFactorize(), PanuaPardisoIterate()]) end -@info algs + + for alg in algs u = solve(prob1, alg; cache_kwargs...).u @test A1 * u ≈ b1 u = solve(prob2, alg; cache_kwargs...).u -# @test eltype(u) <: Complex -# @test A2 * u ≈ b2 + @test eltype(u) <: Complex + @test A2 * u ≈ b2 end Random.seed!(10) From 0c0de1ca88fedf36e58454f2efd2d07ad82aed7c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Wed, 5 Jun 2024 12:33:47 +0200 Subject: [PATCH 3/3] pardiso: format --- ext/LinearSolvePardisoExt.jl | 11 +- src/extension_algs.jl | 17 ++- test/pardiso/pardiso.jl | 203 +++++++++++++++++++---------------- 3 files changed, 122 insertions(+), 109 deletions(-) diff --git a/ext/LinearSolvePardisoExt.jl b/ext/LinearSolvePardisoExt.jl index 43ca1e2f6..0b4cfbb11 100644 --- a/ext/LinearSolvePardisoExt.jl +++ b/ext/LinearSolvePardisoExt.jl @@ -27,9 +27,9 @@ function LinearSolve.init_cacheval(alg::PardisoJL, if isnothing(vendor) if Pardiso.panua_is_available() - vendor=:Panua + vendor = :Panua else - vendor=:MKL + vendor = :MKL end end @@ -42,7 +42,7 @@ function LinearSolve.init_cacheval(alg::PardisoJL, # 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 + transposed_iparm = 2 solver else @@ -53,7 +53,7 @@ function LinearSolve.init_cacheval(alg::PardisoJL, solver = Pardiso.PardisoSolver() Pardiso.pardisoinit(solver) solver_type !== nothing && Pardiso.set_solver!(solver, solver_type) - + solver else error("Panua Pardiso is not available.") @@ -61,7 +61,7 @@ function LinearSolve.init_cacheval(alg::PardisoJL, else error("Pardiso vendor must be either `:MKL` or `:Panua`") end - + if matrix_type !== nothing Pardiso.set_matrixtype!(solver, matrix_type) else @@ -125,7 +125,6 @@ function LinearSolve.init_cacheval(alg::PardisoJL, b) end ->>>>>>> main return solver end diff --git a/src/extension_algs.jl b/src/extension_algs.jl index f4d5f4012..c5e2db955 100644 --- a/src/extension_algs.jl +++ b/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(; vendor=:MKL, solver_type = 0, kwargs...) +MKLPardisoFactorize(; kwargs...) = PardisoJL(; vendor = :MKL, solver_type = 0, kwargs...) """ ```julia @@ -136,8 +136,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. """ -MKLPardisoIterate(; kwargs...) = PardisoJL(; vendor=:MKL, solver_type = 1, kwargs...) - +MKLPardisoIterate(; kwargs...) = PardisoJL(; vendor = :MKL, solver_type = 1, kwargs...) """ ```julia @@ -165,7 +164,8 @@ 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...) +PanuaPardisoFactorize(; kwargs...) = PardisoJL(; + vendor = :Panua, solver_type = 0, kwargs...) """ ```julia @@ -188,8 +188,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. """ -PanuaPardisoIterate(; kwargs...) = PardisoJL(; vendor=:Panua, solver_type = 1, kwargs...) - +PanuaPardisoIterate(; kwargs...) = PardisoJL(; vendor = :Panua, solver_type = 1, kwargs...) """ ```julia @@ -198,7 +197,7 @@ PardisoJL(; nprocs::Union{Int, Nothing} = 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 + vendor::Union{Symbol, Nothing} = nothing ) ``` @@ -225,7 +224,7 @@ struct PardisoJL{T1, T2} <: LinearSolve.SciMLLinearSolveAlgorithm cache_analysis::Bool iparm::Union{Vector{Tuple{Int, Int}}, Nothing} dparm::Union{Vector{Tuple{Int, Int}}, Nothing} - vendor::Union{Symbol,Nothing} + vendor::Union{Symbol, Nothing} function PardisoJL(; nprocs::Union{Int, Nothing} = nothing, solver_type = nothing, @@ -233,7 +232,7 @@ struct PardisoJL{T1, T2} <: LinearSolve.SciMLLinearSolveAlgorithm cache_analysis = false, iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing, dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing, - vendor::Union{Symbol,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`") diff --git a/test/pardiso/pardiso.jl b/test/pardiso/pardiso.jl index 9ffa6ec00..d244c4e33 100644 --- a/test/pardiso/pardiso.jl +++ b/test/pardiso/pardiso.jl @@ -18,18 +18,26 @@ cache_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiter = 30) prob2 = LinearProblem(A2, b2) -algs=[PardisoJL()] +algs=LinearSolve.SciMLLinearSolveAlgorithm[PardisoJL()] +solvers=Pardiso.AbstractPardisoSolver[] +extended_algs=LinearSolve.SciMLLinearSolveAlgorithm[PardisoJL()] if Pardiso.mkl_is_available() - algs=vcat(algs,[MKLPardisoFactorize(), MKLPardisoIterate()]) + push!(algs,MKLPardisoFactorize()) + push!(solvers,Pardiso.MKLPardisoSolver()) + extended_algs=vcat(extended_algs,[MKLPardisoFactorize(), MKLPardisoIterate()]) + @info "Testing MKL Pardiso" end if Pardiso.panua_is_available() - algs=vcat(algs,[PanuaPardisoFactorize(), PanuaPardisoIterate()]) + push!(algs,PanuaPardisoFactorize()) + push!(solvers,Pardiso.PardisoSolver()) + extended_algs=vcat(extended_algs,[PanuaPardisoFactorize(), PanuaPardisoIterate()]) + @info "Testing Panua Pardiso" end -for alg in algs +for alg in extended_algs u = solve(prob1, alg; cache_kwargs...).u @test A1 * u ≈ b1 @@ -38,7 +46,6 @@ for alg in algs @test A2 * u ≈ b2 end -return Random.seed!(10) @@ -48,7 +55,7 @@ b1 = rand(n); b2 = rand(n); prob = LinearProblem(copy(A), copy(b1)) -prob = LinearProblem(copy(A), copy(b1)) + linsolve = init(prob, UMFPACKFactorization()) sol11 = solve(linsolve) linsolve = LinearSolve.set_b(sol11.cache, copy(b2)) @@ -56,16 +63,19 @@ sol12 = solve(linsolve) linsolve = LinearSolve.set_A(sol12.cache, copy(A2)) sol13 = solve(linsolve) -linsolve = init(prob, MKLPardisoFactorize()) -sol31 = solve(linsolve) -linsolve = LinearSolve.set_b(sol31.cache, copy(b2)) -sol32 = solve(linsolve) -linsolve = LinearSolve.set_A(sol32.cache, copy(A2)) -sol33 = solve(linsolve) -@test sol11.u ≈ sol31.u -@test sol12.u ≈ sol32.u -@test sol13.u ≈ sol33.u + +for alg in algs + linsolve = init(prob, alg) + sol31 = solve(linsolve) + linsolve = LinearSolve.set_b(sol31.cache, copy(b2)) + sol32 = solve(linsolve) + linsolve = LinearSolve.set_A(sol32.cache, copy(A2)) + sol33 = solve(linsolve) + @test sol11.u ≈ sol31.u + @test sol12.u ≈ sol32.u + @test sol13.u ≈ sol33.u +end # Test for problem from #497 @@ -78,87 +88,92 @@ function makeA() return(A) end -A=makeA() -u0=fill(0.1,size(A,2)) -linprob = LinearProblem(A, A*u0) -u = LinearSolve.solve(linprob, PardisoJL()) -@test norm(u-u0) < 1.0e-14 +for alg in algs + A=makeA() + u0=fill(0.1,size(A,2)) + linprob = LinearProblem(A, A*u0) + u = LinearSolve.solve(linprob, alg) + @test norm(u-u0) < 1.0e-14 +end -# Testing and demonstrating Pardiso.set_iparm! for MKLPardisoSolver -solver = Pardiso.MKLPardisoSolver() -iparm = [ - (1, 1), - (2, 2), - (3, 0), - (4, 0), - (5, 0), - (6, 0), - (7, 0), - (8, 20), - (9, 0), - (10, 13), - (11, 1), - (12, 1), - (13, 1), - (14, 0), - (15, 0), - (16, 0), - (17, 0), - (18, -1), - (19, -1), - (20, 0), - (21, 0), - (22, 0), - (23, 0), - (24, 10), - (25, 0), - (26, 0), - (27, 1), - (28, 0), - (29, 0), - (30, 0), - (31, 0), - (32, 0), - (33, 0), - (34, 0), - (35, 0), - (36, 0), - (37, 0), - (38, 0), - (39, 0), - (40, 0), - (41, 0), - (42, 0), - (43, 0), - (44, 0), - (45, 0), - (46, 0), - (47, 0), - (48, 0), - (49, 0), - (50, 0), - (51, 0), - (52, 0), - (53, 0), - (54, 0), - (55, 0), - (56, 0), - (57, 0), - (58, 0), - (59, 0), - (60, 0), - (61, 0), - (62, 0), - (63, 0), - (64, 0) -] - -for i in iparm - Pardiso.set_iparm!(solver, i...) -end -for i in Base.OneTo(length(iparm)) - @test Pardiso.get_iparm(solver, i) == iparm[i][2] + +# Testing and demonstrating Pardiso.set_iparm! for MKLPardisoSolver +for solver in solvers + iparm = [ + (1, 1), + (2, 2), + (3, 0), + (4, 0), + (5, 0), + (6, 0), + (7, 0), + (8, 20), + (9, 0), + (10, 13), + (11, 1), + (12, 1), + (13, 1), + (14, 0), + (15, 0), + (16, 0), + (17, 0), + (18, -1), + (19, -1), + (20, 0), + (21, 0), + (22, 0), + (23, 0), + (24, 10), + (25, 0), + (26, 0), + (27, 1), + (28, 0), + (29, 0), + (30, 0), + (31, 0), + (32, 0), + (33, 0), + (34, 0), + (35, 0), + (36, 0), + (37, 0), + (38, 0), + (39, 0), + (40, 0), + (41, 0), + (42, 0), + (43, 0), + (44, 0), + (45, 0), + (46, 0), + (47, 0), + (48, 0), + (49, 0), + (50, 0), + (51, 0), + (52, 0), + (53, 0), + (54, 0), + (55, 0), + (56, 0), + (57, 0), + (58, 0), + (59, 0), + (60, 0), + (61, 0), + (62, 0), + (63, 0), + (64, 0) + ] + + for i in iparm + Pardiso.set_iparm!(solver, i...) + end + + for i in Base.OneTo(length(iparm)) + @test Pardiso.get_iparm(solver, i) == iparm[i][2] + end end