From fe0d7905f5eb70a70c3e805da3cb963b3464279b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 1 Jun 2024 00:14:41 +0200 Subject: [PATCH 1/6] Try to fix #497 Pardiso defaults for highly indefinite matrices. This commit essentially reverts #89 and introduces a new kwarg "cache_analysis" (default `false`) to PardisoJL() which, if true would lead to the behaviour of #89. Also, allow the user to overwrite all iparms modified by the extension besides of 12. --- ext/LinearSolvePardisoExt.jl | 57 ++++++++++++++++++++---------------- src/extension_algs.jl | 5 +++- test/pardiso/pardiso.jl | 35 ++++++++++++++++++++-- 3 files changed, 68 insertions(+), 29 deletions(-) diff --git a/ext/LinearSolvePardisoExt.jl b/ext/LinearSolvePardisoExt.jl index 60312ce66..4cf5d3fcb 100644 --- a/ext/LinearSolvePardisoExt.jl +++ b/ext/LinearSolvePardisoExt.jl @@ -22,7 +22,7 @@ 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, cache_analysis, iparm, dparm = alg A = convert(AbstractMatrix, A) solver = if Pardiso.PARDISO_LOADED[] @@ -52,22 +52,6 @@ function LinearSolve.init_cacheval(alg::PardisoJL, end verbose && Pardiso.set_msglvl!(solver, Pardiso.MESSAGE_LEVEL_ON) - # pass in vector of tuples like [(iparm::Int, key::Int) ...] - if iparm !== nothing - for i in iparm - Pardiso.set_iparm!(solver, i...) - end - end - - if dparm !== nothing - for d in dparm - Pardiso.set_dparm!(solver, d...) - end - end - - # Make sure to say it's transposed because its CSC not CSR - Pardiso.set_iparm!(solver, 12, 1) - #= Note: It is recommended to use IPARM(11)=1 (scaling) and IPARM(13)=1 (matchings) for highly indefinite symmetric matrices e.g. from interior point optimizations or saddle point problems. @@ -79,10 +63,10 @@ function LinearSolve.init_cacheval(alg::PardisoJL, be changed to Pardiso.ANALYSIS_NUM_FACT in the solver loop otherwise instabilities occur in the example https://github.com/SciML/OrdinaryDiffEq.jl/issues/1569 =# - Pardiso.set_iparm!(solver, 11, 0) - Pardiso.set_iparm!(solver, 13, 0) - - Pardiso.set_phase!(solver, Pardiso.ANALYSIS) + if cache_analysis + Pardiso.set_iparm!(solver, 11, 0) + Pardiso.set_iparm!(solver, 13, 0) + end if alg.solver_type == 1 # PARDISO uses a numerical factorization A = LU for the first system and @@ -92,10 +76,30 @@ function LinearSolve.init_cacheval(alg::PardisoJL, Pardiso.set_iparm!(solver, 3, round(Int, abs(log10(reltol)), RoundDown) * 10 + 1) end - Pardiso.pardiso(solver, - u, - SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)), - b) + # pass in vector of tuples like [(iparm::Int, key::Int) ...] + if iparm !== nothing + for i in iparm + Pardiso.set_iparm!(solver, i...) + end + end + + if dparm !== nothing + for d in dparm + Pardiso.set_dparm!(solver, d...) + end + end + + # Make sure to say it's transposed because its CSC not CSR + # This is also the only value which should not be overwritten by users + Pardiso.set_iparm!(solver, 12, 1) + + if cache_analysis + Pardiso.set_phase!(solver, Pardiso.ANALYSIS) + Pardiso.pardiso(solver, + u, + SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)), + b) + end return solver end @@ -105,7 +109,8 @@ function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::PardisoJL; kwargs A = convert(AbstractMatrix, A) if cache.isfresh - Pardiso.set_phase!(cache.cacheval, Pardiso.NUM_FACT) + phase = alg.cache_analysis ? Pardiso.NUM_FACT : Pardiso.ANALYSIS_NUM_FACT + Pardiso.set_phase!(cache.cacheval, phase) Pardiso.pardiso(cache.cacheval, A, eltype(A)[]) cache.isfresh = false end diff --git a/src/extension_algs.jl b/src/extension_algs.jl index 4c2d4fecb..e59c984ac 100644 --- a/src/extension_algs.jl +++ b/src/extension_algs.jl @@ -154,12 +154,14 @@ 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} 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, Int}}, Nothing} = nothing) ext = Base.get_extension(@__MODULE__, :LinearSolvePardisoExt) @@ -170,7 +172,8 @@ 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, cache_analysis, iparm, dparm) end end end diff --git a/test/pardiso/pardiso.jl b/test/pardiso/pardiso.jl index 6753be21b..c84a849c6 100644 --- a/test/pardiso/pardiso.jl +++ b/test/pardiso/pardiso.jl @@ -1,4 +1,4 @@ -using LinearSolve, SparseArrays, Random +using LinearSolve, SparseArrays, Random, LinearAlgebra import Pardiso A1 = sparse([1.0 0 -2 3 @@ -13,12 +13,22 @@ n = 4 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) prob2 = LinearProblem(A2, b2) -for alg in (PardisoJL(), MKLPardisoFactorize(), MKLPardisoIterate()) +for alg in (PardisoJL(), MKLPardisoFactorize()) + 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 +end + +for alg in (MKLPardisoIterate(),) u = solve(prob1, alg; cache_kwargs...).u @test A1 * u ≈ b1 @@ -27,6 +37,8 @@ for alg in (PardisoJL(), MKLPardisoFactorize(), MKLPardisoIterate()) @test_broken A2 * u ≈ b2 end + + Random.seed!(10) A = sprand(n, n, 0.8); A2 = 2.0 .* A; @@ -53,6 +65,25 @@ sol33 = solve(linsolve) @test sol12.u ≈ sol32.u @test sol13.u ≈ sol33.u + +# Test for problem from #497 +function makeA() + n = 60 + colptr = [1, 4, 7, 11, 15, 17, 22, 26, 30, 34, 38, 40, 46, 50, 54, 58, 62, 64, 70, 74, 78, 82, 86, 88, 94, 98, 102, 106, 110, 112, 118, 122, 126, 130, 134, 136, 142, 146, 150, 154, 158, 160, 166, 170, 174, 178, 182, 184, 190, 194, 198, 202, 206, 208, 214, 218, 222, 224, 226, 228, 232] + rowval = [1, 3, 4, 1, 2, 4, 2, 4, 9, 10, 3, 5, 11, 12, 1, 3, 2, 4, 6, 11, 12, 2, 7, 9, 10, 2, 7, 8, 10, 8, 10, 15, 16, 9, 11, 17, 18, 7, 9, 2, 8, 10, 12, 17, 18, 8, 13, 15, 16, 8, 13, 14, 16, 14, 16, 21, 22, 15, 17, 23, 24, 13, 15, 8, 14, 16, 18, 23, 24, 14, 19, 21, 22, 14, 19, 20, 22, 20, 22, 27, 28, 21, 23, 29, 30, 19, 21, 14, 20, 22, 24, 29, 30, 20, 25, 27, 28, 20, 25, 26, 28, 26, 28, 33, 34, 27, 29, 35, 36, 25, 27, 20, 26, 28, 30, 35, 36, 26, 31, 33, 34, 26, 31, 32, 34, 32, 34, 39, 40, 33, 35, 41, 42, 31, 33, 26, 32, 34, 36, 41, 42, 32, 37, 39, 40, 32, 37, 38, 40, 38, 40, 45, 46, 39, 41, 47, 48, 37, 39, 32, 38, 40, 42, 47, 48, 38, 43, 45, 46, 38, 43, 44, 46, 44, 46, 51, 52, 45, 47, 53, 54, 43, 45, 38, 44, 46, 48, 53, 54, 44, 49, 51, 52, 44, 49, 50, 52, 50, 52, 57, 58, 51, 53, 59, 60, 49, 51, 44, 50, 52, 54, 59, 60, 50, 55, 57, 58, 50, 55, 56, 58, 56, 58, 57, 59, 55, 57, 50, 56, 58, 60] + nzval = [-0.64, 1.0, -1.0, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -1.0806825309567203, 1.0, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0] + A = SparseMatrixCSC(n, n, colptr, rowval, nzval) + return(A) +end + +A=makeA() +u0=fill(0.1,size(A,2)) +linprob = LinearProblem(A, A*u0) +u = LinearSolve.solve(linprob, PardisoJL(),verbose=true) +@test norm(u-u0) < 1.0e-14 + + + # Testing and demonstrating Pardiso.set_iparm! for MKLPardisoSolver solver = Pardiso.MKLPardisoSolver() iparm = [ From cab1f15a8e2fe74771fee68c76a7830056f7da0c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 1 Jun 2024 00:20:47 +0200 Subject: [PATCH 2/6] dparms should be Float64 --- src/extension_algs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/extension_algs.jl b/src/extension_algs.jl index e59c984ac..30db6e63d 100644 --- a/src/extension_algs.jl +++ b/src/extension_algs.jl @@ -163,7 +163,7 @@ struct PardisoJL{T1, T2} <: LinearSolve.SciMLLinearSolveAlgorithm matrix_type = nothing, cache_analysis = false, iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing, - dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing) + dparm::Union{Vector{Tuple{Int, Float64}}, Nothing} = nothing) ext = Base.get_extension(@__MODULE__, :LinearSolvePardisoExt) if ext === nothing error("PardisoJL requires that Pardiso is loaded, i.e. `using Pardiso`") From 2b75abb7b73045eaadfdf6847b40a443c963193a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 1 Jun 2024 00:28:19 +0200 Subject: [PATCH 3/6] update docstrings --- src/extension_algs.jl | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/src/extension_algs.jl b/src/extension_algs.jl index 30db6e63d..d60f6895b 100644 --- a/src/extension_algs.jl +++ b/src/extension_algs.jl @@ -86,6 +86,7 @@ end ```julia MKLPardisoFactorize(; 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) ``` @@ -98,7 +99,11 @@ A sparse factorization method using MKL Pardiso. ## Keyword Arguments -For the definition of the keyword arguments, see the Pardiso.jl documentation. +Setting `cache_analysis = true` disables Pardiso's scaling and matching defaults +and caches the result of the initial analysis phase for all further computations +with this solver. + +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. @@ -109,6 +114,7 @@ MKLPardisoFactorize(; kwargs...) = PardisoJL(; solver_type = 0, kwargs...) ```julia MKLPardisoIterate(; 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) ``` @@ -121,7 +127,11 @@ A mixed factorization+iterative method using MKL Pardiso. ## Keyword Arguments -For the definition of the keyword arguments, see the Pardiso.jl documentation. +Setting `cache_analysis = true` disables Pardiso's scaling and matching defaults +and caches the result of the initial analysis phase for all further computations +with this solver. + +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. @@ -133,6 +143,7 @@ MKLPardisoIterate(; kwargs...) = PardisoJL(; solver_type = 1, kwargs...) 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, Int}}, Nothing} = nothing) ``` @@ -145,6 +156,10 @@ A generic method using MKL Pardiso. Specifying `solver_type` is required. ## Keyword Arguments +Setting `cache_analysis = true` disables Pardiso's scaling and matching defaults +and caches the result of the initial analysis phase for all further computations +with this solver. + 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 From 2707d96ed48581018120ac1401f885a534e902b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 4 Jun 2024 18:59:07 +0200 Subject: [PATCH 4/6] Fix MKL Pardiso * iparm 12 -> 2 (transposed vs conjugate) * iparm 4 instead of 3 for reltol --- ext/LinearSolvePardisoExt.jl | 12 ++++++++---- test/pardiso/pardiso.jl | 2 +- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/ext/LinearSolvePardisoExt.jl b/ext/LinearSolvePardisoExt.jl index 4cf5d3fcb..1c323bdf9 100644 --- a/ext/LinearSolvePardisoExt.jl +++ b/ext/LinearSolvePardisoExt.jl @@ -25,7 +25,8 @@ function LinearSolve.init_cacheval(alg::PardisoJL, @unpack nprocs, solver_type, matrix_type, cache_analysis, iparm, dparm = alg A = convert(AbstractMatrix, A) - solver = if Pardiso.PARDISO_LOADED[] + transposed_iparm = 1 + solver = if false && Pardiso.PARDISO_LOADED[] solver = Pardiso.PardisoSolver() solver_type !== nothing && Pardiso.set_solver!(solver, solver_type) @@ -33,7 +34,9 @@ function LinearSolve.init_cacheval(alg::PardisoJL, else solver = Pardiso.MKLPardisoSolver() nprocs !== nothing && Pardiso.set_nprocs!(solver, nprocs) - + # for mkl 1 means conjugated an 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 end @@ -73,7 +76,8 @@ 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) + # Be aware that in the intel docs, iparm indexes are one lower. + Pardiso.set_iparm!(solver, 4, round(Int, abs(log10(reltol)), RoundDown) * 10 + 1) end # pass in vector of tuples like [(iparm::Int, key::Int) ...] @@ -91,7 +95,7 @@ function LinearSolve.init_cacheval(alg::PardisoJL, # Make sure to say it's transposed because its CSC not CSR # This is also the only value which should not be overwritten by users - Pardiso.set_iparm!(solver, 12, 1) + Pardiso.set_iparm!(solver, 12, transposed_iparm) if cache_analysis Pardiso.set_phase!(solver, Pardiso.ANALYSIS) diff --git a/test/pardiso/pardiso.jl b/test/pardiso/pardiso.jl index c84a849c6..aecda7c0a 100644 --- a/test/pardiso/pardiso.jl +++ b/test/pardiso/pardiso.jl @@ -34,7 +34,7 @@ for alg in (MKLPardisoIterate(),) u = solve(prob2, alg; cache_kwargs...).u @test eltype(u) <: Complex - @test_broken A2 * u ≈ b2 + @test A2 * u ≈ b2 end From d843f0188783b8de5200c5850fb2f1f5871c5b59 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 4 Jun 2024 23:51:10 +0200 Subject: [PATCH 5/6] Fix sequence of Pardiso calls. Panua Pardiso seems to depend on calling pardisoinit first. --- ext/LinearSolvePardisoExt.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/ext/LinearSolvePardisoExt.jl b/ext/LinearSolvePardisoExt.jl index 1c323bdf9..98db67a09 100644 --- a/ext/LinearSolvePardisoExt.jl +++ b/ext/LinearSolvePardisoExt.jl @@ -26,22 +26,24 @@ function LinearSolve.init_cacheval(alg::PardisoJL, A = convert(AbstractMatrix, A) transposed_iparm = 1 - solver = if false && Pardiso.PARDISO_LOADED[] + solver = if Pardiso.PARDISO_LOADED[] solver = Pardiso.PardisoSolver() + Pardiso.pardisoinit(solver) solver_type !== nothing && Pardiso.set_solver!(solver, solver_type) solver else solver = Pardiso.MKLPardisoSolver() + Pardiso.pardisoinit(solver) nprocs !== nothing && Pardiso.set_nprocs!(solver, nprocs) - # for mkl 1 means conjugated an 2 means transposed. + + # 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 end - Pardiso.pardisoinit(solver) # default initialization - if matrix_type !== nothing Pardiso.set_matrixtype!(solver, matrix_type) else @@ -118,9 +120,9 @@ function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::PardisoJL; kwargs Pardiso.pardiso(cache.cacheval, A, eltype(A)[]) cache.isfresh = false end - 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 From a33950961099344e19f9b0214b19b94c6277354c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 4 Jun 2024 23:52:04 +0200 Subject: [PATCH 6/6] remove verbosity from pardiso tests Otherwise Panua Pardiso may segfault during iteration for complex system --- test/pardiso/pardiso.jl | 27 +++++++++------------------ 1 file changed, 9 insertions(+), 18 deletions(-) diff --git a/test/pardiso/pardiso.jl b/test/pardiso/pardiso.jl index aecda7c0a..f05cca53b 100644 --- a/test/pardiso/pardiso.jl +++ b/test/pardiso/pardiso.jl @@ -13,30 +13,21 @@ n = 4 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) - prob2 = LinearProblem(A2, b2) -for alg in (PardisoJL(), MKLPardisoFactorize()) - 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 -end +cache_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiter = 30) -for alg in (MKLPardisoIterate(),) - u = solve(prob1, alg; cache_kwargs...).u - @test A1 * u ≈ b1 +for alg in (PardisoJL(), MKLPardisoFactorize(), MKLPardisoIterate()) + 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 + u = solve(prob2, alg; cache_kwargs...).u + @test eltype(u) <: Complex + @test A2 * u ≈ b2 end +return Random.seed!(10) @@ -79,7 +70,7 @@ end A=makeA() u0=fill(0.1,size(A,2)) linprob = LinearProblem(A, A*u0) -u = LinearSolve.solve(linprob, PardisoJL(),verbose=true) +u = LinearSolve.solve(linprob, PardisoJL()) @test norm(u-u0) < 1.0e-14