Add sparspak lu factorization (via require)

Author: j-fu

Project.toml

@@ -1,13 +1,11 @@
 name = "ExtendableSparse"
 uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 authors = ["Juergen Fuhrmann <juergen.fuhrmann@wias-berlin.de>"]
-version = "0.7.0"
+version = "0.8.0"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -16,7 +14,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 DocStringExtensions = "0.8.0,0.9"
-IterativeSolvers = "0.9.2"
 Requires = "1.1.3"
 julia = "1.6"
 

src/ExtendableSparse.jl

@@ -1,7 +1,10 @@
 module ExtendableSparse
 using SparseArrays
 using LinearAlgebra
-using SuiteSparse
+
+if Base.USE_GPL_LIBS
+    using SuiteSparse
+end
 
 using Requires
 
@@ -23,7 +26,7 @@ export AbstractFactorization,LUFactorization, CholeskyFactorization
 export issolver
 export factorize!, update!
 export ILUTPreconditioner, AMGPreconditioner
-export PardisoLU, MKLPardisoLU
+export PardisoLU, MKLPardisoLU,SparspakLU
 
 include("simple_iteration.jl")
 export simple,simple!
@@ -39,6 +42,7 @@ 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")
+    @require Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac" include("sparspak.jl")
 end
 
 

src/sparspak.jl

@@ -0,0 +1,42 @@
+mutable struct SparspakLU{Tv, Ti} <: AbstractLUFactorization{Tv,Ti} 
+    A::Union{ExtendableSparseMatrix{Tv,Ti},Nothing}
+    P::Union{Sparspak.SpkProblem.Problem{Ti,Tv},Nothing}
+    S::Union{Sparspak.SparseSolver.SparseSolver{Ti,Tv},Nothing}
+    phash::UInt64
+end
+
+function SparspakLU{Tv,Ti}() where {Tv,Ti}
+    fact=SparspakLU{Tv,Ti}(nothing,nothing,nothing,0)
+end
+
+
+"""
+```
+SparspakLU(;valuetype=Float64) 
+```
+
+LU factorization based on P.Kryls's re-implementation of SPARSPAK.
+"""
+SparspakLU(;valuetype::Type=Float64, indextype::Type=Int64,kwargs...)=SparspakLU{valuetype,indextype}(;kwargs...)
+
+function update!(lufact::SparspakLU{Tv,Ti}) where {Tv, Ti}
+    flush!(lufact.A)
+    lufact.P=Sparspak.SpkProblem.Problem(size(lufact.A)...)
+    Sparspak.SpkProblem.insparse!(lufact.P,lufact.A)
+    lufact.S=Sparspak.SparseSolver.SparseSolver(lufact.P)
+end
+
+function LinearAlgebra.ldiv!(u::AbstractArray{T,1} where T, lufact::SparspakLU, v::AbstractArray{T,1} where T)
+    u.=ldiv!(lufact,v)
+end
+
+function LinearAlgebra.ldiv!(lufact::SparspakLU, v::AbstractArray{T,1} where T)
+    Sparspak.SpkProblem.infullrhs!(lufact.P,v)
+    Sparspak.SparseSolver.solve!(lufact.S)
+    lufact.P.x
+end
+
+
+@eval begin
+    @makefrommatrix SparspakLU
+end

test/Project.toml

@@ -7,9 +7,12 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 AlgebraicMultigrid = "0.4,0.5"
 IncompleteLU = "^0.2"
+IterativeSolvers = "0.9"
 Pardiso = "^0.5.1"
+Sparspak = "0.2"

test/runtests.jl

@@ -9,6 +9,7 @@ using Pardiso
 using AlgebraicMultigrid
 using IncompleteLU
 using IterativeSolvers
+using Sparspak
 
 
 ##############################################################
@@ -403,26 +404,21 @@ end
 end
 
 
-if Base.USE_GPL_LIBS
-#requires SuiteSparse which is not available on non-GPL builds
 @testset "mkl-pardiso" begin
-if !Sys.isapple()
-    @test test_lu1(10,10,10,lufac=MKLPardisoLU())
-    @test test_lu1(25,40,1,lufac=MKLPardisoLU())
-    @test test_lu1(1000,1,1,lufac=MKLPardisoLU())
-
-    @test test_lu2(10,10,10,lufac=MKLPardisoLU())
-    @test test_lu2(25,40,1,lufac=MKLPardisoLU())
-    @test test_lu2(1000,1,1,lufac=MKLPardisoLU())
-end
+    if !Sys.isapple()
+        @test test_lu1(10,10,10,lufac=MKLPardisoLU())
+        @test test_lu1(25,40,1,lufac=MKLPardisoLU())
+        @test test_lu1(1000,1,1,lufac=MKLPardisoLU())
+        
+        @test test_lu2(10,10,10,lufac=MKLPardisoLU())
+        @test test_lu2(25,40,1,lufac=MKLPardisoLU())
+        @test test_lu2(1000,1,1,lufac=MKLPardisoLU())
+    end
     # @test test_lu2(10,10,10,lufac=MKLPardisoLU(mtype=2))
     # @test test_lu2(25,40,1,lufac=MKLPardisoLU(mtype=2))
     # @test test_lu2(1000,1,1,lufac=MKLPardisoLU(mtype=2))
 end
-end
 
-if Base.USE_GPL_LIBS
-#requires SuiteSparse which is not available on non-GPL builds
 if Pardiso.PARDISO_LOADED[]
     @testset "pardiso" begin
         @test test_lu1(10,10,10,lufac=PardisoLU())
@@ -439,6 +435,17 @@ if Pardiso.PARDISO_LOADED[]
 
     end
 end
+
+
+@testset "sparspak" begin
+    @test test_lu1(10,10,10,lufac=SparspakLU())
+    @test test_lu1(25,40,1,lufac=SparspakLU())
+    @test test_lu1(1000,1,1,lufac=SparspakLU())
+    
+    @test test_lu2(10,10,10,lufac=SparspakLU())
+    @test test_lu2(25,40,1,lufac=SparspakLU())
+    @test test_lu2(1000,1,1,lufac=SparspakLU())
+    
 end