-
-
Notifications
You must be signed in to change notification settings - Fork 10
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Lazy concatenation / block diagonal forms of operators #161
Comments
No, we don't have that functionality yet. It is on my to-do list. Would make writing vector-calculus operators really convenient. And will be very useful for solving coupled PDE systems. |
Potential API: abstract type Abstract BlockOperator{T} <: AbstractSciMLOperator{T} end
###
# Block Operator
###
struct BlockOperator{T, Ti, Tb}
Mblocks:Ti
Nblocks::Ti
blocks::Tb
size_info
function BlockOperator(Mblocks::Integer, Mblocks::Integer, blocks)
@assert Mblocks * Nblocks == length(blocks)
size_consistent = true
for
# loop over rows, cols, to ensure size consistency.
end
@assert size_consistent
T = promote_type(eltype.(blocks)...)
BlockOperator{T, typeof(M), typeof(blocks)}(Mblocks, Nblocks, blocks)
end
end
function BlockOperator(A::AbstractVecOrMat{<:AbstractSciMLOperator})
BlockOperator(size(A)..., vec(A))
end
function Base.:*(L::BlockOperator, u::AbstractVecOrMat)
end
L = BlockOperator([A1, A2])
###
# Block Diagonal
###
struct BlockDiagonal{T, Tb}
blocks::Tb
function BlockDiagonal(blocks... :: AbstractSciMLOperator)
T = promote_type(eltype.(blocks)...)
BlockDiagonal{T, typeof(blocks)}(blocks)
end
end
function Base.:*(L::BlockOperator, u::AbstractVecOrMat)
vcat(Tuple(B * u for B in L.blocks)...)
end |
if we support sparse block patterns, then block diagonal becomes a special case |
Yeah for sparse patterns, we can have a constructor that accepts pairs |
It might make sense to implement |
Can always do |
Also want to be mindful of code dup with kronecker products. |
I remember putting in short-circuits in |
G = BlockOperator([Dx, Dy])
T = DiagonalBlock([Tx, Ty])
Lapl = G' * T * G # assume sizes match
Lapl == Dx' * Tx * Dx + Dy' * Ty * Dy # true so we need to overload Base.:*(A::AbstractBlockOperator, B::AbstractBlockOperator) |
linking old block matrix issue here for completeness. #45 |
Does there exist functionality analogous to those listed in https://julialinearalgebra.github.io/LinearMaps.jl/stable/types/#BlockMap-and-BlockDiagonalMap?
In particular, lazy hcat/vcat/hvcat/diagonal op formation, etc.
We should also revisit this line of code:
SciMLOperators.jl/src/matrix.jl
Line 116 in 31275d4
The text was updated successfully, but these errors were encountered: