Linear Combinations of Operators #55
Conversation
Now the coefficient is also taken into account. Also adds test set for `DiffEqArrayOperator`.
Codecov Report
@@ Coverage Diff @@
## master #55 +/- ##
=========================================
Coverage ? 68.98%
=========================================
Files ? 9
Lines ? 790
Branches ? 0
=========================================
Hits ? 545
Misses ? 245
Partials ? 0
Continue to review full report at Codecov.
|
MSeeker1340
changed the title
|
Seems I still haven't covered everything thoroughly... To calculate caching phi in expA = expm(A*dt)
phi1 = ((expA-I)/A)and we cannot divide by a derivative operator at the moment. So there's more work to be done. |
src/operator_combination.jl
Outdated
| Base.expm(A::LinearCombination) = expm(full(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) |
There was a problem hiding this comment.
There's a lazy way to do normest2 which will be required when A doesn't fit into memory, but we can handle that later. Probably worth a comment and an issue though.
|
I'll drop the WIP flag once things check out on the OrdinaryDiffEq end. Also, should we explicitly call expA = expm(A*dt)
phi1 = (expA - I) / Ause fullA = full(A)
expA = expm(fullA*dt)
phi1 = (expA - I) / fullAThis way we will avoid a lot of potential problems for future caching methods. |
Yes, for the cached ones we know we have to |
MSeeker1340
changed the title
|
All set. I'll do a PR on OrdinaryDiffEq after this is updated. |
src/array_operator.jl
Outdated
| end | ||
| out | ||
| end | ||
| Base.full(L::DiffEqArrayOperator) = full(L.A) * L.α.coeff |
There was a problem hiding this comment.
.* this so that way it's compatible with coefficient arrays. That'll come up in space-dependent coefficients and quasi-linear PDEs.
| 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.
IIRC isn't / defined in terms of \ with a transpose? It may be better to do that one directly. I would use @which to dive into the Base code and see how / is called.
There was a problem hiding this comment.
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:457The definition is
\(x,y) = (y'/x')'So Julia does in fact define \ in terms of / for generic types. However for regular matrices the relationship is inverted:
# In linalg\generic.jl:820
(/)(A::AbstractVecOrMat, B::AbstractVecOrMat) = (B' \ A')'Do you think I should define / using this convention, or just go with the generic one?
There was a problem hiding this comment.
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 n x n. Transpose on the dense matrix will always work though.
|
Small changes I missed last time, then this looks good to go. |
| 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')' |
There was a problem hiding this comment.
These transposes might be incorrect because they are just implemented as a no-op. Sorry if I was unclear before. I think you should just full and /, or does the generic version but full before transpose.
Adds
expmandnormmethods for linear combinations ofDiffEqArrayOperatorand derivative operators. Also includes anormboundmethod that computes a upper bound to the norm using triangle inequality.normbound(A,Inf)can serve as an efficient estimate for the operator norm in applications where the exact norm is not desired (for example in error estimation in Krylov methods), because it is easy to calculate the Inf norm for derivative operators.In addition, fixes an error in
normforDiffEqArrayOperator(it now takes into account the coefficient too) and adds a test set forDiffEqArrayOperator.The efficient
fullmethod for linear combinations is currently bugged as Julia's type dispatching will always default to thefullmethod defined in LinearMaps.jl.Not sure where to put the functions for the linear combination though, since it handles both
DiffEqArrayOperatorand derivative operators.