Merge pull request #36 from Wimmerer/UMFPACK_fix

Correct UMFPACK and test proper variables

Author: ChrisRackauckas

src/default.jl

@@ -30,7 +30,7 @@ function SciMLBase.solve(cache::LinearCache, alg::Nothing,
         alg = GenericFactorization(;fact_alg=ldlt!)
         SciMLBase.solve(cache, alg, args...; kwargs...)
     elseif A isa SparseMatrixCSC
-        alg = LUFactorization()
+        alg = UMFPACKFactorization()
         SciMLBase.solve(cache, alg, args...; kwargs...)
 
     # This catches the cases where a factorization overload could exist

src/factorization.jl

@@ -7,7 +7,7 @@ function SciMLBase.solve(cache::LinearCache, alg::AbstractFactorization)
     ldiv!(cache.u, cache.cacheval, cache.b)
 end
 
-## LUFactorization
+## LU Factorizations
 
 struct LUFactorization{P} <: AbstractFactorization
     pivot::P
@@ -33,6 +33,21 @@ function init_cacheval(alg::LUFactorization, A, b, u)
     return fact
 end
 
+# This could be a GenericFactorization perhaps?
+struct UMFPACKFactorization <: AbstractFactorization
+end
+
+function init_cacheval(alg::UMFPACKFactorization, A, b, u)
+    if A isa AbstractDiffEqOperator
+        A = A.A
+    end
+    if A isa SparseMatrixCSC
+        return lu(A)
+    else
+        error("Sparse LU is not defined for $(typeof(A))")
+    end
+end
+
 ## QRFactorization
 
 struct QRFactorization{P} <: AbstractFactorization

test/runtests.jl

@@ -42,23 +42,23 @@ end
     @test A1 *  y  ≈ b1
 
     _prob = LinearProblem(SymTridiagonal(A1.A), b1; u0=x1)
-    y = solve(prob1)
+    y = solve(_prob)
     @test A1 *  y  ≈ b1
 
     _prob = LinearProblem(Tridiagonal(A1.A), b1; u0=x1)
-    y = solve(prob1)
+    y = solve(_prob)
     @test A1 *  y  ≈ b1
 
     _prob = LinearProblem(Symmetric(A1.A), b1; u0=x1)
-    y = solve(prob1)
+    y = solve(_prob)
     @test A1 *  y  ≈ b1
 
     _prob = LinearProblem(Hermitian(A1.A), b1; u0=x1)
-    y = solve(prob1)
+    y = solve(_prob)
     @test A1 *  y  ≈ b1
 
     _prob = LinearProblem(sparse(A1.A), b1; u0=x1)
-    y = solve(prob1)
+    y = solve(_prob)
     @test A1 *  y  ≈ b1
 end