Interface to LinearSolve

Author: j-fu

Project.toml

@@ -1,21 +1,30 @@
 name = "ExtendableSparse"
 uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 authors = ["Juergen Fuhrmann <juergen.fuhrmann@wias-berlin.de>"]
-version = "0.6.8"
+version = "0.7.0"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
 DocStringExtensions = "0.8.0,0.9"
+LinearSolve = "1.23"
 Requires = "^1.1.3"
-julia = "^1.5"
+Setfield = "0.7, 0.8, 1"
+SciMLBase="1.49"
+UnPack = "1"
+julia = "^1.6"
 
 [extras]
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

src/ExtendableSparse.jl

@@ -2,19 +2,28 @@ module ExtendableSparse
 using SparseArrays
 using LinearAlgebra
 using SuiteSparse
+using Setfield: @set!
+using UnPack: @unpack
+import LinearSolve,SciMLBase
+
 using Requires
 
 using DocStringExtensions
 
 import SparseArrays: rowvals, getcolptr, nonzeros
-import Base: copy
+
 
 include("sparsematrixcsc.jl")
 include("sparsematrixlnk.jl")
 include("extendable.jl")
 
 export SparseMatrixLNK,ExtendableSparseMatrix,flush!,nnz, updateindex!, rawupdateindex!, colptrs
 
+
+include("linearsolve.jl")
+
+
+
 include("factorizations.jl")
 export JacobiPreconditioner, ILU0Preconditioner, ParallelJacobiPreconditioner, ParallelILU0Preconditioner, reorderlinsys
 export AbstractFactorization,LUFactorization, CholeskyFactorization

src/extendable.jl

@@ -62,6 +62,7 @@ Create similar but emtpy extendableSparseMatrix
 """
 Base.similar(m::ExtendableSparseMatrix{Tv,Ti}) where {Tv,Ti}=ExtendableSparseMatrix{Tv,Ti}(size(m)...)
 
+Base.similar(m::ExtendableSparseMatrix{Tv,Ti},::Type{T}) where {Tv,Ti,T}=ExtendableSparseMatrix{T,Ti}(size(m)...)
 
 """
 $(SIGNATURES)
@@ -394,10 +395,11 @@ function SparseArrays.dropzeros!(ext::ExtendableSparseMatrix)
 end
 
 
-function copy(S::ExtendableSparseMatrix)
+function Base.copy(S::ExtendableSparseMatrix)
     if isnothing(S.lnkmatrix)
         ExtendableSparseMatrix(copy(S.cscmatrix), nothing, S.phash)
     else
         ExtendableSparseMatrix(copy(S.cscmatrix), copy(S.lnkmatrix), S.phash)
     end
-end
+end
+

src/linearsolve.jl

@@ -0,0 +1,31 @@
+"""
+    LinearSolve.LinearProblem(A::ExtendableSparseMatrix,x,p=SciMLBase.NullParameters();u0=nothing,kwargs...)
+
+Create linear problem from ExtendableSparseMatrix. This uses the internal SparseMatrixCSC, and thus allows to access
+the functionality of [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)
+"""
+function LinearSolve.LinearProblem(A::ExtendableSparseMatrix,x,p=SciMLBase.NullParameters();u0=nothing,kwargs...)
+    flush!(A)
+    LinearSolve.LinearProblem{false}(A.cscmatrix,x,p;u0,kwargs...)
+end
+
+"""
+    LinearSolve.set_A(cache::LinearSolve.LinearCache, Aext::ExtendableSparseMatrix)
+
+Update the linear solve cache from an ExtendableSparseMatrix. Note that this update allows
+to take into account changes of the matrix pattern.
+"""
+function LinearSolve.set_A(cache::LinearSolve.LinearCache, Aext::ExtendableSparseMatrix)
+    flush!(Aext)
+
+    @unpack alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose = cache
+
+    if !pattern_equal(Aext.cscmatrix,A)
+        cacheval = LinearSolve.init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
+        @set! cache.cacheval=cacheval
+        @set! cache.isfresh = false
+    end
+    @set! cache.A = Aext.cscmatrix
+    @set! cache.isfresh = true
+    return cache
+end

src/sparsematrixcsc.jl

@@ -75,3 +75,11 @@ 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.
 
+
+"""
+    pattern_equal(a::SparseMatrixCSC,b::SparseMatrixCSC)
+
+Check if sparsity patterns of two SparseMatrixCSC objects are equal.
+This is generally faster than comparing hashes.
+"""
+pattern_equal(a::SparseMatrixCSC,b::SparseMatrixCSC) = a.colptr==b.colptr && a.rowval==b.rowval

