Skip to content

Suppport Diagonal block times padded #66

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.

Sign up for GitHub

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

Jump to bottom
Merged
merged 8 commits into from
Feb 4, 2021
Merged

Suppport Diagonal block times padded #66

merged 8 commits into from
Feb 4, 2021

Conversation

dlfivefifty
Copy link
Member

@TSGut with this the following works:

julia> D * [[1,2,3]; zeros(∞)]
∞-element LazyArrays.CachedArray{Float64,1,Array{Float64,1},Zeros{Float64,1,Tuple{BlockedUnitRange{RangeCumsum{Int64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}} with indices 1:1:∞:
  0.0
 ────
 -4.0
 -6.0
  0.0
 ────
  0.0
  0.0
  0.0
  0.0
  0.0
 ────
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
 ────
  0.0
  0.0
  ⋮

I can't do your example of D * ones(∞) because of another bug. Do you really need this right now or was it just a test problem?

@TSGut
Copy link
Contributor

TSGut commented Feb 4, 2021 •
edited
Loading

The ones() thing was just an example. For the application I have in mind the multiplication would be between a block diagonal operator and a block vector, in the language of HarmonicOrthogonalPolynomials this is basically what Δ*S*(S\f.(xyz)) comes down to if we want to do it adaptively.

@dlfivefifty
Copy link
Member Author

Ok so that seems close to my example, so hopefully it’ll work 🤞

@TSGut
Copy link
Contributor

TSGut commented Feb 4, 2021

I can work with it as it is now with with a minor workaround, so there is no rush to get this completely ironed out.

As it is, I'm seeing what appears to be two different kinds of errors now. Your example of D * [[1,2,3]; zeros(∞)] works for me as well. Running D*v where v is a BlockVector - e.g. running D*mortar(Fill.((0:∞), 1:2:∞)) - appears to get stuck in an infinite loop somewhere, so it never finishes. The ones(∞) example throws an actual error though:

julia> D*ones(∞)
ERROR: LoadError: ArgumentError: Cannot create infinite Array
Stacktrace:
  [1] Vector{Vector{Float64}}(#unused#::UndefInitializer, #unused#::Tuple{InfiniteArrays.Infinity})
    @ InfiniteArrays ~/.julia/packages/InfiniteArrays/YDweO/src/infarrays.jl:19
  [2] _BlockArray(#unused#::Type{Vector{Vector{Float64}}}, baxes::Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.InfStepRange{Int64, Int64}}}})
    @ BlockArrays ~/.julia/packages/BlockArrays/Ssbrg/src/blockarray.jl:92
  [3] macro expansion
    @ ~/.julia/packages/BlockArrays/Ssbrg/src/blockarray.jl:126 [inlined]
  [4] initialized_blocks_BlockArray(#unused#::Type{Vector{Vector{Float64}}}, baxes::Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.InfStepRange{Int64, Int64}}}})
    @ BlockArrays ~/.julia/packages/BlockArrays/Ssbrg/src/blockarray.jl:124
  [5] BlockArray
    @ ~/.julia/packages/BlockArrays/Ssbrg/src/blockarray.jl:140 [inlined]
  [6] similar
    @ ~/.julia/packages/BlockArrays/Ssbrg/src/blockarray.jl:354 [inlined]
  [7] similar(a::BlockVector{Int64, BroadcastVector{Fill{Int64, 1, Tuple{OneTo{Int64}}}, Type{Fill}, Tuple{BroadcastVector{Int64, typeof(-), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, BroadcastVector{Int64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, InfiniteArrays.InfUnitRange{Int64}, Base.RefValue{Val{2}}}}}}, InfiniteArrays.InfStepRange{Int64, Int64}}}, Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.InfStepRange{Int64, Int64}}}}}, #unused#::Type{Float64})
    @ BlockArrays ~/.julia/packages/BlockArrays/Ssbrg/src/abstractblockarray.jl:37
  [8] AbstractVector{Float64}(A::BlockVector{Int64, BroadcastVector{Fill{Int64, 1, Tuple{OneTo{Int64}}}, Type{Fill}, Tuple{BroadcastVector{Int64, typeof(-), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, BroadcastVector{Int64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, InfiniteArrays.InfUnitRange{Int64}, Base.RefValue{Val{2}}}}}}, InfiniteArrays.InfStepRange{Int64, Int64}}}, Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.InfStepRange{Int64, Int64}}}}})
    @ Base ./array.jl:541
  [9] convert(#unused#::Type{AbstractVector{Float64}}, a::BlockVector{Int64, BroadcastVector{Fill{Int64, 1, Tuple{OneTo{Int64}}}, Type{Fill}, Tuple{BroadcastVector{Int64, typeof(-), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, BroadcastVector{Int64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, InfiniteArrays.InfUnitRange{Int64}, Base.RefValue{Val{2}}}}}}, InfiniteArrays.InfStepRange{Int64, Int64}}}, Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.InfStepRange{Int64, Int64}}}}})
    @ Base ./abstractarray.jl:16
 [10] copy_oftype(A::BlockVector{Int64, BroadcastVector{Fill{Int64, 1, Tuple{OneTo{Int64}}}, Type{Fill}, Tuple{BroadcastVector{Int64, typeof(-), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, BroadcastVector{Int64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, InfiniteArrays.InfUnitRange{Int64}, Base.RefValue{Val{2}}}}}}, InfiniteArrays.InfStepRange{Int64, Int64}}}, Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.InfStepRange{Int64, Int64}}}}}, #unused#::Type{Float64})
    @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/LinearAlgebra.jl:350
 [11] broadcasted(#unused#::BlockArrays.BlockStyle{1}, #unused#::typeof(*), a::BlockVector{Int64, BroadcastVector{Fill{Int64, 1, Tuple{OneTo{Int64}}}, Type{Fill}, Tuple{BroadcastVector{Int64, typeof(-), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, BroadcastVector{Int64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, InfiniteArrays.InfUnitRange{Int64}, Base.RefValue{Val{2}}}}}}, InfiniteArrays.InfStepRange{Int64, Int64}}}, Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.InfStepRange{Int64, Int64}}}}}, b::Ones{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}})
    @ BlockArrays ~/.julia/packages/BlockArrays/Ssbrg/src/blockbroadcast.jl:292
 [12] broadcasted
    @ ./broadcast.jl:1313 [inlined]
 [13] broadcast(::typeof(*), ::BlockVector{Int64, BroadcastVector{Fill{Int64, 1, Tuple{OneTo{Int64}}}, Type{Fill}, Tuple{BroadcastVector{Int64, typeof(-), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, BroadcastVector{Int64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, InfiniteArrays.InfUnitRange{Int64}, Base.RefValue{Val{2}}}}}}, InfiniteArrays.InfStepRange{Int64, Int64}}}, Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.InfStepRange{Int64, Int64}}}}}, ::Ones{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}})
    @ Base.Broadcast ./broadcast.jl:821
 [14] layout_broadcasted(#unused#::BlockArrays.BlockLayout{LazyArrays.BroadcastLayout{Type{Fill}}, ArrayLayouts.FillLayout}, #unused#::LazyArrays.CachedLayout{ArrayLayouts.DenseColumnMajor, ArrayLayouts.OnesLayout}, op::Function, A::BlockVector{Int64, BroadcastVector{Fill{Int64, 1, Tuple{OneTo{Int64}}}, Type{Fill}, Tuple{BroadcastVector{Int64, typeof(-), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, BroadcastVector{Int64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, InfiniteArrays.InfUnitRange{Int64}, Base.RefValue{Val{2}}}}}}, InfiniteArrays.InfStepRange{Int64, Int64}}}, Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.InfStepRange{Int64, Int64}}}}}, B::LazyArrays.CachedArray{Float64, 1, Vector{Float64}, Ones{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}})
    @ LazyArrays ~/.julia/packages/LazyArrays/6zNYO/src/cache.jl:280
 [15] layout_broadcasted(op::Function, A::BlockVector{Int64, BroadcastVector{Fill{Int64, 1, Tuple{OneTo{Int64}}}, Type{Fill}, Tuple{BroadcastVector{Int64, typeof(-), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, BroadcastVector{Int64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, InfiniteArrays.InfUnitRange{Int64}, Base.RefValue{Val{2}}}}}}, InfiniteArrays.InfStepRange{Int64, Int64}}}, Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.InfStepRange{Int64, Int64}}}}}, B::LazyArrays.CachedArray{Float64, 1, Vector{Float64}, Ones{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}})
    @ LazyArrays ~/.julia/packages/LazyArrays/6zNYO/src/cache.jl:297
 [16] broadcasted(#unused#::LazyArrays.LazyArrayStyle{1}, op::Function, A::BlockVector{Int64, BroadcastVector{Fill{Int64, 1, Tuple{OneTo{Int64}}}, Type{Fill}, Tuple{BroadcastVector{Int64, typeof(-), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, BroadcastVector{Int64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, InfiniteArrays.InfUnitRange{Int64}, Base.RefValue{Val{2}}}}}}, InfiniteArrays.InfStepRange{Int64, Int64}}}, Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.InfStepRange{Int64, Int64}}}}}, B::LazyArrays.CachedArray{Float64, 1, Vector{Float64}, Ones{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}})
    @ LazyArrays ~/.julia/packages/LazyArrays/6zNYO/src/cache.jl:301
 [17] broadcasted
    @ ./broadcast.jl:1313 [inlined]
 [18] copy(M::ArrayLayouts.Lmul{ArrayLayouts.DiagonalLayout{BlockArrays.BlockLayout{LazyArrays.BroadcastLayout{Type{Fill}}, ArrayLayouts.FillLayout}}, LazyArrays.CachedLayout{ArrayLayouts.DenseColumnMajor, ArrayLayouts.OnesLayout}, Diagonal{Int64, BlockVector{Int64, BroadcastVector{Fill{Int64, 1, Tuple{OneTo{Int64}}}, Type{Fill}, Tuple{BroadcastVector{Int64, typeof(-), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, BroadcastVector{Int64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, InfiniteArrays.InfUnitRange{Int64}, Base.RefValue{Val{2}}}}}}, InfiniteArrays.InfStepRange{Int64, Int64}}}, Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.InfStepRange{Int64, Int64}}}}}}, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, Ones{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/ATRic/src/diagonal.jl:20
 [19] copy(M::Mul{ArrayLayouts.DiagonalLayout{BlockArrays.BlockLayout{LazyArrays.BroadcastLayout{Type{Fill}}, ArrayLayouts.FillLayout}}, LazyArrays.CachedLayout{ArrayLayouts.DenseColumnMajor, ArrayLayouts.OnesLayout}, Diagonal{Int64, BlockVector{Int64, BroadcastVector{Fill{Int64, 1, Tuple{OneTo{Int64}}}, Type{Fill}, Tuple{BroadcastVector{Int64, typeof(-), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, BroadcastVector{Int64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, InfiniteArrays.InfUnitRange{Int64}, Base.RefValue{Val{2}}}}}}, InfiniteArrays.InfStepRange{Int64, Int64}}}, Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.InfStepRange{Int64, Int64}}}}}}, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, Ones{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/ATRic/src/mul.jl:108
 [20] materialize(M::Mul{ArrayLayouts.DiagonalLayout{BlockArrays.BlockLayout{LazyArrays.BroadcastLayout{Type{Fill}}, ArrayLayouts.FillLayout}}, LazyArrays.CachedLayout{ArrayLayouts.DenseColumnMajor, ArrayLayouts.OnesLayout}, Diagonal{Int64, BlockVector{Int64, BroadcastVector{Fill{Int64, 1, Tuple{OneTo{Int64}}}, Type{Fill}, Tuple{BroadcastVector{Int64, typeof(-), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, BroadcastVector{Int64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, InfiniteArrays.InfUnitRange{Int64}, Base.RefValue{Val{2}}}}}}, InfiniteArrays.InfStepRange{Int64, Int64}}}, Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.InfStepRange{Int64, Int64}}}}}}, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, Ones{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/ATRic/src/mul.jl:105
 [21] mul
    @ ~/.julia/packages/ArrayLayouts/ATRic/src/mul.jl:106 [inlined]
 [22] *(A::Diagonal{Int64, BlockVector{Int64, BroadcastVector{Fill{Int64, 1, Tuple{OneTo{Int64}}}, Type{Fill}, Tuple{BroadcastVector{Int64, typeof(-), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, BroadcastVector{Int64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, InfiniteArrays.InfUnitRange{Int64}, Base.RefValue{Val{2}}}}}}, InfiniteArrays.InfStepRange{Int64, Int64}}}, Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.InfStepRange{Int64, Int64}}}}}}, B::LazyArrays.CachedArray{Float64, 1, Vector{Float64}, Ones{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/ATRic/src/mul.jl:135
 [23] top-level scope

@codecov
Copy link

codecov bot commented Feb 4, 2021 •
edited
Loading

Codecov Report

Merging #66 (51e8fb8) into master (b90447f) will increase coverage by 0.31%.
The diff coverage is 100.00%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master      #66      +/-   ##
==========================================
+ Coverage   66.38%   66.70%   +0.31%     
==========================================
  Files           9        9              
  Lines         940      949       +9     
==========================================
+ Hits          624      633       +9     
  Misses        316      316
Impacted Files Coverage Δ
src/blockbanded/blockbanded.jl 60.60% <100.00%> (+14.77%) ⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 7a730b5...51e8fb8. Read the comment docs.

@dlfivefifty
Copy link
Member Author

@TSGut this is done, for now. There is an issue with

x̃ = PseudoBlockArray(x, (axes(D,1),))
D * x̃

that it is left lazy, whereas D*x hits a special case of broadcasted(::LazyArrayStyle, op, A::AbstractVector, B::CachedVector) that applies it to the cache.

There's an easy work around by overloading broadcasted(::LazyArrayStyle, op, A::AbstractVector, B::PseudoBlockVector) to apply the operator to the blocks of B. But I'm not sure this is optimal in the application, so I'll leave it lazy for the time being.

@dlfivefifty dlfivefifty merged commit 288d44e into master Feb 4, 2021
@dlfivefifty dlfivefifty deleted the dl/blockdiagmultiply branch February 4, 2021 22:55
@TSGut TSGut mentioned this pull request Feb 5, 2021
Merged
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers
No reviews
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.
2 participants
@dlfivefifty @TSGut