Linear Combinations of Operators

src/DiffEqOperators.jl

@@ -3,6 +3,7 @@ __precompile__()
 module DiffEqOperators
 
 import LinearMaps: LinearMap, AbstractLinearMap
+using LinearMaps: LinearCombination, IdentityMap
 import Base: *, getindex
 using DiffEqBase, StaticArrays
 import DiffEqBase: update_coefficients, update_coefficients!
@@ -21,6 +22,10 @@ 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")
+
 export DiffEqScalar, DiffEqArrayOperator
 export AbstractDerivativeOperator, DerivativeOperator, UpwindOperator, FiniteDifference
+export normbound
 end # module

src/array_operator.jl

@@ -52,18 +52,11 @@ Base.issymmetric(L::DiffEqArrayOperator) = L._issymmetric
 Base.ishermitian(L::DiffEqArrayOperator) = L._ishermitian
 Base.isposdef(L::DiffEqArrayOperator) = L._isposdef
 DiffEqBase.is_constant(L::DiffEqArrayOperator) = L.update_func == DEFAULT_UPDATE_FUNC
-function Base.expm(L::DiffEqArrayOperator)
-    tmp = full(L.A) # If not lazy then this is a no-op
-    tmp .*= L.α.coeff
-    out = expm(tmp)
-    if tmp === L.A
-        L.A ./= L.α.coeff # Undo change if not lazy
-    end
-    out
-end
+Base.full(L::DiffEqArrayOperator) = full(L.A) .* L.α.coeff
+Base.expm(L::DiffEqArrayOperator) = expm(full(L))
 DiffEqBase.has_expm(L::DiffEqArrayOperator) = true
 Base.size(L::DiffEqArrayOperator) = size(L.A)
-Base.norm(L::DiffEqArrayOperator, p::Real=2) = norm(L.A, p)
+Base.norm(L::DiffEqArrayOperator, p::Real=2) = norm(L.A, p) * abs(L.α.coeff)
 DiffEqBase.update_coefficients!(L::DiffEqArrayOperator,t,u) = (L.update_func(L.A,t,u); L.α = L.α(t); nothing)
 DiffEqBase.update_coefficients(L::DiffEqArrayOperator,t,u)  = (L.update_func(L.A,t,u); L.α = L.α(t); L)
 

src/derivative_operators/abstract_operator_functions.jl

@@ -209,3 +209,16 @@ function Base.sparse(A::AbstractDerivativeOperator{T}) where T
     end
     return mat
 end
+
+#=
+    Fallback methods that use the full representation of the operator
+
+    As with the convention for regular matrices, right division is defined in 
+    terms of left division. (The result may not be correct for some BC as the 
+    transpose of a derivative operator is currently defined to be a no-op)
+=#
+Base.expm(A::AbstractDerivativeOperator{T}) where T = expm(full(A))
+Base.:\(A::AbstractVecOrMat, B::AbstractDerivativeOperator) = A \ full(B)
+Base.:\(A::AbstractDerivativeOperator, B::AbstractVecOrMat) = full(A) \ B
+Base.:/(A::AbstractVecOrMat, B::AbstractDerivativeOperator) = (B' \ A')'
+Base.:/(A::AbstractDerivativeOperator, B::AbstractVecOrMat) = (B' \ A')'

src/operator_combination.jl