src/sparsematrixlnk.jl

@@ -415,4 +415,4 @@ rowvals(S::SparseMatrixLNK) = getfield(S, :rowval)
 getcolptr(S::SparseMatrixLNK) = getfield(S, :colptr)
 nonzeros(S::SparseMatrixLNK) = getfield(S, :nzval)
 
-copy(S::SparseMatrixLNK) = SparseMatrixLNK(size(S, 1), size(S, 2), S.nnz, S.nentries, copy(getcolptr(S)), copy(rowvals(S)), copy(nonzeros(S)))
+Base.copy(S::SparseMatrixLNK) = SparseMatrixLNK(size(S, 1), size(S, 2), S.nnz, S.nentries, copy(getcolptr(S)), copy(rowvals(S)), copy(nonzeros(S)))

test/Project.toml

@@ -4,12 +4,13 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 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"
+AlgebraicMultigrid = "0.4,0.5"
 IncompleteLU = "^0.2"
 Pardiso = "^0.5.1"

test/runtests.jl

@@ -10,6 +10,8 @@ using AlgebraicMultigrid
 using IncompleteLU
 using IterativeSolvers
 
+using LinearSolve
+
 ##############################################################
 @testset "Constructors" begin
     function test_constructors()
@@ -320,7 +322,7 @@ end
     @test test_hermitian(300,:L)
 end
 
-function test_lu1(k,l,m; lufac=LUFactorization())
+function test_lu1(k,l,m; lufac=ExtendableSparse.LUFactorization())
     Acsc=fdrand(k,l,m,rand=()->1,matrixtype=SparseMatrixCSC)
     b=rand(k*l*m)
     LUcsc=lu(Acsc)
@@ -342,7 +344,7 @@ function test_lu1(k,l,m; lufac=LUFactorization())
     x1csc≈x1ext && x2csc ≈ x2ext
 end
 
-function test_lu2(k,l,m;lufac=LUFactorization())
+function test_lu2(k,l,m;lufac=ExtendableSparse.LUFactorization())
     Aext=fdrand(k,l,m,rand=()->1,matrixtype=ExtendableSparseMatrix)
     b=rand(k*l*m)
     lu!(lufac,Aext)
@@ -357,7 +359,7 @@ function test_lu2(k,l,m;lufac=LUFactorization())
 end
 
 
-@testset "LUFactorization" begin
+@testset "ExtendableSparse.LUFactorization" begin
     @test test_lu1(10,10,10)
     @test test_lu1(25,40,1)
     @test test_lu1(1000,1,1)
@@ -432,3 +434,30 @@ end
     @test test_parilu0(1000)
 end
 
+
+function test_linearsolve(n)
+ 
+    A=fdrand(n,1,1, matrixtype=ExtendableSparseMatrix)
+    b=rand(n)
+    c=A\b
+    @test c ≈ LinearSolve.solve(LinearProblem(A,b)).u
+    
+    A=fdrand(n,n,1, matrixtype=ExtendableSparseMatrix)
+    b=rand(n*n)
+    c=A\b
+    @test c ≈ LinearSolve.solve(LinearProblem(A,b)).u
+
+    @test c ≈ LinearSolve.solve(LinearProblem(A,b),UMFPACKFactorization()).u
+    @test c ≈ LinearSolve.solve(LinearProblem(A,b),KLUFactorization()).u
+    @test c ≈ LinearSolve.solve(LinearProblem(A,b),IterativeSolversJL_CG(),Pl=ILU0Preconditioner(A)).u
+    @test c ≈ LinearSolve.solve(LinearProblem(A,b),IterativeSolversJL_CG(),Pl=JacobiPreconditioner(A)).u
+    @test c ≈ LinearSolve.solve(LinearProblem(A,b),IterativeSolversJL_CG(),Pl=ParallelJacobiPreconditioner(A)).u
+    @test c ≈ LinearSolve.solve(LinearProblem(A,b),IterativeSolversJL_CG(),Pl=ILUTPreconditioner(A)).u
+    @test c ≈ LinearSolve.solve(LinearProblem(A,b),IterativeSolversJL_CG(),Pl=AMGPreconditioner(A)).u
+
+    
+end
+
+@testset "LinearSolve" begin
+    test_linearsolve(20)
+end