Merge remote-tracking branch 'origin/master'

Author: ChrisRackauckas

docs/KdV.md

@@ -1,4 +1,4 @@
-# Solving the Heat Equation using DiffEqOperators
+# Solving the KdV Equation using DiffEqOperators
 
 In this tutorial we will try to solve the famous **KdV equation** which describes the motion of waves on shallow water surfaces.
 The equation is commonly written as 
@@ -29,4 +29,4 @@ Now call the ODE solver as follows:-
     sol1 = solve(prob1, dense=false, tstops=0:0.01:10);
 
 and plot the solutions to see the waveform evolve with time.
-**Note:** The waveform being solved for here is non-directional unlike the many waves you might see like traveling solitons. In that case you might need to use the Upwind operators.
+**Note:** The waveform being solved for here is non-directional unlike the many waves you might see like traveling solitons. In that case you might need to use the Upwind operators.

src/DiffEqOperators.jl

@@ -8,13 +8,17 @@ using SparseArrays
 
 abstract type AbstractDerivativeOperator{T} <: AbstractDiffEqLinearOperator{T} end
 abstract type AbstractDiffEqCompositeOperator{T} <: AbstractDiffEqLinearOperator{T} end
+abstract type AbstractMatrixFreeOperator{T} <: AbstractDiffEqLinearOperator{T} end
 
 ### Common default methods for the operators
 include("common_defaults.jl")
 
 ### Basic Operators
 include("basic_operators.jl")
 
+### Matrix-free Operators
+include("matrixfree_operators.jl")
+
 ### Derivative Operators
 include("derivative_operators/fornberg.jl")
 include("derivative_operators/upwind_operator.jl")
@@ -33,6 +37,7 @@ for T in [DiffEqScalar, DiffEqArrayOperator, FactorizedDiffEqArrayOperator, Diff
   (L::T)(du,u,p,t) = (update_coefficients!(L,u,p,t); mul!(du,L,u))
 end
 
+export MatrixFreeOperator
 export DiffEqScalar, DiffEqArrayOperator, DiffEqIdentity, getops
 export AbstractDerivativeOperator, DerivativeOperator, UpwindOperator, FiniteDifference
 end # module

src/matrixfree_operators.jl

@@ -0,0 +1,64 @@
+import LinearAlgebra: mul!
+struct MatrixFreeOperator{F,N} <: AbstractMatrixFreeOperator{F}
+  f::F
+  args::N
+  function MatrixFreeOperator(f::F, args::N) where {F,N}
+    @assert N === Nothing || (N <: Tuple && length(args) == 2) "Arguments of a "*
+    "MatrixFreeOperator must be nothing or a tuple with two elements"
+    return new{F,N}(f, args)
+  end
+end
+MatrixFreeOperator(f) = MatrixFreeOperator(f, nothing)
+
+function Base.getproperty(M::MatrixFreeOperator{F,N}, sym::Symbol) where {F,N}
+  if sym === :update_func
+    return N === Nothing ? DEFAULT_UPDATE_FUNC : M.f
+  else
+    return getfield(M, sym)
+  end
+end
+
+function (M::MatrixFreeOperator{F,N})(du, u, p, t) where {F,N}
+  if N === Nothing
+    M.f(du, u)
+  else
+    M.f(du, u, p, t)
+  end
+  du
+end
+
+function (M::MatrixFreeOperator{F,N})(u, p, t) where {F,N}
+  du = similar(u)
+  if N === nothing
+    M.f(du, u)
+  else
+    M.f(du, u, p, t)
+  end
+  du
+end
+
+@inline function mul!(y::AbstractVector, A::MatrixFreeOperator{F,N}, x::AbstractVector) where {F,N}
+  if N === Nothing
+    A.f(y, x)
+  else
+    A.f(y, x, A.args[1], A.args[2])
+  end
+  y
+end
+
+@inline function mul!(Y::AbstractMatrix, A::MatrixFreeOperator{F,N}, X::AbstractMatrix) where {F,N}
+  m,  n = size(Y)
+  k, _n = size(X)
+  if n != _n
+    throw(DimensionMismatch("the second dimension of Y, $N, does not match the second dimension of X, $_n"))
+  end
+  for i in 1:n
+    y = view(Y, :, i)
+    x = view(X, :, i)
+    mul!(y, A, x)
+  end
+  Y
+end
+@inline Base.:*(A::MatrixFreeOperator, X::AbstractVecOrMat) = mul!(similar(X), A, X)
+
+DiffEqBase.numargs(::MatrixFreeOperator) = 4

test/matrixfree.jl

@@ -0,0 +1,23 @@
+using Test, LinearAlgebra
+using DiffEqOperators, OrdinaryDiffEq
+
+@testset "Matrix Free Operator" begin
+  f = cumsum!
+  A = MatrixFreeOperator(f)
+  b = rand(5)
+  prob = ODEProblem(A, b, (0,1.))
+  @test_nowarn solve(prob, Tsit5())
+  prob = SplitODEProblem(A, (du,u,p,t)->0, b, (0,1.))
+  Base.size(::typeof(A), n) = n==1 || n == 2 ? 5 : 1
+  LinearAlgebra.opnorm(::typeof(A), n::Real) = n==Inf ? 5 : nothing
+  LinearAlgebra.ishermitian(::typeof(A)) = false
+  @test_nowarn solve(prob, LawsonEuler(krylov=true, m=5), dt=0.1)
+  f = (du, u, p, t) -> cumsum!(du, u)
+  args = (1, 1)
+  A = MatrixFreeOperator(f,args)
+  b = rand(5)
+  Base.size(::typeof(A), n) = n==1 || n == 2 ? 5 : 1
+  LinearAlgebra.opnorm(::typeof(A), n::Real) = n==Inf ? 5 : nothing
+  LinearAlgebra.ishermitian(::typeof(A)) = false
+  @test_nowarn solve(prob, LawsonEuler(krylov=true, m=5), dt=0.1)
+end

test/runtests.jl

@@ -14,6 +14,7 @@ import Base: isapprox
 @time @testset "Finite Difference Operator" begin include("generic_operator_check.jl") end
 #@time @testset "KdV" begin include("KdV.jl") end # KdV times out and all fails
 @time @testset "Heat Equation" begin include("heat_eqn.jl") end
+@time @testset "Matrix-Free Operators" begin include("matrixfree.jl") end