From 6961555acbb2ccc7e515d58a2dea2cca33480ff6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Wed, 8 Feb 2023 21:03:18 +0100 Subject: [PATCH 1/4] Use AbstractSparseMatrixCSC for sparse factorizations --- .../src/LinearSolvePardiso.jl | 2 +- src/LinearSolve.jl | 1 + src/default.jl | 4 +- src/factorization.jl | 43 ++++++++++++------- 4 files changed, 31 insertions(+), 19 deletions(-) diff --git a/lib/LinearSolvePardiso/src/LinearSolvePardiso.jl b/lib/LinearSolvePardiso/src/LinearSolvePardiso.jl index 30ca59146..c2d8983b3 100644 --- a/lib/LinearSolvePardiso/src/LinearSolvePardiso.jl +++ b/lib/LinearSolvePardiso/src/LinearSolvePardiso.jl @@ -90,7 +90,7 @@ function LinearSolve.init_cacheval(alg::PardisoJL, A, b, u, Pl, Pr, maxiters::In Pardiso.set_iparm!(solver, 3, round(Int, abs(log10(reltol)), RoundDown) * 10 + 1) end - Pardiso.pardiso(solver, u, A, b) + Pardiso.pardiso(solver, u, SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A), b) return solver end diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index badd2a1ed..968667d20 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -9,6 +9,7 @@ using Base: cache_dependencies, Bool using LinearAlgebra using IterativeSolvers: Identity using SparseArrays +using SparseArrays: AbstractSparseMatrixCSC, nonzeros, rowvals, getcolptr using SciMLBase: AbstractLinearAlgorithm using SciMLOperators using SciMLOperators: AbstractSciMLOperator, IdentityOperator diff --git a/src/default.jl b/src/default.jl index ca529cd89..1c81cb3c2 100644 --- a/src/default.jl +++ b/src/default.jl @@ -50,13 +50,13 @@ function defaultalg(A::Diagonal, b, ::OperatorAssumptions{Nothing}) DiagonalFactorization() end -function defaultalg(A::SparseMatrixCSC{Tv, Ti}, b, +function defaultalg(A::AbstractSparseMatrixCSC{Tv, Ti}, b, ::OperatorAssumptions{true}) where {Tv, Ti} SparspakFactorization() end @static if INCLUDE_SPARSE - function defaultalg(A::SparseMatrixCSC{<:Union{Float64, ComplexF64}, Ti}, b, + function defaultalg(A::AbstractSparseMatrixCSC{<:Union{Float64, ComplexF64}, Ti}, b, ::OperatorAssumptions{true}) where {Ti} if length(b) <= 10_000 KLUFactorization() diff --git a/src/factorization.jl b/src/factorization.jl index 2349f45d1..eae6ea288 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -52,8 +52,8 @@ end function do_factorization(alg::LUFactorization, A, b, u) A = convert(AbstractMatrix, A) - if A isa SparseMatrixCSC - return lu(A) + if A isa AbstractSparseMatrixCSC + return lu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A))) else fact = lu!(A, alg.pivot) end @@ -277,7 +277,8 @@ function init_cacheval(alg::UMFPACKFactorization, A, b, u, Pl, Pr, maxiters::Int copy(nonzeros(A)), 0) finalizer(SuiteSparse.UMFPACK.umfpack_free_symbolic, res) else - return SuiteSparse.UMFPACK.UmfpackLU(A) + return SuiteSparse.UMFPACK.UmfpackLU(SparseMatrixCSC(size(A)..., getcolptr(A), + rowvals(A), nonzeros(A))) end end @@ -290,12 +291,15 @@ function SciMLBase.solve(cache::LinearCache, alg::UMFPACKFactorization; kwargs.. if alg.check_pattern && !(SuiteSparse.decrement(SparseArrays.getcolptr(A)) == cache.cacheval.colptr && SuiteSparse.decrement(SparseArrays.getrowval(A)) == cache.cacheval.rowval) - fact = lu(A) + fact = lu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), + nonzeros(A))) else - fact = lu!(cache.cacheval, A) + fact = lu!(cache.cacheval, + SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), + nonzeros(A))) end else - fact = lu(A) + fact = lu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A))) end cache = set_cacheval(cache, fact) end @@ -312,7 +316,8 @@ end function init_cacheval(alg::KLUFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) - return KLU.KLUFactorization(convert(AbstractMatrix, A)) # this takes care of the copy internally. + return KLU.KLUFactorization(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), + nonzeros(A))) end function SciMLBase.solve(cache::LinearCache, alg::KLUFactorization; kwargs...) @@ -323,21 +328,25 @@ function SciMLBase.solve(cache::LinearCache, alg::KLUFactorization; kwargs...) if alg.check_pattern && !(SuiteSparse.decrement(SparseArrays.getcolptr(A)) == cache.cacheval.colptr && SuiteSparse.decrement(SparseArrays.getrowval(A)) == cache.cacheval.rowval) - fact = KLU.klu(A) + fact = KLU.klu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), + nonzeros(A))) else # If we have a cacheval already, run umfpack_symbolic to ensure the symbolic factorization exists # This won't recompute if it does. KLU.klu_analyze!(cache.cacheval) - copyto!(cache.cacheval.nzval, A.nzval) + copyto!(cache.cacheval.nzval, nonzeros(A)) if cache.cacheval._numeric === C_NULL # We MUST have a numeric factorization for reuse, unlike UMFPACK. KLU.klu_factor!(cache.cacheval) end - fact = KLU.klu!(cache.cacheval, A) + fact = KLU.klu!(cache.cacheval, + SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), + nonzeros(A))) end else # New fact each time since the sparsity pattern can change # and thus it needs to reallocate - fact = KLU.klu(A) + fact = KLU.klu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), + nonzeros(A))) end cache = set_cacheval(cache, fact) end @@ -510,18 +519,20 @@ end function init_cacheval(::SparspakFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) - A = convert(AbstractMatrix, A) - return sparspaklu(A, factorize = false) + return sparspaklu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)), + factorize = false) end function SciMLBase.solve(cache::LinearCache, alg::SparspakFactorization; kwargs...) A = cache.A - A = convert(AbstractMatrix, A) if cache.isfresh if cache.cacheval !== nothing && alg.reuse_symbolic - fact = sparspaklu!(cache.cacheval, A) + fact = sparspaklu!(cache.cacheval, + SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), + nonzeros(A))) else - fact = sparspaklu(A) + fact = sparspaklu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), + nonzeros(A))) end cache = set_cacheval(cache, fact) end From bb3c690ef04a5b4018ff2c547ec0bb75ecb620d8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Wed, 8 Feb 2023 21:18:29 +0100 Subject: [PATCH 2/4] fix typo for pardiso --- lib/LinearSolvePardiso/Project.toml | 5 +-- .../src/LinearSolvePardiso.jl | 36 ++++++++++++++----- lib/LinearSolvePardiso/test/runtests.jl | 14 ++++---- 3 files changed, 37 insertions(+), 18 deletions(-) diff --git a/lib/LinearSolvePardiso/Project.toml b/lib/LinearSolvePardiso/Project.toml index ec899b0a7..dcdeffe81 100644 --- a/lib/LinearSolvePardiso/Project.toml +++ b/lib/LinearSolvePardiso/Project.toml @@ -8,12 +8,13 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [compat] -SciMLBase = "1.25" LinearSolve = "1.24" -Pardiso = "0.5" +Pardiso = "0.5" +SciMLBase = "1.25" UnPack = "1" julia = "1.6" diff --git a/lib/LinearSolvePardiso/src/LinearSolvePardiso.jl b/lib/LinearSolvePardiso/src/LinearSolvePardiso.jl index c2d8983b3..cb65b7fe9 100644 --- a/lib/LinearSolvePardiso/src/LinearSolvePardiso.jl +++ b/lib/LinearSolvePardiso/src/LinearSolvePardiso.jl @@ -1,14 +1,17 @@ module LinearSolvePardiso using Pardiso, LinearSolve, SciMLBase +using SparseArrays +using SparseArrays: nonzeros, rowvals, getcolptr + using UnPack Base.@kwdef struct PardisoJL <: LinearSolve.SciMLLinearSolveAlgorithm - nprocs::Union{Int, Nothing} = nothing - solver_type::Union{Int, Pardiso.Solver, Nothing} = nothing - matrix_type::Union{Int, Pardiso.MatrixType, Nothing} = nothing - iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing - dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing + nprocs::Union{Int,Nothing} = nothing + solver_type::Union{Int,Pardiso.Solver,Nothing} = nothing + matrix_type::Union{Int,Pardiso.MatrixType,Nothing} = nothing + iparm::Union{Vector{Tuple{Int,Int}},Nothing} = nothing + dparm::Union{Vector{Tuple{Int,Int}},Nothing} = nothing end MKLPardisoFactorize(; kwargs...) = PardisoJL(; solver_type = 0, kwargs...) @@ -17,9 +20,19 @@ LinearSolve.needs_concrete_A(alg::PardisoJL) = true # TODO schur complement functionality -function LinearSolve.init_cacheval(alg::PardisoJL, A, b, u, Pl, Pr, maxiters::Int, abstol, - reltol, verbose::Bool, - assumptions::LinearSolve.OperatorAssumptions) +function LinearSolve.init_cacheval( + alg::PardisoJL, + A, + b, + u, + Pl, + Pr, + maxiters::Int, + abstol, + reltol, + verbose::Bool, + assumptions::LinearSolve.OperatorAssumptions, +) @unpack nprocs, solver_type, matrix_type, iparm, dparm = alg A = convert(AbstractMatrix, A) @@ -90,7 +103,12 @@ function LinearSolve.init_cacheval(alg::PardisoJL, A, b, u, Pl, Pr, maxiters::In 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) + Pardiso.pardiso( + solver, + u, + SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)), + b, + ) return solver end diff --git a/lib/LinearSolvePardiso/test/runtests.jl b/lib/LinearSolvePardiso/test/runtests.jl index 0e26171e3..aa4a64f18 100644 --- a/lib/LinearSolvePardiso/test/runtests.jl +++ b/lib/LinearSolvePardiso/test/runtests.jl @@ -1,9 +1,11 @@ using LinearSolve, LinearSolvePardiso, SparseArrays, Random -A1 = sparse([1.0 0 -2 3 - 0 5 1 2 - -2 1 4 -7 - 3 2 -7 5]) +A1 = sparse([ + 1.0 0 -2 3 + 0 5 1 2 + -2 1 4 -7 + 3 2 -7 5 +]) b1 = rand(4) prob1 = LinearProblem(A1, b1) @@ -17,9 +19,7 @@ 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(), MKLPardisoIterate()) u = solve(prob1, alg; cache_kwargs...).u @test A1 * u ≈ b1 From 25559fe83e621b4f2c2f6e8af5702a391b69a66c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Wed, 8 Feb 2023 21:28:14 +0100 Subject: [PATCH 3/4] fix format --- .../src/LinearSolvePardiso.jl | 44 +++++++++---------- lib/LinearSolvePardiso/test/runtests.jl | 10 ++--- 2 files changed, 24 insertions(+), 30 deletions(-) diff --git a/lib/LinearSolvePardiso/src/LinearSolvePardiso.jl b/lib/LinearSolvePardiso/src/LinearSolvePardiso.jl index cb65b7fe9..ce4eb5fb0 100644 --- a/lib/LinearSolvePardiso/src/LinearSolvePardiso.jl +++ b/lib/LinearSolvePardiso/src/LinearSolvePardiso.jl @@ -7,11 +7,11 @@ using SparseArrays: nonzeros, rowvals, getcolptr using UnPack Base.@kwdef struct PardisoJL <: LinearSolve.SciMLLinearSolveAlgorithm - nprocs::Union{Int,Nothing} = nothing - solver_type::Union{Int,Pardiso.Solver,Nothing} = nothing - matrix_type::Union{Int,Pardiso.MatrixType,Nothing} = nothing - iparm::Union{Vector{Tuple{Int,Int}},Nothing} = nothing - dparm::Union{Vector{Tuple{Int,Int}},Nothing} = nothing + nprocs::Union{Int, Nothing} = nothing + solver_type::Union{Int, Pardiso.Solver, Nothing} = nothing + matrix_type::Union{Int, Pardiso.MatrixType, Nothing} = nothing + iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing + dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing end MKLPardisoFactorize(; kwargs...) = PardisoJL(; solver_type = 0, kwargs...) @@ -20,19 +20,17 @@ LinearSolve.needs_concrete_A(alg::PardisoJL) = true # TODO schur complement functionality -function LinearSolve.init_cacheval( - alg::PardisoJL, - A, - b, - u, - Pl, - Pr, - maxiters::Int, - abstol, - reltol, - verbose::Bool, - assumptions::LinearSolve.OperatorAssumptions, -) +function LinearSolve.init_cacheval(alg::PardisoJL, + A, + b, + u, + Pl, + Pr, + maxiters::Int, + abstol, + reltol, + verbose::Bool, + assumptions::LinearSolve.OperatorAssumptions) @unpack nprocs, solver_type, matrix_type, iparm, dparm = alg A = convert(AbstractMatrix, A) @@ -103,12 +101,10 @@ function LinearSolve.init_cacheval( 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, - ) + Pardiso.pardiso(solver, + u, + SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)), + b) return solver end diff --git a/lib/LinearSolvePardiso/test/runtests.jl b/lib/LinearSolvePardiso/test/runtests.jl index aa4a64f18..9681d0599 100644 --- a/lib/LinearSolvePardiso/test/runtests.jl +++ b/lib/LinearSolvePardiso/test/runtests.jl @@ -1,11 +1,9 @@ using LinearSolve, LinearSolvePardiso, SparseArrays, Random -A1 = sparse([ - 1.0 0 -2 3 - 0 5 1 2 - -2 1 4 -7 - 3 2 -7 5 -]) +A1 = sparse([1.0 0 -2 3 + 0 5 1 2 + -2 1 4 -7 + 3 2 -7 5]) b1 = rand(4) prob1 = LinearProblem(A1, b1) From 291b193c7d344dfe2450ec58b61a4dcd1149ba0a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Wed, 8 Feb 2023 22:18:19 +0100 Subject: [PATCH 4/4] re-establish all convert calls --- src/factorization.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/factorization.jl b/src/factorization.jl index eae6ea288..60a800411 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -316,6 +316,7 @@ end function init_cacheval(alg::KLUFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + A = convert(AbstractMatrix, A) return KLU.KLUFactorization(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A))) end @@ -519,6 +520,7 @@ end function init_cacheval(::SparspakFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + A = convert(AbstractMatrix, A) return sparspaklu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)), factorize = false) end