Basic operator composition

src/DiffEqOperators.jl

@@ -2,15 +2,19 @@ __precompile__()
 
 module DiffEqOperators
 
-import Base: *, getindex
+import Base: +, -, *, /, \, size, getindex, setindex!, Matrix, convert
 using DiffEqBase, StaticArrays, LinearAlgebra
-import DiffEqBase: update_coefficients, update_coefficients!
+import LinearAlgebra: mul!, ldiv!, lmul!, rmul!, axpy!, opnorm, factorize
+import DiffEqBase: AbstractDiffEqLinearOperator, update_coefficients!, is_constant
 
-abstract type AbstractDerivativeOperator{T} <: DiffEqBase.AbstractDiffEqLinearOperator{T} end
+abstract type AbstractDerivativeOperator{T} <: AbstractDiffEqLinearOperator{T} end
+abstract type AbstractDiffEqCompositeOperator{T} <: AbstractDiffEqLinearOperator{T} end
+
+### Common default methods for the operators
+include("common_defaults.jl")
 
 ### Basic Operators
-include("diffeqscalar.jl")
-include("array_operator.jl")
+include("basic_operators.jl")
 
 ### Derivative Operators
 include("derivative_operators/fornberg.jl")
@@ -20,10 +24,16 @@ include("derivative_operators/derivative_operator.jl")
 include("derivative_operators/abstract_operator_functions.jl")
 include("derivative_operators/boundary_operators.jl")
 
-### Linear Combination of Operators
-#include("operator_combination.jl")
+### Composite Operators
+include("composite_operators.jl")
+
+# The (u,p,t) and (du,u,p,t) interface
+for T in [DiffEqScalar, DiffEqArrayOperator, FactorizedDiffEqArrayOperator, DiffEqIdentity,
+  DiffEqScaledOperator, DiffEqOperatorCombination, DiffEqOperatorComposition]
+  (L::T)(u,p,t) = (update_coefficients!(L,u,p,t); L * u)
+  (L::T)(du,u,p,t) = (update_coefficients!(L,u,p,t); mul!(du,L,u))
+end
 
-export DiffEqScalar, DiffEqArrayOperator
+export DiffEqScalar, DiffEqArrayOperator, DiffEqIdentity
 export AbstractDerivativeOperator, DerivativeOperator, UpwindOperator, FiniteDifference
-export opnormbound
 end # module

src/basic_operators.jl

@@ -0,0 +1,94 @@
+struct DiffEqIdentity{T,N} <: AbstractDiffEqLinearOperator{T} end
+DiffEqIdentity(u) = DiffEqIdentity{eltype(u),length(u)}()
+size(::DiffEqIdentity{T,N}) where {T,N} = (N,N)
+size(::DiffEqIdentity{T,N}, m::Integer) where {T,N} = (m == 1 || m == 2) ? N : 1
+opnorm(::DiffEqIdentity{T,N}, p::Real=2) where {T,N} = one(T)
+convert(::Type{AbstractMatrix}, ::DiffEqIdentity{T,N}) where {T,N} = Diagonal(ones(T,N))
+for op in (:*, :/, :\)
+  @eval $op(::DiffEqIdentity{T,N}, x::AbstractVecOrMat) where {T,N} = $op(I, x)
+  @eval $op(x::AbstractVecOrMat, ::DiffEqIdentity{T,N}) where {T,N} = $op(x, I)
+end
+mul!(Y::AbstractVecOrMat, ::DiffEqIdentity, B::AbstractVecOrMat) = Y .= B
+ldiv!(Y::AbstractVecOrMat, ::DiffEqIdentity, B::AbstractVecOrMat) = Y .= B
+for pred in (:isreal, :issymmetric, :ishermitian, :isposdef)
+  @eval LinearAlgebra.$pred(::DiffEqIdentity) = true
+end
+
+"""
+    DiffEqScalar(val[; update_func])
+
+Represents a time-dependent scalar/scaling operator. The update function
+is called by `update_coefficients!` and is assumed to have the following
+signature:
+
+    update_func(oldval,u,p,t) -> newval
+
+You can also use `setval!(α,val)` to bypass the `update_coefficients!`
+interface and directly mutate the scalar's value.
+"""
+mutable struct DiffEqScalar{T<:Number,F} <: AbstractDiffEqLinearOperator{T}
+  val::T
+  update_func::F
+  DiffEqScalar(val::T; update_func=DEFAULT_UPDATE_FUNC) where {T} =
+    new{T,typeof(update_func)}(val, update_func)
+end
+
+size(::DiffEqScalar) = ()
+size(::DiffEqScalar, ::Integer) = 1
+update_coefficients!(α::DiffEqScalar,u,p,t) = (α.val = α.update_func(α.val,u,p,t); α)
+setval!(α::DiffEqScalar, val) = (α.val = val; α)
+is_constant(α::DiffEqScalar) = α.update_func == DEFAULT_UPDATE_FUNC
+
+for op in (:*, :/, :\)
+  @eval $op(α::DiffEqScalar, x::AbstractVecOrMat) = $op(α.val, x)
+  @eval $op(x::AbstractVecOrMat, α::DiffEqScalar) = $op(x, α.val)
+end
+lmul!(α::DiffEqScalar, B::AbstractVecOrMat) = lmul!(α.val, B)
+rmul!(B::AbstractVecOrMat, α::DiffEqScalar) = rmul!(B, α.val)
+mul!(Y::AbstractVecOrMat, α::DiffEqScalar, B::AbstractVecOrMat) = mul!(Y, α.val, B)
+axpy!(α::DiffEqScalar, X::AbstractVecOrMat, Y::AbstractVecOrMat) = axpy!(α.val, X, Y)
+Base.abs(α::DiffEqScalar) = abs(α.val)
+
+"""
+    DiffEqArrayOperator(A[; update_func])
+
+Represents a time-dependent linear operator given by an AbstractMatrix. The
+update function is called by `update_coefficients!` and is assumed to have
+the following signature:
+
+    update_func(A::AbstractMatrix,u,p,t) -> [modifies A]
+
+You can also use `setval!(α,A)` to bypass the `update_coefficients!` interface
+and directly mutate the array's value.
+"""
+mutable struct DiffEqArrayOperator{T,AType<:AbstractMatrix{T},F} <: AbstractDiffEqLinearOperator{T}
+  A::AType
+  update_func::F
+  DiffEqArrayOperator(A::AType; update_func=DEFAULT_UPDATE_FUNC) where {AType} = 
+    new{eltype(A),AType,typeof(update_func)}(A, update_func)
+end
+
+update_coefficients!(L::DiffEqArrayOperator,u,p,t) = (L.update_func(L.A,u,p,t); L)
+setval!(L::DiffEqArrayOperator, A) = (L.A = A; L)
+is_constant(L::DiffEqArrayOperator) = L.update_func == DEFAULT_UPDATE_FUNC
+
+convert(::Type{AbstractMatrix}, L::DiffEqArrayOperator) = L.A
+setindex!(L::DiffEqArrayOperator, v, i::Int) = (L.A[i] = v)
+setindex!(L::DiffEqArrayOperator, v, I::Vararg{Int, N}) where {N} = (L.A[I...] = v)
+
+"""
+    FactorizedDiffEqArrayOperator(F)
+
+Like DiffEqArrayOperator, but stores a Factorization instead.
+
+Supports left division and `ldiv!` when applied to an array.
+"""
+struct FactorizedDiffEqArrayOperator{T<:Number,FType<:Factorization{T}} <: AbstractDiffEqLinearOperator{T}
+  F::FType
+end
+
+Matrix(L::FactorizedDiffEqArrayOperator) = Matrix(L.F)
+convert(::Type{AbstractMatrix}, L::FactorizedDiffEqArrayOperator) = convert(AbstractMatrix, L.F)
+size(L::FactorizedDiffEqArrayOperator, args...) = size(L.F, args...)
+ldiv!(Y::AbstractVecOrMat, L::FactorizedDiffEqArrayOperator, B::AbstractVecOrMat) = ldiv!(Y, L.F, B)
+\(L::FactorizedDiffEqArrayOperator, x::AbstractVecOrMat) = L.F \ x

