Basic operator composition

src/DiffEqOperators.jl

@@ -10,11 +10,11 @@ import DiffEqBase: AbstractDiffEqLinearOperator, update_coefficients!, is_consta
 abstract type AbstractDerivativeOperator{T} <: AbstractDiffEqLinearOperator{T} end
 abstract type AbstractDiffEqCompositeOperator{T} <: AbstractDiffEqLinearOperator{T} end
 
-DEFAULT_UPDATE_FUNC(A,u,p,t) = A
+### 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")
@@ -25,9 +25,15 @@ include("derivative_operators/abstract_operator_functions.jl")
 include("derivative_operators/boundary_operators.jl")
 
 ### Composite Operators
-include("composite_operator.jl")
+include("composite_operators.jl")
 
-export DiffEqScalar, DiffEqArrayOperator
+# 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, 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))

src/composite_operator.jl → src/composite_operators.jl

@@ -2,43 +2,14 @@
 # derivative) using arithmetic or other operator compositions. The composite
 # operator types are lazy and maintains the structure used to build them.
 
-# Common defaults
-## Recursive routines that use `getops`
+# Recursive routines that use `getops`
 function update_coefficients!(L::AbstractDiffEqCompositeOperator,u,p,t)
   for op in getops(L)
     update_coefficients!(op,u,p,t)
   end
   L
 end
 is_constant(L::AbstractDiffEqCompositeOperator) = all(is_constant, getops(L))
-## Routines that use the AbstractMatrix representation
-size(L::AbstractDiffEqCompositeOperator, args...) = size(convert(AbstractMatrix,L), args...)
-opnorm(L::AbstractDiffEqCompositeOperator, p::Real=2) = opnorm(convert(AbstractMatrix,L), p)
-getindex(L::AbstractDiffEqCompositeOperator, i::Int) = convert(AbstractMatrix,L)[i]
-getindex(L::AbstractDiffEqCompositeOperator, I::Vararg{Int, N}) where {N} = 
-  convert(AbstractMatrix,L)[I...]
-for op in (:*, :/, :\)
-  @eval $op(L::AbstractDiffEqCompositeOperator{T}, x::AbstractVecOrMat{T}) where {T} =
-    $op(convert(AbstractMatrix,L), x)
-  @eval $op(x::AbstractVecOrMat{T}, L::AbstractDiffEqCompositeOperator{T}) where {T} =
-    $op(x, convert(AbstractMatrix,L))
-end
-mul!(Y::AbstractVecOrMat{T}, L::AbstractDiffEqCompositeOperator{T},
-  B::AbstractVecOrMat{T}) where {T} = mul!(Y, convert(AbstractMatrix,L), B)
-ldiv!(Y::AbstractVecOrMat{T}, L::AbstractDiffEqCompositeOperator{T},
-  B::AbstractVecOrMat{T}) where {T} = ldiv!(Y, convert(AbstractMatrix,L), B)
-for pred in (:isreal, :issymmetric, :ishermitian, :isposdef)
-  @eval LinearAlgebra.$pred(L::AbstractDiffEqCompositeOperator) = $pred(convert(AbstractArray, L))
-end
-factorize(L::AbstractDiffEqCompositeOperator) = 
-  FactorizedDiffEqArrayOperator(factorize(convert(AbstractArray, L)))
-for fact in (:lu, :lu!, :qr, :qr!, :chol, :chol!, :ldlt, :ldlt!, 
-  :bkfact, :bkfact!, :lq, :lq!, :svd, :svd!)
-  @eval LinearAlgebra.$fact(L::AbstractDiffEqCompositeOperator, args...) = 
-    FactorizedDiffEqArrayOperator($fact(convert(AbstractArray, L), args...))
-end
-## Routines that use the full matrix representation
-LinearAlgebra.exp(L::AbstractDiffEqCompositeOperator) = exp(Matrix(L))
 
 # Scaled operator (α * A)
 struct DiffEqScaledOperator{T,F,OpType<:AbstractDiffEqLinearOperator{T}} <: AbstractDiffEqCompositeOperator{T}
@@ -56,16 +27,16 @@ opnorm(L::DiffEqScaledOperator, p::Real=2) = abs(L.coeff) * opnorm(L.op, p)
 getindex(L::DiffEqScaledOperator, i::Int) = L.coeff * L.op[i]
 getindex(L::DiffEqScaledOperator, I::Vararg{Int, N}) where {N} = 
   L.coeff * L.op[I...]
-*(L::DiffEqScaledOperator{T,F}, x::AbstractVecOrMat{T}) where {T,F} = L.coeff * (L.op * x)
-*(x::AbstractVecOrMat{T}, L::DiffEqScaledOperator{T,F}) where {T,F} = (L.op * x) * L.coeff
-/(L::DiffEqScaledOperator{T,F}, x::AbstractVecOrMat{T}) where {T,F} = L.coeff * (L.op / x)
-/(x::AbstractVecOrMat{T}, L::DiffEqScaledOperator{T,F}) where {T,F} = 1/L.coeff * (x / L.op)
-\(L::DiffEqScaledOperator{T,F}, x::AbstractVecOrMat{T}) where {T,F} = 1/L.coeff * (L.op \ x)
-\(x::AbstractVecOrMat{T}, L::DiffEqScaledOperator{T,F}) where {T,F} = L.coeff * (x \ L)
-mul!(Y::AbstractVecOrMat{T}, L::DiffEqScaledOperator{T,F},
-  B::AbstractVecOrMat{T}) where {T,F} = lmul!(L.coeff, mul!(Y, L.op, B))
-ldiv!(Y::AbstractVecOrMat{T}, L::DiffEqScaledOperator{T,F},
-  B::AbstractVecOrMat{T}) where {T,F} = lmul!(1/L.coeff, ldiv!(Y, L.op, B))
+*(L::DiffEqScaledOperator, x::AbstractVecOrMat) = L.coeff * (L.op * x)
+*(x::AbstractVecOrMat, L::DiffEqScaledOperator) = (L.op * x) * L.coeff
+/(L::DiffEqScaledOperator, x::AbstractVecOrMat) = L.coeff * (L.op / x)
+/(x::AbstractVecOrMat, L::DiffEqScaledOperator) = 1/L.coeff * (x / L.op)
+\(L::DiffEqScaledOperator, x::AbstractVecOrMat) = 1/L.coeff * (L.op \ x)
+\(x::AbstractVecOrMat, L::DiffEqScaledOperator) = L.coeff * (x \ L)
+mul!(Y::AbstractVecOrMat, L::DiffEqScaledOperator, B::AbstractVecOrMat) =
+  lmul!(L.coeff, mul!(Y, L.op, B))
+ldiv!(Y::AbstractVecOrMat, L::DiffEqScaledOperator, B::AbstractVecOrMat) =
+  lmul!(1/L.coeff, ldiv!(Y, L.op, B))
 factorize(L::DiffEqScaledOperator) = L.coeff * factorize(L.op)
 for fact in (:lu, :lu!, :qr, :qr!, :chol, :chol!, :ldlt, :ldlt!, 
   :bkfact, :bkfact!, :lq, :lq!, :svd, :svd!)
@@ -174,9 +145,3 @@ for fact in (:lu, :lu!, :qr, :qr!, :chol, :chol!, :ldlt, :ldlt!,
   @eval LinearAlgebra.$fact(L::DiffEqOperatorComposition, args...) = 
     prod(op -> $fact(op, args...), reverse(L.ops))
 end
-
-# The (u,p,t) and (du,u,p,t) interface
-for T in [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