Merge pull request #43 from Wimmerer/savesymbolic

Initial pass at storing symbolic factorization

Author: ChrisRackauckas

Project.toml

@@ -15,6 +15,7 @@ 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"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
@@ -35,4 +36,4 @@ Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test","Pardiso"]
+test = ["Test", "Pardiso"]

src/LinearSolve.jl

@@ -11,7 +11,7 @@ using SciMLBase: AbstractDiffEqOperator, AbstractLinearAlgorithm
 using Setfield
 using UnPack
 using Requires
-
+using SuiteSparse
 # wrap
 import Krylov
 import KrylovKit # TODO
@@ -44,7 +44,7 @@ function __init__()
 end
 
 export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization,
-       RFLUFactorizaation
+       RFLUFactorizaation, UMFPACKFactorization
 export KrylovJL, KrylovJL_CG, KrylovJL_GMRES, KrylovJL_BICGSTAB, KrylovJL_MINRES,
        IterativeSolversJL, IterativeSolversJL_CG, IterativeSolversJL_GMRES,
        IterativeSolversJL_BICGSTAB, IterativeSolversJL_MINRES

src/factorization.jl

@@ -35,10 +35,11 @@ function init_cacheval(alg::LUFactorization, A, b, u)
 end
 
 # This could be a GenericFactorization perhaps?
-struct UMFPACKFactorization <: AbstractFactorization
+Base.@kwdef struct UMFPACKFactorization <: AbstractFactorization
+    reuse_symbolic::Bool = true
 end
 
-function init_cacheval(alg::UMFPACKFactorization, A, b, u)
+function init_cacheval(::UMFPACKFactorization, A, b, u)
     if A isa AbstractDiffEqOperator
         A = A.A
     end
@@ -49,6 +50,27 @@ function init_cacheval(alg::UMFPACKFactorization, A, b, u)
     end
 end
 
+function SciMLBase.solve(cache::LinearCache, alg::UMFPACKFactorization)
+    A = cache.A
+    if A isa AbstractDiffEqOperator
+        A = A.A
+    end
+    if cache.isfresh
+        if cache.cacheval !== nothing && alg.reuse_symbolic
+            # If we have a cacheval already, run umfpack_symbolic to ensure the symbolic factorization exists
+            # This won't recompute if it does.
+            SuiteSparse.UMFPACK.umfpack_symbolic!(cache.cacheval)
+            fact = lu!(cache.cacheval, A)
+        else
+            fact = init_cacheval(alg, A, cache.b, cache.u)
+        end
+        cache = set_cacheval(cache, fact)
+    end
+
+    y = ldiv!(cache.u, cache.cacheval, cache.b)
+    SciMLBase.build_linear_solution(alg,y,nothing)
+end
+
 ## QRFactorization
 
 struct QRFactorization{P} <: AbstractFactorization

test/runtests.jl

@@ -59,17 +59,33 @@ end
     y = solve(_prob)
     @test A1 *  y  ≈ b1
 
+    
     _prob = LinearProblem(sparse(A1.A), b1; u0=x1)
     y = solve(_prob)
     @test A1 *  y  ≈ b1
 end
 
+@testset "UMFPACK Factorization" begin
+    A1 = A/1; b1 = rand(n); x1 = zero(b)
+    A2 = A/2; b2 = rand(n); x2 = zero(b)
+
+    prob1 = LinearProblem(sparse(A1), b1; u0=x1)
+    prob2 = LinearProblem(sparse(A2), b2; u0=x2)
+    test_interface(UMFPACKFactorization(), prob1, prob2)
+
+    # Test that refactoring wrong throws.
+    cache = SciMLBase.init(prob1,UMFPACKFactorization(reuse_symbolic=true); cache_kwargs...) # initialize cache
+    y = solve(cache)
+    cache = LinearSolve.set_A(cache,sprand(n, n, 0.8))
+    @test_throws ArgumentError solve(cache)
+end
+
 @testset "Concrete Factorizations" begin
     for alg in (
                 LUFactorization(),
                 QRFactorization(),
                 SVDFactorization(),
-                RFLUFactorizaation(),
+                RFLUFactorizaation()
                )
         @testset "$alg" begin
             test_interface(alg, prob1, prob2)