src/common_defaults.jl

@@ -0,0 +1,34 @@
+# The `update_coefficients!` interface
+DEFAULT_UPDATE_FUNC(A,u,p,t) = A # no-op used by the basic operators
+# is_constant(::AbstractDiffEqLinearOperator) = true # already defined in DiffEqBase
+update_coefficients!(L::AbstractDiffEqLinearOperator,u,p,t) = L
+
+# Routines that use the AbstractMatrix representation
+convert(::Type{AbstractArray}, L::AbstractDiffEqLinearOperator) = convert(AbstractMatrix, L)
+size(L::AbstractDiffEqLinearOperator, args...) = size(convert(AbstractMatrix,L), args...)
+opnorm(L::AbstractDiffEqLinearOperator, p::Real=2) = opnorm(convert(AbstractMatrix,L), p)
+getindex(L::AbstractDiffEqLinearOperator, i::Int) = convert(AbstractMatrix,L)[i]
+getindex(L::AbstractDiffEqLinearOperator, I::Vararg{Int, N}) where {N} = 
+  convert(AbstractMatrix,L)[I...]
+for op in (:*, :/, :\)
+  @eval $op(L::AbstractDiffEqLinearOperator, x::AbstractVecOrMat) = $op(convert(AbstractMatrix,L), x)
+  @eval $op(x::AbstractVecOrMat, L::AbstractDiffEqLinearOperator) = $op(x, convert(AbstractMatrix,L))
+end
+mul!(Y::AbstractVecOrMat, L::AbstractDiffEqLinearOperator, B::AbstractVecOrMat) =
+  mul!(Y, convert(AbstractMatrix,L), B)
+ldiv!(Y::AbstractVecOrMat, L::AbstractDiffEqLinearOperator, B::AbstractVecOrMat) =
+  ldiv!(Y, convert(AbstractMatrix,L), B)
+for pred in (:isreal, :issymmetric, :ishermitian, :isposdef)
+  @eval LinearAlgebra.$pred(L::AbstractDiffEqLinearOperator) = $pred(convert(AbstractArray, L))
+end
+factorize(L::AbstractDiffEqLinearOperator) = 
+  FactorizedDiffEqArrayOperator(factorize(convert(AbstractMatrix, L)))
+for fact in (:lu, :lu!, :qr, :qr!, :chol, :chol!, :ldlt, :ldlt!, 
+  :bkfact, :bkfact!, :lq, :lq!, :svd, :svd!)
+  @eval LinearAlgebra.$fact(L::AbstractDiffEqLinearOperator, args...) = 
+    FactorizedDiffEqArrayOperator($fact(convert(AbstractMatrix, L), args...))
+end
+
+# Routines that use the full matrix representation
+Matrix(L::AbstractDiffEqLinearOperator) = Matrix(convert(AbstractMatrix, L))
+LinearAlgebra.exp(L::AbstractDiffEqLinearOperator) = exp(Matrix(L))