From 3129487a39ea180497c2b8aecfc25528c2af99c8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sun, 11 Apr 2021 01:22:48 +0200 Subject: [PATCH] Infrastructure for factorizations including Pardiso Infrastructure for factorizations including Pardiso * Reusable LU factorization * infrastructure for Pardiso, AlgebraicMultigrid, IncompleteLU etc. based on Require.jl * minmal julia version -> 1.5 (because of lu!) * Unittests for packages loaded via require * Documentation update * Define generic default methods for factorizations, specialize only if necessary. --- .github/workflows/ci.yml | 2 +- Project.toml | 7 +- docs/Project.toml | 4 +- docs/make.jl | 4 +- docs/src/changes.md | 9 +- docs/src/extsparse.md | 9 +- docs/src/internal.md | 15 +++ docs/src/iter.md | 26 ++++- src/ExtendableSparse.jl | 26 ++++- src/amg.jl | 26 +++++ src/extendable.jl | 61 +++------- src/factorizations.jl | 239 +++++++++++++++++++++++++++++++++++++++ src/ilu0.jl | 41 +++---- src/ilut.jl | 27 +++++ src/jacobi.jl | 29 +---- src/parallel_jacobi.jl | 30 +---- src/pardiso_lu.jl | 52 +++++++++ src/preconditioners.jl | 29 ----- src/sparsematrixcsc.jl | 11 ++ src/umfpack_lu.jl | 32 ++++++ test/Project.toml | 8 ++ test/runtests.jl | 87 +++++++++++++- 22 files changed, 609 insertions(+), 165 deletions(-) create mode 100644 docs/src/internal.md create mode 100644 src/amg.jl create mode 100644 src/factorizations.jl create mode 100644 src/ilut.jl create mode 100644 src/pardiso_lu.jl delete mode 100644 src/preconditioners.jl create mode 100644 src/umfpack_lu.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 5e44b6f..679e40c 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -15,7 +15,7 @@ jobs: fail-fast: false matrix: version: - - '1.3' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'. + - '1.5' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'. - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia. # - 'nightly' os: diff --git a/Project.toml b/Project.toml index 2e8ecaf..c666da2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,18 +1,21 @@ name = "ExtendableSparse" uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3" authors = ["Juergen Fuhrmann "] -version = "0.4.2" +version = "0.5.0" [deps] DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Requires = "ae029012-a4dd-5104-9daa-d747884805df" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] DocStringExtensions = "^0.8.0" -julia = "1" +Requires = "^1.1.3" +julia = "^1.5" [extras] Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" diff --git a/docs/Project.toml b/docs/Project.toml index db38c92..64bac97 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,7 +1,9 @@ [deps] +AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3" - +IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" +Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" diff --git a/docs/make.jl b/docs/make.jl index 151a499..1a74746 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,6 +1,5 @@ push!(LOAD_PATH,"../src/") -using Documenter, ExtendableSparse - +using Documenter, ExtendableSparse,Pardiso,AlgebraicMultigrid,IncompleteLU makedocs(sitename="ExtendableSparse.jl", modules = [ExtendableSparse], @@ -13,6 +12,7 @@ makedocs(sitename="ExtendableSparse.jl", "example.md", "extsparse.md", "iter.md", + "internal.md", "changes.md", ]) diff --git a/docs/src/changes.md b/docs/src/changes.md index 0f0d3b9..7200292 100644 --- a/docs/src/changes.md +++ b/docs/src/changes.md @@ -1,5 +1,12 @@ # Changes - +## v0.5, April 10, 2021 +- Introduce lu/lu! , factorize/factorize!, unifying LU factorizations and preconditioners +- Interface packages: Pardiso, AlgebraicMultigrid, IncompleteLU via Requires.jl +## v0.4, March 2021 +- Fix handling of Symmetrix matrices +- `rawupdateindex` does not check for entering zeros +- Compare with `COO` method +- Benchmarks in documentation ## v0.3.7, March 20, 2021 - Added parallel jacobi preconditioner (thanks, @jkr) - Fixes ldiv diff --git a/docs/src/extsparse.md b/docs/src/extsparse.md index 93dcf2a..a5d4f0d 100644 --- a/docs/src/extsparse.md +++ b/docs/src/extsparse.md @@ -1,7 +1,14 @@ # Sparse matrix handling +## Matrix creation and update API ```@autodocs Modules = [ExtendableSparse] -Pages = ["sparsematrixlnk.jl","sparsematrixcsc.jl","extendable.jl", "sprand.jl"] +Pages = ["extendable.jl"] +``` + +## Test matrix creation +```@autodocs +Modules = [ExtendableSparse] +Pages = ["sprand.jl"] ``` diff --git a/docs/src/internal.md b/docs/src/internal.md new file mode 100644 index 0000000..4d57068 --- /dev/null +++ b/docs/src/internal.md @@ -0,0 +1,15 @@ +# Internal API + + +## Linked List Sparse Matrix format +```@autodocs +Modules = [ExtendableSparse] +Pages = ["sparsematrixlnk.jl"] +``` + +## Some methods for SparseMatrixCSC +```@autodocs +Modules = [ExtendableSparse] +Pages = ["sparsematrixcsc.jl"] +``` + diff --git a/docs/src/iter.md b/docs/src/iter.md index cf67ddb..bf23225 100644 --- a/docs/src/iter.md +++ b/docs/src/iter.md @@ -1,5 +1,27 @@ -# Preconditioners +# Factorizations & Preconditioners + +## Factorizations + +In this package, preconditioners and LU factorizations are subcategories are both seen +as complete or approximate _factorizations_. Correspondingly there is a common API for +their creation. + +Factorizations from these package know the matrices which have been factorized. + ```@autodocs Modules = [ExtendableSparse] -Pages = ["preconditioners.jl", "ilu0.jl", "jacobi.jl", "parallel_jacobi.jl","simple_iteration.jl"] +Pages = ["factorizations.jl"] ``` + +## LU Factorizations +```@autodocs +Modules = [ExtendableSparse] +Pages = ["umfpack_lu.jl", "pardiso_lu.jl"] +``` + +## Preconditioners +```@autodocs +Modules = [ExtendableSparse] +Pages = ["jacobi.jl","ilu0.jl","parallel_jacobi.jl","ilut.jl","amg.jl"] +``` + diff --git a/src/ExtendableSparse.jl b/src/ExtendableSparse.jl index e6eb76d..f8bc062 100644 --- a/src/ExtendableSparse.jl +++ b/src/ExtendableSparse.jl @@ -1,19 +1,41 @@ module ExtendableSparse -using DocStringExtensions using SparseArrays using LinearAlgebra +using SuiteSparse +using Requires + +using DocStringExtensions + + include("sparsematrixcsc.jl") include("sparsematrixlnk.jl") include("extendable.jl") + export SparseMatrixLNK,ExtendableSparseMatrix,flush!,nnz, updateindex!, rawupdateindex!, colptrs -include("preconditioners.jl") +include("factorizations.jl") export JacobiPreconditioner, ILU0Preconditioner, ParallelJacobiPreconditioner +export issolver +export factorize,factorize!, update! +export ILUTPreconditioner, AMGPreconditioner include("simple_iteration.jl") export simple,simple! + + + include("sprand.jl") export sprand!,sprand_sdd!, fdrand,fdrand!,fdrand_coo + + +function __init__() + @require Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" include("pardiso_lu.jl") + @require IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" include("ilut.jl") + @require AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" include("amg.jl") +end + + + end # module diff --git a/src/amg.jl b/src/amg.jl new file mode 100644 index 0000000..316284c --- /dev/null +++ b/src/amg.jl @@ -0,0 +1,26 @@ +""" +$(TYPEDEF) + +""" +mutable struct AMGPreconditioner{Tv, Ti} <: AbstractExtendableSparsePreconditioner{Tv,Ti} + A::ExtendableSparseMatrix{Tv,Ti} + fact +end + +""" +``` +AMGPreconditioner(extmatrix) +``` +""" +function AMGPreconditioner(A::ExtendableSparseMatrix{Tv,Ti}) where {Tv,Ti} + @inbounds flush!(A) + p=AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(A.cscmatrix)) + AMGPreconditioner{Tv,Ti}(A,p) +end + +function update!(precon::AMGPreconditioner{Tv,Ti}) where {Tv,Ti} + @inbounds flush!(precon.A) + precon.fact=AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(precon.A.cscmatrix)) +end + + diff --git a/src/extendable.jl b/src/extendable.jl index 9b967c9..d54e48f 100644 --- a/src/extendable.jl +++ b/src/extendable.jl @@ -19,60 +19,46 @@ mutable struct ExtendableSparseMatrix{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv lnkmatrix::Union{SparseMatrixLNK{Tv,Ti},Nothing} """ - Time stamp of last pattern update + Pattern hash """ - pattern_timestamp::Float64 + phash::UInt64 end """ -$(SIGNATURES) +``` +ExtendableSparseMatrix(Tv,Ti,m,n) +ExtendableSparseMatrix(Tv,m,n) +ExtendableSparseMatrix(m,n) +``` +Create empty ExtendableSparseMatrix. This is equivalent to `spzeros(m,n)` for +`SparseMartrixCSC`. -Create empty ExtendableSparseMatrix. """ function ExtendableSparseMatrix{Tv,Ti}(m, n) where{Tv,Ti<:Integer} - ExtendableSparseMatrix{Tv,Ti}(spzeros(Tv,Ti,m,n),nothing,time()) + ExtendableSparseMatrix{Tv,Ti}(spzeros(Tv,Ti,m,n),nothing,0) end -""" -$(SIGNATURES) - -Create empty ExtendableSparseMatrix. -""" function ExtendableSparseMatrix(valuetype::Type{Tv},indextype::Type{Ti},m, n) where{Tv,Ti<:Integer} ExtendableSparseMatrix{Tv,Ti}(m,n) end -""" -$(SIGNATURES) - -Create empty ExtendablSparseMatrix. -This is a pendant to spzeros. -""" ExtendableSparseMatrix(valuetype::Type{Tv},m, n) where{Tv}=ExtendableSparseMatrix{Tv,Int}(m,n) - - -""" -$(SIGNATURES) - -Create empty ExtendableSparseMatrix. -This is a pendant to spzeros. -""" ExtendableSparseMatrix(m, n)=ExtendableSparseMatrix{Float64,Int}(m,n) - """ $(SIGNATURES) - Create ExtendableSparseMatrix from SparseMatrixCSC + Create ExtendableSparseMatrix from SparseMatrixCSC """ function ExtendableSparseMatrix(csc::SparseMatrixCSC{Tv,Ti}) where{Tv,Ti<:Integer} - return ExtendableSparseMatrix{Tv,Ti}(csc, nothing, time()) + return ExtendableSparseMatrix{Tv,Ti}(csc, nothing, phash(csc)) end """ $(SIGNATURES) - Create similar extendableSparseMatrix + +Create similar but emtpy extendableSparseMatrix """ Base.similar(m::ExtendableSparseMatrix{Tv,Ti}) where {Tv,Ti}=ExtendableSparseMatrix{Tv,Ti}(size(m)...) @@ -119,7 +105,6 @@ end $(SIGNATURES) Like [`updateindex!`](@ref) but without checking if v is zero. - """ function rawupdateindex!(ext::ExtendableSparseMatrix{Tv,Ti}, op,v, i,j) where {Tv,Ti<:Integer} k=findindex(ext.cscmatrix,i,j) @@ -141,7 +126,6 @@ $(SIGNATURES) Find index in CSC matrix and set value if it exists. Otherwise, set index in extension if `v` is nonzero. """ - function Base.setindex!(ext::ExtendableSparseMatrix{Tv,Ti}, v, i,j) where {Tv,Ti} k=findindex(ext.cscmatrix,i,j) if k>0 @@ -199,8 +183,6 @@ function Base.show(io::IO,::MIME"text/plain",ext::ExtendableSparseMatrix) end end - - """ $(SIGNATURES) @@ -211,7 +193,7 @@ function flush!(ext::ExtendableSparseMatrix) if ext.lnkmatrix!=nothing && nnz(ext.lnkmatrix)>0 ext.cscmatrix=ext.lnkmatrix+ext.cscmatrix ext.lnkmatrix=nothing - ext.pattern_timestamp=time() + ext.phash=phash(ext.cscmatrix) end return ext end @@ -273,17 +255,6 @@ function SparseArrays.findnz(ext::ExtendableSparseMatrix) end -""" -$(SIGNATURES) - -[`flush!`](@ref) and return LU factorization of ext.cscmatrix -""" -function LinearAlgebra.lu(ext::ExtendableSparseMatrix) - @inbounds flush!(ext) - return LinearAlgebra.lu(ext.cscmatrix) -end - - """ $(SIGNATURES) @@ -347,7 +318,7 @@ $(SIGNATURES) [`flush!`](@ref) and calculate norm from cscmatrix """ function LinearAlgebra.norm(A::ExtendableSparseMatrix, p::Real=2) - @time @inbounds flush!(A) + @inbounds flush!(A) return LinearAlgebra.norm(A.cscmatrix,p) end diff --git a/src/factorizations.jl b/src/factorizations.jl new file mode 100644 index 0000000..323a431 --- /dev/null +++ b/src/factorizations.jl @@ -0,0 +1,239 @@ +""" + $(TYPEDEF) + Abstract type for a factorization with ExtandableSparseMatrix. + + Any such preconditioner should have the following fields +```` + A + fact + phash +```` +and methods +```` + factorize(A; kwargs) + ldiv!(u, precon,v) + ldiv!(precon,v) + issolver(precon) + update!(precon) + factorize!(precon,A; kwargs) +```` + The idea is that, depending if the matrix pattern has changed, + different steps are needed to update the preconditioner. + + Moreover, they have the ExtendableSparseMatrix as a field, ensuring + consistency after construction. +""" +abstract type AbstractExtendableSparseFactorization{Tv, Ti} end + +""" +$(TYPEDEF) +Abstract subtype for preconditioners +""" +abstract type AbstractExtendableSparsePreconditioner{Tv, Ti} <:AbstractExtendableSparseFactorization{Tv, Ti} end + +""" +$(TYPEDEF) +Abstract subtype for (full) LU factorizations +""" +abstract type AbstractExtendableSparseLU{Tv, Ti} <:AbstractExtendableSparseFactorization{Tv, Ti} end + + + +""" +``` +issolver(factorization) +``` + +Determine if factorization is a solver or not +""" +issolver(::AbstractExtendableSparseLU)=true +issolver(::AbstractExtendableSparsePreconditioner)=false + + + +# +# Print default dict for interpolation into docstrings +# +function _myprint(dict::Dict{Symbol,Any}) + lines_out=IOBuffer() + for (k,v) in dict + println(lines_out," - $(k): $(v)") + end + String(take!(lines_out)) +end + + +""" +``` +const default_options +``` + +Default options for various factorizations: + +$(_myprint(default_options)) + +""" +const default_options=Dict{Symbol,Any}( + :kind => :umfpacklu, + :droptol => 1.0e-3, + :ensurelu => false, +) + + + +""" +``` +options(;kwargs...) +``` +Set default options and blend them with `kwargs` +""" +function options(;kwargs...) + opt=copy(default_options) + for (k,v) in kwargs + if haskey(opt,Symbol(k)) + opt[Symbol(k)]=v + end + end + opt +end + + +""" +``` +factorize(matrix) +factorize(matrix; kind=:umfpacklu) +``` +Default Julia LU factorization based on UMFPACK. + +``` +factorize(matrix; kind=:pardiso) +factorize(matrix; kind=:mklpardiso) +``` +LU factorization based on pardiso. For using this, you need to issue +`using Pardiso`. `:pardiso` uses the solver from [pardiso-project.org](https://pardiso-project.org), +while `:mklpardiso` uses the early 2000's fork in Intel's MKL library. + +``` +factorize(matrix; kind=:ilu0) +``` +Create the [`ILU0Preconditioner`](@ref) from this package. + +``` +factorize(matrix; kind=:jacobi) +``` +Create the [`JacobiPreconditioner`](@ref) from this package. + +``` +factorize(matrix; kind=:pjacobi) +``` +Create the [`ParallelJacobiPreconditioner`](@ref) from this package. + + +``` +factorize(matrix; kind=:ilut, droptol=1.0e-3) +``` +Create the [`ILUTPreconditioner`](@ref) wrapping the one +from [IncompleteLU.jl](https://github.com/haampie/IncompleteLU.jl) +For using this, you need to issue `using IncompleteLU` + + +``` +factorize(matrix; kind=:rsamg) +``` +Create the [`AMGPreconditioner`](@ref) wrapping the Ruge-Stüben AMG preconditioner from [AlgebraicMultigrid.jl](https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl) +""" +function factorize(A::ExtendableSparseMatrix; kwargs...) + opt=options(;kwargs...) + opt[:kind]==:umfpacklu && return ExtendableSparseUmfpackLU(A) + opt[:kind]==:pardiso && return PardisoLU(A,ps=Pardiso.PardisoSolver()) + opt[:kind]==:mklpardiso && return PardisoLU(A,ps=Pardiso.MKLPardisoSolver()) + if opt[:ensurelu] + error("Factorization $(opt[:kind]) is not an lu factorization") + end + opt[:kind]==:ilu0 && return ILU0Preconditioner(A) + opt[:kind]==:jacobi && return JacobiPreconditioner(A) + opt[:kind]==:pjacobi && return ParallelJacobiPreconditioner(A) + opt[:kind]==:ilut && return ILUTPreconditioner(A,droptol=opt[:droptol]) + opt[:kind]==:rsamg && return AMGPreconditioner(A) + error("Unknown factorization kind: $(opt[:kind])") +end + +""" +``` +factorize!(factorization_or_nothing, matrix; kwargs...) +``` + +Update factorization, possibly reusing information from the current state. +This method is aware of pattern changes. + +If `nothing` is passed as first parameter, [`factorize`](@ref) is called. +""" +factorize!(::Nothing, A::ExtendableSparseMatrix; kwargs...) = factorize(A; kwargs...) +function factorize!(p::AbstractExtendableSparseFactorization, A::ExtendableSparseMatrix; kwargs...) + p.A=A + update!(p) + p +end + + +""" +``` +lu(matrix) +lu(matrix,kind=:pardiso) +lu(matrix,kind=:mklpardiso) +``` +Wrapper for [`factorize`](@ref) restricted to lu factorizations. +""" +LinearAlgebra.lu(A::ExtendableSparseMatrix; kwargs...)= factorize(A; ensurelu=true, kwargs...) + + +""" +``` +lu!(factorization_or_nothing, matrix; kwargs...) +``` + +Update LU factorization, possibly reusing information from the current state. +This method is aware of pattern changes. + +If `nothing` is passed as first parameter, [`factorize`](@ref) is called. +""" +LinearAlgebra.lu!(::Nothing, A::ExtendableSparseMatrix; kwargs...)= factorize(A; kwargs...) +LinearAlgebra.lu!(lufact::AbstractExtendableSparseFactorization, A::ExtendableSparseMatrix; kwargs...)=factorize!(lufact,A;kwargs...) + +""" +``` + lufact\rhs +``` + +Solve LU factorization problem. +""" +Base.:\(lufact::AbstractExtendableSparseLU, v::AbstractArray{T,1} where T)=ldiv!(similar(v), lufact,v) + + +""" +``` +update!(factorization) +``` +Update factorization after matrix update. +""" +update!(::AbstractExtendableSparseFactorization) + + +""" +``` +ldiv!(u,factorization,v) +ldiv!(factorization,v) +``` + +Solve factorization. +""" +LinearAlgebra.ldiv!(u,fact::AbstractExtendableSparseFactorization, v)=ldiv!(u, fact.fact, v) +LinearAlgebra.ldiv!(fact::AbstractExtendableSparseFactorization, v)=ldiv!(fact.fact,v) + + + + +include("jacobi.jl") +include("ilu0.jl") +include("parallel_jacobi.jl") +include("umfpack_lu.jl") diff --git a/src/ilu0.jl b/src/ilu0.jl index 11a16c8..ee8555a 100644 --- a/src/ilu0.jl +++ b/src/ilu0.jl @@ -3,17 +3,19 @@ $(TYPEDEF) ILU(0) Preconditioner """ -mutable struct ILU0Preconditioner{Tv, Ti} <: AbstractExtendablePreconditioner{Tv,Ti} +mutable struct ILU0Preconditioner{Tv, Ti} <: AbstractExtendableSparsePreconditioner{Tv,Ti} extmatrix::ExtendableSparseMatrix{Tv,Ti} xdiag::Array{Tv,1} idiag::Array{Ti,1} - pattern_timestamp::Float64 + phash::UInt64 end -""" -$(SIGNATURES) -Constructor for ILU(0) preconditioner +""" +``` +ILU0Preconditioner(extsparse) +ILU0Preconditioner(cscmatrix) +``` """ function ILU0Preconditioner(extmatrix::ExtendableSparseMatrix{Tv,Ti}) where {Tv,Ti} @assert size(extmatrix,1)==size(extmatrix,2) @@ -25,19 +27,8 @@ function ILU0Preconditioner(extmatrix::ExtendableSparseMatrix{Tv,Ti}) where {Tv, update!(precon) end -""" -$(SIGNATURES) - -Constructor for ILU(0) preconditioner -""" ILU0Preconditioner(cscmatrix::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}=ILU0Preconditioner(ExtendableSparseMatrix(cscmatrix)) - -""" -$(SIGNATURES) - -Update ILU(0) preconditioner -""" function update!(precon::ILU0Preconditioner{Tv,Ti}) where {Tv,Ti} cscmatrix=precon.extmatrix.cscmatrix colptr=cscmatrix.colptr @@ -49,7 +40,7 @@ function update!(precon::ILU0Preconditioner{Tv,Ti}) where {Tv,Ti} # Find main diagonal index and # copy main diagonal values - if need_symbolic_update(precon) + if precon.phash != precon.extmatrix.phash @inbounds for j=1:n @inbounds for k=colptr[j]:colptr[j+1]-1 i=rowval[k] @@ -59,7 +50,7 @@ function update!(precon::ILU0Preconditioner{Tv,Ti}) where {Tv,Ti} end end end - timestamp!(precon) + precon.phash=precon.extmatrix.phash end @inbounds for j=1:n @@ -77,12 +68,14 @@ function update!(precon::ILU0Preconditioner{Tv,Ti}) where {Tv,Ti} precon end +function factorize!(precon::ILU0Preconditioner, A::ExtendableSparseMatrix; kwargs...) + flush!(A) + precon.extmatrix=A + update!(precon) + precon +end -""" -$(SIGNATURES) -Solve preconditioning system for ILU(0) -""" function LinearAlgebra.ldiv!(u::AbstractArray{T,1}, precon::ILU0Preconditioner, v::AbstractArray{T,1}) where T cscmatrix=precon.extmatrix.cscmatrix colptr=cscmatrix.colptr @@ -109,11 +102,7 @@ function LinearAlgebra.ldiv!(u::AbstractArray{T,1}, precon::ILU0Preconditioner, end end -""" -$(SIGNATURES) -Inplace solve of preconditioning system for ILU(0) -""" function LinearAlgebra.ldiv!(precon::ILU0Preconditioner, v::AbstractArray{T,1} where T) ldiv!(v, precon, v) end diff --git a/src/ilut.jl b/src/ilut.jl new file mode 100644 index 0000000..a71edca --- /dev/null +++ b/src/ilut.jl @@ -0,0 +1,27 @@ +""" +$(TYPEDEF) + +ILU(T) preconditioner +""" +mutable struct ILUTPreconditioner{Tv, Ti} <: AbstractExtendableSparsePreconditioner{Tv,Ti} + A::ExtendableSparseMatrix + fact::IncompleteLU.ILUFactorization{Tv,Ti} + droptol::Float64 +end + +""" +``` +ILUTPreconditioner(matrix; droptol=1.0e-3) +``` +""" +function ILUTPreconditioner(A::ExtendableSparseMatrix; droptol=1.0e-3) + @inbounds flush!(A) + ILUTPreconditioner(A,IncompleteLU.ilu(A.cscmatrix,τ=droptol),droptol) +end + +function update!(precon::ILUTPreconditioner, A::ExtendableSparseMatrix) + A=precon.A + @inbounds flush!(A) + precon.fact=IncompleteLU.ilu(A.cscmatrix,τ=precon.droptol) +end + diff --git a/src/jacobi.jl b/src/jacobi.jl index cb04693..b5c83ef 100644 --- a/src/jacobi.jl +++ b/src/jacobi.jl @@ -3,16 +3,11 @@ $(TYPEDEF) Jacobi preconditoner """ -struct JacobiPreconditioner{Tv, Ti} <: AbstractExtendablePreconditioner{Tv,Ti} +struct JacobiPreconditioner{Tv, Ti} <: AbstractExtendableSparsePreconditioner{Tv,Ti} extmatrix::ExtendableSparseMatrix{Tv,Ti} invdiag::Array{Tv,1} end -""" -$(SIGNATURES) - -Update Jacobi preconditoner -""" function update!(precon::JacobiPreconditioner) cscmatrix=precon.extmatrix.cscmatrix invdiag=precon.invdiag @@ -24,9 +19,10 @@ function update!(precon::JacobiPreconditioner) end """ -$(SIGNATURES) - -Construct Jacobi preconditoner +``` +JacobiPreconditioner(extmatrix) +JacobiPreconditioner(cscmatrix) +``` """ function JacobiPreconditioner(extmatrix::ExtendableSparseMatrix{Tv,Ti}) where {Tv,Ti} @assert size(extmatrix,1)==size(extmatrix,2) @@ -36,19 +32,9 @@ function JacobiPreconditioner(extmatrix::ExtendableSparseMatrix{Tv,Ti}) where {T update!(precon) end -""" -$(SIGNATURES) - -Construct Jacobi preconditoner -""" JacobiPreconditioner(cscmatrix::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}=JacobiPreconditioner(ExtendableSparseMatrix(cscmatrix)) -""" -$(SIGNATURES) - -Solve Jacobi preconditioning system -""" function LinearAlgebra.ldiv!(u::AbstractArray{T,1} where T, precon::JacobiPreconditioner, v::AbstractArray{T,1} where T) invdiag=precon.invdiag n=length(invdiag) @@ -57,11 +43,6 @@ function LinearAlgebra.ldiv!(u::AbstractArray{T,1} where T, precon::JacobiPreco end end -""" -$(SIGNATURES) - -Inplace solve Jacobi preconditioning system -""" function LinearAlgebra.ldiv!(precon::JacobiPreconditioner, v::AbstractArray{T,1} where T) ldiv!(v, precon, v) end diff --git a/src/parallel_jacobi.jl b/src/parallel_jacobi.jl index a7737c1..1af2a0f 100644 --- a/src/parallel_jacobi.jl +++ b/src/parallel_jacobi.jl @@ -3,16 +3,11 @@ $(TYPEDEF) Parallel Jacobi preconditioner """ -struct ParallelJacobiPreconditioner{Tv, Ti} <: AbstractExtendablePreconditioner{Tv,Ti} +struct ParallelJacobiPreconditioner{Tv, Ti} <: AbstractExtendableSparsePreconditioner{Tv,Ti} extmatrix::ExtendableSparseMatrix{Tv,Ti} invdiag::Array{Tv,1} end -""" -$(SIGNATURES) - -Parallel Jacobi preconditioner update -""" function update!(precon::ParallelJacobiPreconditioner) cscmatrix=precon.extmatrix.cscmatrix invdiag=precon.invdiag @@ -24,9 +19,10 @@ function update!(precon::ParallelJacobiPreconditioner) end """ -$(SIGNATURES) - -Parallel Jacobi preconditioner construtor +``` +ParallelJacobiPreconditioner(extmatrix) +ParallelJacobiPreconditioner(cscmatrix) +``` """ function ParallelJacobiPreconditioner(extmatrix::ExtendableSparseMatrix{Tv,Ti}) where {Tv,Ti} @assert size(extmatrix,1)==size(extmatrix,2) @@ -36,19 +32,8 @@ function ParallelJacobiPreconditioner(extmatrix::ExtendableSparseMatrix{Tv,Ti}) update!(precon) end -""" -$(SIGNATURES) - -Parallel Jacobi preconditioner construtor -""" ParallelJacobiPreconditioner(cscmatrix::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}=ParallelJacobiPreconditioner(ExtendableSparseMatrix(cscmatrix)) - -""" -$(SIGNATURES) - -Parallel Jacobi preconditioner solve -""" function LinearAlgebra.ldiv!(u::AbstractArray{T,1} where T, precon::ParallelJacobiPreconditioner, v::AbstractArray{T,1} where T) invdiag=precon.invdiag n=length(invdiag) @@ -57,11 +42,6 @@ function LinearAlgebra.ldiv!(u::AbstractArray{T,1} where T, precon::ParallelJac end end -""" -$(SIGNATURES) - -Parallel Jacobi preconditioner inplace solve -""" function LinearAlgebra.ldiv!(precon::ParallelJacobiPreconditioner, v::AbstractArray{T,1} where T) ldiv!(v, precon, v) end diff --git a/src/pardiso_lu.jl b/src/pardiso_lu.jl new file mode 100644 index 0000000..9851a49 --- /dev/null +++ b/src/pardiso_lu.jl @@ -0,0 +1,52 @@ +""" +$(TYPEDEF) + +LU Factorization +""" +mutable struct PardisoLU{Tv, Ti} <: AbstractExtendableSparseLU{Tv,Ti} + A::ExtendableSparseMatrix{Tv,Ti} + ps::Pardiso.AbstractPardisoSolver + phash::UInt64 +end + +""" +``` +PardisoLU(A; ps=Pardiso.MKLPardisoSolver) +``` +""" +function PardisoLU(A::ExtendableSparseMatrix{Tv,Ti};ps::Pardiso.AbstractPardisoSolver=Pardiso.MKLPardisoSolver()) where {Tv,Ti} + @inbounds flush!(A) + Acsc=A.cscmatrix + eltype(Acsc) == Float64 ? Pardiso.set_matrixtype!(ps, Pardiso.REAL_NONSYM) : Pardiso.set_matrixtype!(ps, Pardiso.COMPLEX_NONSYM) + Pardiso.pardisoinit(ps) + Pardiso.fix_iparm!(ps, :N) + Pardiso.set_phase!(ps, Pardiso.ANALYSIS_NUM_FACT) + Pardiso.pardiso(ps, Tv[], Acsc, Tv[]) + PardisoLU(A,ps,A.phash) +end + +function update!(lufact::PardisoLU{Tv,Ti}) where {Tv, Ti} + ps=lufact.ps + flush!(lufact.A) + Acsc=lufact.A.cscmatrix + if lufact.phash!=lufact.A.phash + Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL) + Pardiso.pardiso(ps, Tv[], Acsc, Tv[]) + Pardiso.set_phase!(ps, Pardiso.ANALYSIS_NUM_FACT) + lufact.phash=lufact.A.phash + else + Pardiso.set_phase!(ps, Pardiso.NUM_FACT) + end + Pardiso.pardiso(ps, Tv[], Acsc, Tv[]) + lufact +end + +function LinearAlgebra.ldiv!(u::AbstractArray{T,1} where T, lufact::PardisoLU, v::AbstractArray{T,1} where T) + ps=lufact.ps + Acsc=lufact.A.cscmatrix + Pardiso.set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE) + Pardiso.pardiso(ps, u, Acsc, v) + u +end + +LinearAlgebra.ldiv!(fact::PardisoLU, v::AbstractArray{T,1} where T)=ldiv!(v,fact,copy(v)) diff --git a/src/preconditioners.jl b/src/preconditioners.jl deleted file mode 100644 index 6f46159..0000000 --- a/src/preconditioners.jl +++ /dev/null @@ -1,29 +0,0 @@ -""" - Abstract type for an extendabe preconditoioner working together - with ExtandableSparseMatrix. - Still beta, so not in the published documentation. - - Any such preconditioner should have the following fields -```` - extmatrix - pattern_timestamp -```` -and methods -```` - update!(precon) -```` - The idea is that, depending if the matrix pattern has changed, - different steps are needed to update the preconditioner. - - Moreover, they have the ExtendableSparseMatrix as a field, ensuring - consistency after construction. -""" -abstract type AbstractExtendablePreconditioner{Tv,Ti} end - -need_symbolic_update(precon::AbstractExtendablePreconditioner)=precon.extmatrix.pattern_timestamp>precon.pattern_timestamp -timestamp!(precon::AbstractExtendablePreconditioner)= precon.pattern_timestamp=time() - -include("jacobi.jl") -include("ilu0.jl") -include("parallel_jacobi.jl") - diff --git a/src/sparsematrixcsc.jl b/src/sparsematrixcsc.jl index fade7b8..2397358 100644 --- a/src/sparsematrixcsc.jl +++ b/src/sparsematrixcsc.jl @@ -64,3 +64,14 @@ Trival flush! method for allowing to run the code with both `ExtendableSparseMat `SparseMatrixCSC`. """ flush!(csc::SparseMatrixCSC)=csc + + + +""" +$(SIGNATURES) + +Hash of csc matrix pattern. +""" +phash(csc::SparseMatrixCSC)=hash((hash(csc.colptr),hash(csc.rowval))) +# probably no good idea to add two hashes, so we hash them together. + diff --git a/src/umfpack_lu.jl b/src/umfpack_lu.jl new file mode 100644 index 0000000..7a5d8d6 --- /dev/null +++ b/src/umfpack_lu.jl @@ -0,0 +1,32 @@ +""" +$(TYPEDEF) + +Default Julia LU Factorization based on umfpack. +""" +mutable struct ExtendableSparseUmfpackLU{Tv, Ti} <: AbstractExtendableSparseLU{Tv,Ti} + A::ExtendableSparseMatrix{Tv,Ti} + fact::SuiteSparse.UMFPACK.UmfpackLU{Tv,Ti} + phash::UInt64 +end + +""" +``` +ExtendableSparseUmfpackLU(A) +``` +""" +function ExtendableSparseUmfpackLU(A::ExtendableSparseMatrix{Tv,Ti}) where {Tv,Ti} + flush!(A) + ExtendableSparseUmfpackLU(A,lu(A.cscmatrix),A.phash) +end + +function update!(lufact::ExtendableSparseUmfpackLU) + flush!(lufact.A) + if lufact.A.phash!=lufact.phash + lufact.fact=lu(lufact.A.cscmatrix) + lufact.phash=lufact.A.phash + else + lufact.fact=lu!(lufact.fact,lufact.A.cscmatrix) + end + lufact +end + diff --git a/test/Project.toml b/test/Project.toml index 17d9c79..623b4f3 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,14 @@ [deps] +AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[compat] +AlgebraicMultigrid = "^0.4" +IncompleteLU = "^0.2" +Pardiso = "^0.5.1" diff --git a/test/runtests.jl b/test/runtests.jl index 838a7f2..c458874 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,6 +5,10 @@ using ExtendableSparse using Printf using BenchmarkTools +using Pardiso +using AlgebraicMultigrid +using IncompleteLU + ############################################################## @testset "Constructors" begin function test_constructors() @@ -41,9 +45,9 @@ end end ################################################################# function test_timing(k,l,m) - t1=@belapsed fdrand($k,$l,$m,matrixtype=$SparseMatrixCSC) seconds=1 - t2=@belapsed fdrand($k,$l,$m,matrixtype=$ExtendableSparseMatrix) seconds=1 - t3=@belapsed fdrand($k,$l,$m,matrixtype=$SparseMatrixLNK) seconds=1 + t1=@belapsed fdrand($k,$l,$m,matrixtype=$SparseMatrixCSC) seconds=0.1 + t2=@belapsed fdrand($k,$l,$m,matrixtype=$ExtendableSparseMatrix) seconds=0.1 + t3=@belapsed fdrand($k,$l,$m,matrixtype=$SparseMatrixLNK) seconds=0.1 @printf("CSC: %.4f EXT: %.4f LNK: %.4f\n",t1*1000,t2*1000,t3*1000 ) t3 x<1,r[end-100:end]./r[end-101:end-1]),norm(it-exact) + nr=length(r) + tail=min(100,length(r)÷2) + all(x-> x<1,r[end-tail:end]./r[end-tail-1:end-1]),norm(it-exact) end @@ -226,6 +232,8 @@ end @test all(isapprox.(test_precon(ILU0Preconditioner,20,20,20), (true, 1.3535160424212675e-5), rtol=1.0e-5)) @test all(isapprox.(test_precon(JacobiPreconditioner,20,20,20), (true, 2.0406032775945658e-5), rtol=1.0e-5)) @test all(isapprox.(test_precon(ParallelJacobiPreconditioner,20,20,20), (true, 2.0406032775945658e-5), rtol=1.0e-5)) + @test all(isapprox.(test_precon(ILUTPreconditioner,20,20,20), (true, 1.2719511868322086e-5), rtol=1.0e-5)) + @test all(isapprox.(test_precon(AMGPreconditioner,20,20,20), (true, 6.863753664354144e-7), rtol=1.0e-2)) end @@ -277,3 +285,74 @@ end +function test_lu(k,l,m; kind=:umfpacklu) + Acsc=fdrand(k,l,m,rand=()->1,matrixtype=SparseMatrixCSC) + b=rand(k*l*m) + LUcsc=lu(Acsc) + x1csc=LUcsc\b + for i=1:k*l*m + Acsc[i,i]+=1.0 + end + LUcsc=lu!(LUcsc,Acsc) + x2csc=LUcsc\b + + Aext=fdrand(k,l,m,rand=()->1,matrixtype=ExtendableSparseMatrix) + LUext=lu(Aext,kind=kind) + x1ext=LUext\b + for i=1:k*l*m + Aext[i,i]+=1.0 + end + update!(LUext) + x2ext=LUext\b + x1csc≈x1ext && x2csc ≈ x2ext +end + +function test_lupattern1(k,l,m; kind=:umfpacklu) + Aext=fdrand(k,l,m,rand=()->1,matrixtype=ExtendableSparseMatrix) + b=rand(k*l*m) + LUext=lu(Aext,kind=kind) + x1ext=LUext\b + for i=1:k*l*m-3 + Aext[i,i+3]-=1.0e-4 + end + LUext=lu!(LUext,Aext) + x2ext=LUext\b + all(x1ext.< x2ext) +end + +function test_lupattern2(k,l,m) + Aext=fdrand(k,l,m,rand=()->1,matrixtype=ExtendableSparseMatrix) + b=rand(k*l*m) + LUext=lu(Aext) + x1ext=LUext\b + for i=1:k*l*m-3 + Aext[i,i+3]-=1.0e-4 + end + update!(LUext) + x2ext=LUext\b + all(x1ext.< x2ext) +end + +@testset "lu!+update!" begin + @test test_lu(10,10,10) + @test test_lu(25,40,1) + @test test_lu(1000,1,1) + + @test test_lupattern1(10,10,10) + @test test_lupattern1(25,40,1) + @test test_lupattern1(1000,1,1) + + @test test_lupattern2(10,10,10) + @test test_lupattern2(25,40,1) + @test test_lupattern2(1000,1,1) +end + +@testset "pardiso" begin + @test test_lu(10,10,10,kind=:mklpardiso) + @test test_lu(25,40,1,kind=:mklpardiso) + @test test_lu(1000,1,1,kind=:mklpardiso) + + @test test_lupattern1(10,10,10,kind=:mklpardiso) + @test test_lupattern1(25,40,1,kind=:mklpardiso) + @test test_lupattern1(1000,1,1,kind=:mklpardiso) +end