-
-
Notifications
You must be signed in to change notification settings - Fork 74
Linear Combinations of Operators #55
Changes from 7 commits
673bacc
aeda250
11835e8
67c81cd
8018256
dc74ec0
79897f4
ea7e198
db2838b
590acb6
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -209,3 +209,11 @@ function Base.sparse(A::AbstractDerivativeOperator{T}) where T | |
end | ||
return mat | ||
end | ||
|
||
#= | ||
Fallback methods that use the full representation of the operator | ||
=# | ||
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.:\ is also defined | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. IIRC isn't There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I commented out this portion and found that using DiffEqOperators
N = 20
A = DerivativeOperator{Float64}(2,2,1.0,N,:Dirichlet0,:Dirichlet0)
mat = rand(N,N)
@which A / mat
no method found for the specified argument types
@which mat' \ A'
\(x, y) at operators.jl:457 The definition is \(x,y) = (y'/x')' So Julia does in fact define # In linalg\generic.jl:820
(/)(A::AbstractVecOrMat, B::AbstractVecOrMat) = (B' \ A')' Do you think I should define There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I think you might need both. The transpose operator will be hard to do in general, and right now I think it's a no-op which only holds for some BCs and since we have the operator as |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,52 @@ | ||
#= | ||
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.norm(A::IdentityMap{T}, p::Real=2) where T = real(one(T)) | ||
Base.norm(A::LinearCombination, p::Real=2) = norm(full(A), p) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. There's a lazy way to do |
||
#= | ||
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)) |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
.*
this so that way it's compatible with coefficient arrays. That'll come up in space-dependent coefficients and quasi-linear PDEs.