Merge pull request #44 from Wimmerer/KLU

Add KLU factorization

Author: ChrisRackauckas

Project.toml

@@ -6,6 +6,7 @@ version = "0.1.4"
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
+KLU = "ef3ab10e-7fda-4108-b977-705223b18434"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -21,6 +22,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 [compat]
 ArrayInterface = "3"
 IterativeSolvers = "0.9.2"
+KLU = "0.2.1"
 Krylov = "0.7.9"
 KrylovKit = "0.5"
 RecursiveFactorization = "0.2"

src/LinearSolve.jl

@@ -12,6 +12,7 @@ using Setfield
 using UnPack
 using Requires
 using SuiteSparse
+using KLU
 # wrap
 import Krylov
 import KrylovKit # TODO
@@ -44,7 +45,7 @@ function __init__()
 end
 
 export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization,
-       RFLUFactorizaation, UMFPACKFactorization
+       RFLUFactorizaation, UMFPACKFactorization, KLUFactorization
 export KrylovJL, KrylovJL_CG, KrylovJL_GMRES, KrylovJL_BICGSTAB, KrylovJL_MINRES,
        IterativeSolversJL, IterativeSolversJL_CG, IterativeSolversJL_GMRES,
        IterativeSolversJL_BICGSTAB, IterativeSolversJL_MINRES

src/factorization.jl

@@ -71,6 +71,42 @@ function SciMLBase.solve(cache::LinearCache, alg::UMFPACKFactorization)
     SciMLBase.build_linear_solution(alg,y,nothing)
 end
 
+Base.@kwdef struct KLUFactorization <: AbstractFactorization
+    reuse_symbolic::Bool = true
+end
+
+function init_cacheval(::KLUFactorization, A, b, u)
+    if A isa AbstractDiffEqOperator
+        A = A.A
+    end
+    if A isa SparseMatrixCSC
+        return klu(A)
+    else
+        error("KLU is not defined for $(typeof(A))")
+    end
+end
+
+function SciMLBase.solve(cache::LinearCache, alg::KLUFactorization)
+    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.
+            KLU.klu_analyze!(cache.cacheval)
+            fact = klu!(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

@@ -74,12 +74,30 @@ end
     test_interface(UMFPACKFactorization(), prob1, prob2)
 
     # Test that refactoring wrong throws.
-    cache = SciMLBase.init(prob1,UMFPACKFactorization(reuse_symbolic=true); cache_kwargs...) # initialize cache
+    cache = SciMLBase.init(prob1,UMFPACKFactorization(); cache_kwargs...) # initialize cache
     y = solve(cache)
     cache = LinearSolve.set_A(cache,sprand(n, n, 0.8))
     @test_throws ArgumentError solve(cache)
 end
 
+@testset "KLU 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(KLUFactorization(), prob1, prob2)
+
+    # Test that refactoring wrong throws.
+    cache = SciMLBase.init(prob1,KLUFactorization(); cache_kwargs...) # initialize cache
+    y = solve(cache)
+    X = copy(A1)
+    X[8,8] = 0.0
+    X[7,8] = 1.0
+    cache = LinearSolve.set_A(cache,sparse(X))
+    @test_throws ArgumentError solve(cache)
+end
+
 @testset "Concrete Factorizations" begin
     for alg in (
                 LUFactorization(),