Testing lazy W operator for in-place ImplictEuler

testW.jl

using LinearAlgebra, Random, OrdinaryDiffEq, DiffEqOperators, Test
@testset "In-place ImplicitEuler" begin
  A = [-1.0 0.0; 0.0 -0.5]; u0 = [1.0, 1.0]; tspan = (0.0,1.0)
  srand(0); mm = [2.0 0.0; 0.0 1.0]; mmA = [-0.5 0.0; 0.0 -0.5]

  fun = ODEFunction((du,u,p,t) -> mul!(du,A,u);
                    analytic=(u0,p,t) -> exp(mmA*t)*u0)
  prob = ODEProblem(fun,u0,tspan; mass_matrix=mm)
  sol = solve(prob, ImplicitEuler(); adaptive=false, dt=0.01)
  err = norm(sol(1.0) - fun.analytic(u0,nothing,1.0))
  println("Dense W error = $err")

  fun = ODEFunction((du,u,p,t) -> mul!(du,A,u);
                    jac_prototype=DiffEqArrayOperator(A),
                    analytic=(u0,p,t) -> exp(mmA*t)*u0)
  prob = ODEProblem(fun,u0,tspan; mass_matrix=mm)

  sol = solve(prob, ImplicitEuler(); adaptive=false, dt=0.01)
  err = norm(sol(1.0) - fun.analytic(u0,nothing,1.0))
  println("Lazy W error = $err")

  # Check that WOperator handling is correct
  integrator = init(prob, ImplicitEuler(); adaptive=false, dt=0.01)
  # Check cache initialization
  @test isa(integrator.cache.J, DiffEqArrayOperator)
  @test isa(integrator.cache.W, OrdinaryDiffEq.WOperator)
  step!(integrator)
  # Check that W is indeed updated
  @test integrator.cache.W.mass_matrix == mm
end