@@ -0,0 +1,54 @@
+#=
+    Most of the functionality for linear combination of operators is already
+    covered in LinearMaps.jl.
+=#
+
+(L::LinearCombination)(u,p,t) = L*u
+(L::LinearCombination)(du,u,p,t) = A_mul_B!(du,L,u)
+
+#=
+    The fallback implementation in LinearMaps.jl effectively computes A*eye(N), 
+    which is very inefficient.
+
+    Instead, build up the full matrix for each operator iteratively.
+=#
+# TODO: Type dispatch for this is incorrect at the moment
+# function Base.full(A::LinearCombination{T,Tuple{Vararg{O}},Ts}) where {T,O<:Union{AbstractDiffEqLinearOperator,IdentityMap},Ts}
+#     out = zeros(T,size(A))
+#     for i = 1:length(A.maps)
+#         c = A.coeffs[i]
+#         op = A.maps[i]
+#         if isa(op, IdentityMap)
+#             @. out += c * eye(size(A,1))
+#         else 
+#             @. out += c * full(op)
+#         end
+#     end
+#     return out
+# end
+
+#= 
+    Fallback methods that use the full representation
+=#
+Base.expm(A::LinearCombination) = expm(full(A))
+Base.:\(A::AbstractVecOrMat, B::LinearCombination) = A \ full(B)
+Base.:\(A::LinearCombination, B::AbstractVecOrMat) = full(A) \ B
+Base.:/(A::AbstractVecOrMat, B::LinearCombination) = (B' \ A')'
+Base.:/(A::LinearCombination, B::AbstractVecOrMat) = (B' \ A')'
+
+Base.norm(A::IdentityMap{T}, p::Real=2) where T = real(one(T))
+Base.norm(A::LinearCombination, p::Real=2) = norm(full(A), p)
+#=
+    The norm of A+B is difficult to calculate, but in many applications we only 
+    need an estimate of the norm (e.g. for error analysis) so it makes sense to 
+    compute the upper bound given by the triangle inequality
+
+        |A + B| <= |A| + |B|
+
+    For derivative operators A and B, their Inf norm can be calculated easily 
+    and thus so is the Inf norm bound of A + B.
+=#
+normbound(a::Number, p::Real=2) = abs(a)
+normbound(A::AbstractArray, p::Real=2) = norm(A, p)
+normbound(A::Union{AbstractDiffEqLinearOperator,IdentityMap}, p::Real=2) = norm(A, p)
+normbound(A::LinearCombination, p::Real=2) = sum(abs.(A.coeffs) .* normbound.(A.maps, p))

test/array_operators_interface.jl

@@ -0,0 +1,20 @@
+using DiffEqOperators
+using Base.Test
+
+N = 5
+srand(0); A = rand(N,N); u = rand(N)
+L = DiffEqArrayOperator(A)
+a = 3.5
+La = L * a
+
+@test La * u ≈ (a*A) * u
+@test lufact(La) \ u ≈ (a*A) \ u
+@test norm(La) ≈ norm(a*A)
+@test expm(La) ≈ expm(a*A)
+@test La[2,3] ≈ A[2,3] # should this be La[2,3] == a*A[2,3]?
+
+update_func = (_A,t,u) -> _A .= t * A
+t = 3.0
+Atmp = zeros(N,N)
+Lt = DiffEqArrayOperator(Atmp, a, update_func) 
+@test Lt(u,nothing,t) ≈ (a*t*A) * u

test/derivative_operators_interface.jl

@@ -87,3 +87,25 @@ end
     @test G*B ≈ 8*ones(N,M) atol=1e-2
     @test A*G*B ≈ zeros(N,M) atol=1e-2
 end
+
+@testset "Linear combinations of operators" begin
+    # Only tests the additional functionality defined in "operator_combination.jl"
+    N = 10
+    srand(0); LA = DiffEqArrayOperator(rand(N,N))
+    LD = DerivativeOperator{Float64}(2,2,1.0,N,:Dirichlet0,:Dirichlet0)
+    L = 1.1*LA - 2.2*LD + 3.3*I
+    # Builds full(L) the brute-force way
+    fullL = zeros(N,N)
+    v = zeros(N)
+    for i = 1:N
+        v[i] = 1.0
+        fullL[:,i] = L*v
+        v[i] = 0.0
+    end
+    @test full(L) ≈ fullL
+    @test expm(L) ≈ expm(fullL)
+    for p in [1,2,Inf]
+        @test norm(L,p) ≈ norm(fullL,p)
+        @test normbound(L,p) ≈ 1.1*norm(LA,p) + 2.2*norm(LD,p) + 3.3
+    end
+end

test/runtests.jl

@@ -3,6 +3,7 @@ using Base.Test
 using SpecialMatrices, SpecialFunctions
 
 tic()
+@time @testset "Array Operators Interface" begin include("array_operators_interface.jl") end
 @time @testset "Derivative Operators Interface" begin include("derivative_operators_interface.jl") end
 @time @testset "Dirichlet BCs" begin include("dirichlet.jl") end
 @time @testset "Periodic BCs" begin include("periodic.jl") end