Using ArrayInterface.lu_instance and added tests

Author: mauro3

Project.toml

@@ -22,7 +22,6 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 
 [compat]
 ArrayInterface = "1.1, 2.0"

src/OrdinaryDiffEq.jl

@@ -5,7 +5,7 @@ module OrdinaryDiffEq
 
   using Logging
 
-  using MuladdMacro, SparseArrays, SuiteSparse
+  using MuladdMacro, SparseArrays
 
   using LinearAlgebra
 

src/derivative_utils.jl

@@ -538,23 +538,12 @@ function build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,::Val{IIP}) where IIP
     elseif IIP
       similar(J)
     else
-      W = if u isa StaticArray
-      lu(J)
-      elseif u isa Number
-        u
+      if f.jac_prototype===nothing
+        LU{LinearAlgebra.lutype(uEltypeNoUnits)}(Matrix{uEltypeNoUnits}(undef, 0, 0),
+                                                 Vector{LinearAlgebra.BlasInt}(undef, 0),
+                                                 zero(LinearAlgebra.BlasInt))
       else
-        # TODO: make this more general, maybe by running calc_W and use its returned value
-        if f.jac_prototype isa SparseMatrixCSC
-          SuiteSparse.UMFPACK.UmfpackLU(Ptr{Cvoid}(), Ptr{Cvoid}(), 1, 1,
-                                        f.jac_prototype.colptr[1:1],
-                                        f.jac_prototype.rowval[1:1],
-                                        f.jac_prototype.nzval[1:1],
-                                        0)
-        else
-          LU{LinearAlgebra.lutype(uEltypeNoUnits)}(Matrix{uEltypeNoUnits}(undef, 0, 0),
-                                                   Vector{LinearAlgebra.BlasInt}(undef, 0),
-                                                   zero(LinearAlgebra.BlasInt))
-        end
+        ArrayInterface.lu_instance(f.jac_prototype)
       end
     end # end W
   end

test/interface/sparsediff_tests.jl

@@ -1,17 +1,30 @@
-using Test
+  using Test
 using OrdinaryDiffEq
 using SparseArrays
 using LinearAlgebra
 
+## in-place
 #https://github.com/JuliaDiffEq/SparseDiffTools.jl/blob/master/test/test_integration.jl
-function f(dx,x,p,t)
-    for i in 2:length(x)-1
-      dx[i] = x[i-1] - 2x[i] + x[i+1]
-    end
-    dx[1] = -2x[1] + x[2]
-    dx[end] = x[end-1] - 2x[end]
-    nothing
+function f_ip(dx,x,p,t)
+  for i in 2:length(x)-1
+    dx[i] = x[i-1] - 2x[i] + x[i+1]
+  end
+  dx[1] = -2x[1] + x[2]
+  dx[end] = x[end-1] - 2x[end]
+  nothing
+end
+
+## out-of-place
+function f_oop(x,p,t)
+  dx = similar(x)
+  for i in 2:length(x)-1
+    dx[i] = x[i-1] - 2x[i] + x[i+1]
   end
+  dx[1] = -2x[1] + x[2]
+  dx[end] = x[end-1] - 2x[end]
+  return dx
+end
+
 
 function second_derivative_stencil(N)
   A = zeros(N,N)
@@ -34,22 +47,29 @@ jac_sp = sparse(generate_sparsity_pattern(10))
 colors = repeat(1:3,10)[1:10]
 u0=[1.,2.,3,4,5,5,4,3,2,1]
 tspan=(0.,10.)
-odefun_sp= ODEFunction(f,colorvec=colors,jac_prototype=jac_sp)
-prob_sp = ODEProblem(odefun_sp,u0,tspan)
-prob_std = ODEProblem(f,u0,tspan)
-
-sol_sp=solve(prob_sp,Rodas5(autodiff=false),abstol=1e-10,reltol=1e-10)
-@test sol_sp.retcode==:Success#test sparse finitediff
-sol=solve(prob_std,Rodas5(autodiff=false),abstol=1e-10,reltol=1e-10)
-@test sol_sp.u[end]≈sol.u[end] atol=1e-10
-@test length(sol_sp.t)==length(sol.t)
-
-sol_sp=solve(prob_sp,Rodas5(autodiff=false))
-sol=solve(prob_std,Rodas5(autodiff=false))
-@test sol_sp.u[end]≈sol.u[end]
-@test length(sol_sp.t)==length(sol.t)
-
-sol_sp=solve(prob_sp,Rodas5())
-sol=solve(prob_std,Rodas5())
-@test sol_sp.u[end]≈sol.u[end]
-@test length(sol_sp.t)==length(sol.t)
+
+for f in [f_ip, f_oop]
+  odefun_std = ODEFunction(f)
+  prob_std = ODEProblem(odefun_std,u0,tspan)
+
+  for ad in [true, false]
+    for tol in [nothing, 1e-10]
+      sol_std=solve(prob_std,Rodas5(autodiff=ad),reltol=tol,abstol=tol)
+      @test sol_std.retcode==:Success
+      for prob in map(f->ODEProblem(f,u0,tspan),
+                      [ODEFunction(f,colorvec=colors,jac_prototype=jac_sp),
+                       ODEFunction(f,jac_prototype=jac_sp),
+                       #ODEFunction(f,colorvec=colors) # this one fails both the u[end] and length tests
+                       ])
+        sol=solve(prob,Rodas5(autodiff=ad),reltol=tol,abstol=tol)
+        @test sol.retcode==:Success
+        if tol !=nothing
+          @test sol_std.u[end]≈sol.u[end] atol=tol
+        else
+          @test sol_std.u[end]≈sol.u[end]
+        end
+        @test length(sol_std.t)==length(sol.t)
+      end
+    end
+  end
+end