fix inf-broadcasting with blocked fills

Project.toml

@@ -1,6 +1,6 @@
 name = "InfiniteLinearAlgebra"
 uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
-version = "0.4.4"
+version = "0.4.5"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -16,22 +16,23 @@ MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 
 [compat]
-ArrayLayouts = "0.4.10"
-BandedMatrices = "0.15.23"
-BlockArrays = "0.12.12"
-BlockBandedMatrices = "0.9"
-FillArrays = "0.10"
-InfiniteArrays = "0.8.3"
-LazyArrays = "0.19"
-LazyBandedMatrices = "0.3.2"
-MatrixFactorizations = "0.6, 0.7, 0.8"
-SemiseparableMatrices = "0.2"
+ArrayLayouts = "0.5.1"
+BandedMatrices = "0.16"
+BlockArrays = "0.14"
+BlockBandedMatrices = "0.10"
+FillArrays = "0.11"
+InfiniteArrays = "0.9"
+LazyArrays = "0.20"
+LazyBandedMatrices = "0.4"
+MatrixFactorizations = "0.7.1, 0.8"
+SemiseparableMatrices = "0.2.2"
 julia = "1.5"
 
 [extras]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Random", "SpecialFunctions"]
+test = ["Test", "Random", "SpecialFunctions", "StaticArrays"]

src/InfiniteLinearAlgebra.jl

@@ -9,11 +9,11 @@ import Base.Broadcast: BroadcastStyle, Broadcasted, broadcasted
 
 import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, TridiagonalLayout, 
                         sublayout, _qr, __qr, MatLmulVec, MatLmulMat, AbstractQLayout, materialize!, diagonaldata, subdiagonaldata, supdiagonaldata,
-                        _bidiag_forwardsub!, mulreduce
+                        _bidiag_forwardsub!, mulreduce, RangeCumsum
 import BandedMatrices: BandedMatrix, _BandedMatrix, AbstractBandedMatrix, bandeddata, bandwidths, BandedColumns, bandedcolumns,
                         _default_banded_broadcast, banded_similar
 import FillArrays: AbstractFill, getindex_value, axes_print_matrix_row
-import InfiniteArrays: OneToInf, InfUnitRange, Infinity, InfStepRange, AbstractInfUnitRange, InfAxes
+import InfiniteArrays: OneToInf, InfUnitRange, Infinity, InfStepRange, AbstractInfUnitRange, InfAxes, InfRanges
 import LinearAlgebra: lmul!,  ldiv!, matprod, qr, AbstractTriangular, AbstractQ, adjoint, transpose
 import LazyArrays: applybroadcaststyle, CachedArray, CachedMatrix, CachedVector, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes, ApplyStyle, LazyArrayApplyStyle, LazyArrayStyle,
                     resizedata!, MemoryLayout,
@@ -23,11 +23,11 @@ import LazyArrays: applybroadcaststyle, CachedArray, CachedMatrix, CachedVector,
 import MatrixFactorizations: ul, ul!, _ul, ql, ql!, _ql, QLPackedQ, getL, getR, getU, reflector!, reflectorApply!, QL, QR, QRPackedQ,
                             QRPackedQLayout, AdjQRPackedQLayout, QLPackedQLayout, AdjQLPackedQLayout, LayoutQ
 
-import BlockArrays: AbstractBlockVecOrMat, sizes_from_blocks, _length, BlockedUnitRange, blockcolsupport, BlockLayout, AbstractBlockLayout
+import BlockArrays: AbstractBlockVecOrMat, sizes_from_blocks, _length, BlockedUnitRange, blockcolsupport, BlockLayout, AbstractBlockLayout, BlockSlice
 
 import BandedMatrices: BandedMatrix, bandwidths, AbstractBandedLayout, _banded_qr!, _banded_qr, _BandedMatrix
 
-import LazyBandedMatrices: ApplyBandedLayout, BroadcastBandedLayout, _krontrav_axes
+import LazyBandedMatrices: ApplyBandedLayout, BroadcastBandedLayout, _krontrav_axes, _block_interlace_axes
 
 import BlockBandedMatrices: _BlockSkylineMatrix, _BandedMatrix, _BlockSkylineMatrix, blockstart, blockstride,
         BlockSkylineSizes, BlockSkylineMatrix, BlockBandedMatrix, _BlockBandedMatrix, BlockTridiagonal,
@@ -38,7 +38,7 @@ import SemiseparableMatrices: AbstractAlmostBandedLayout, _almostbanded_qr!
 
 # BroadcastStyle(::Type{<:BandedMatrix{<:Any,<:Any,<:OneToInf}}) = LazyArrayStyle{2}()
 
-function ^(A::BandedMatrix{T,<:Any,<:OneToInf}, p::Integer) where T
+function ArrayLayouts._power_by_squaring(_, ::NTuple{2,Infinity}, A::AbstractMatrix{T}, p::Integer) where T
     if p < 0
         inv(A)^(-p)
     elseif p == 0
@@ -61,38 +61,5 @@ include("infql.jl")
 include("infqr.jl")
 include("inful.jl")
 
-#######
-# block broadcasted
-######
-
-const CumsumOneToInf2 = Cumsum{<:Any,1,<:OneToInf}
-BlockArrays.sortedunion(a::CumsumOneToInf2, ::CumsumOneToInf2) = a
-
-function BlockArrays.sortedunion(a::Vcat{Int,1,<:Tuple{<:AbstractVector{Int},InfStepRange{Int,Int}}},
-                                 b::Vcat{Int,1,<:Tuple{<:AbstractVector{Int},InfStepRange{Int,Int}}})
-    @assert a == b
-    a
-end
-
-
-map(::typeof(length), A::BroadcastArray{OneTo{Int},1,Type{OneTo}}) = A.args[1]
-map(::typeof(length), A::BroadcastArray{<:Fill,1,Type{Fill}}) = A.args[2]
-map(::typeof(length), A::BroadcastArray{<:Zeros,1,Type{Zeros}}) = A.args[1]
-map(::typeof(length), A::BroadcastArray{<:Vcat,1,Type{Vcat}}) = broadcast(+,map.(length,A.args)...)
-broadcasted(::LazyArrayStyle{1}, ::typeof(length), A::BroadcastArray{OneTo{Int},1,Type{OneTo}}) =
-    A.args[1]
-broadcasted(::LazyArrayStyle{1}, ::typeof(length), A::BroadcastArray{<:Fill,1,Type{Fill}}) =
-    A.args[2]
-
-BlockArrays._length(::BlockedUnitRange, ::OneToInf) = ∞
-BlockArrays._last(::BlockedUnitRange, ::OneToInf) = ∞
-
-###
-# KronTrav
-###
-
-_krontrav_axes(A::NTuple{N,OneToInf{Int}}, B::NTuple{N,OneToInf{Int}}) where N =
-     @. blockedrange(OneTo(length(A)))
-
 
 end # module

src/banded/infbanded.jl

@@ -373,14 +373,19 @@ _BandedMatrix(::PertToeplitzLayout, A::AbstractMatrix) =
 # UniformScaling
 ##
 
-for op in (:+, :-), Typ in (:(BandedMatrix{<:Any,<:Any,OneToInf{Int}}), 
-                            :(Adjoint{<:Any,<:BandedMatrix{<:Any,<:Any,OneToInf{Int}}}),
-                            :(Transpose{<:Any,<:BandedMatrix{<:Any,<:Any,OneToInf{Int}}}))
-    @eval begin
-        $op(A::$Typ, λ::UniformScaling) = $op(A, Diagonal(Fill(λ.λ,∞)))
-        $op(λ::UniformScaling, A::$Typ) = $op(Diagonal(Fill(λ.λ,∞)), A)
-    end
-end
+# for op in (:+, :-), Typ in (:(BandedMatrix{<:Any,<:Any,OneToInf{Int}}), 
+#                             :(Adjoint{<:Any,<:BandedMatrix{<:Any,<:Any,OneToInf{Int}}}),
+#                             :(Transpose{<:Any,<:BandedMatrix{<:Any,<:Any,OneToInf{Int}}}))
+#     @eval begin
+#         $op(A::$Typ, λ::UniformScaling) = $op(A, Diagonal(Fill(λ.λ,∞)))
+#         $op(λ::UniformScaling, A::$Typ) = $op(Diagonal(Fill(λ.λ,∞)), A)
+#     end
+# end
+
+
+
+ArrayLayouts._apply(_, ::NTuple{2,Infinity}, op, Λ::UniformScaling, A::AbstractMatrix) = op(Diagonal(Fill(Λ.λ,∞)), A)
+ArrayLayouts._apply(_, ::NTuple{2,Infinity}, op, A::AbstractMatrix, Λ::UniformScaling) = op(A, Diagonal(Fill(Λ.λ,∞)))
 
 _default_banded_broadcast(bc::Broadcasted, ::Tuple{<:OneToInf,<:Any}) = copy(Broadcasted{LazyArrayStyle{2}}(bc.f, bc.args))
 

src/blockbanded/blockbanded.jl

@@ -1,15 +1,67 @@
+const OneToInfCumsum = InfiniteArrays.RangeCumsum{Int,OneToInf{Int}}
+const OneToCumsum = InfiniteArrays.RangeCumsum{Int,OneTo{Int}}
+
+BlockArrays.sortedunion(a::OneToInfCumsum, ::OneToInfCumsum) = a
+BlockArrays.sortedunion(a::OneToCumsum, ::OneToCumsum) = a
+
+function BlockArrays.sortedunion(a::Vcat{Int,1,<:Tuple{<:AbstractVector{Int},InfStepRange{Int,Int}}},
+                                 b::Vcat{Int,1,<:Tuple{<:AbstractVector{Int},InfStepRange{Int,Int}}})
+    @assert a == b
+    a
+end
+
 sizes_from_blocks(A::AbstractVector, ::Tuple{OneToInf{Int}}) = (map(length,A),)
 
+const OneToInfBlocks = BlockedUnitRange{OneToInfCumsum}
+const OneToBlocks = BlockedUnitRange{OneToCumsum}
+
+axes(a::OneToInfBlocks) = (a,)
+axes(a::OneToBlocks) = (a,)
+
+
+function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{Ones{T,1,Tuple{OneToInfBlocks}},AbstractArray{V,N}}}) where {N,T,V}
+    a,b = bc.args
+    @assert bc.axes == axes(b)
+    convert(AbstractArray{promote_type(T,V),N}, b)
+end
+
+function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractArray{T,N},Ones{V,1,Tuple{OneToInfBlocks}}}}) where {N,T,V}
+    a,b = bc.args
+    @assert bc.axes == axes(a)
+    convert(AbstractArray{promote_type(T,V),N}, a)
+end
 
-# for LazyLay in (:(BlockLayout{LazyLayout}), :(TriangularLayout{UPLO,UNIT,BlockLayout{LazyLayout}} where {UPLO,UNIT}))
-#     @eval begin
-#         combine_mul_styles(::$LazyLay) = LazyArrayApplyStyle()
-#         mulapplystyle(::AbstractQLayout, ::$LazyLay) = LazyArrayApplyStyle()
-#     end
-# end
+_block_interlace_axes(::Int, ax::Tuple{BlockedUnitRange{OneToInf{Int}}}...) = (blockedrange(Fill(length(ax), ∞)),)
 
-# BlockArrays.blockbroadcaststyle(::LazyArrayStyle{N}) where N = LazyArrayStyle{N}()
+_block_interlace_axes(nbc::Int, ax::NTuple{2,BlockedUnitRange{OneToInf{Int}}}...) =
+    (blockedrange(Fill(length(ax) ÷ nbc, ∞)),blockedrange(Fill(mod1(length(ax),nbc), ∞)))
 
 
 include("infblocktridiagonal.jl")
 
+
+#######
+# block broadcasted
+######
+
+
+BroadcastStyle(::Type{<:SubArray{T,N,Arr,<:NTuple{N,BlockSlice{BlockRange{1,Tuple{II}}}},false}}) where {T,N,Arr<:BlockArray,II<:InfRanges} =
+    LazyArrayStyle{N}()
+
+# TODO: generalise following
+for Ax in (:(RangeCumsum{Int,OneToInf{Int}}), :(OneToInf{Int}))
+    @eval begin
+        BroadcastStyle(::Type{BlockArray{T,N,Arr,NTuple{N,BlockedUnitRange{$Ax}}}}) where {T,N,Arr} = LazyArrayStyle{N}()
+        BroadcastStyle(::Type{PseudoBlockArray{T,N,Arr,NTuple{N,BlockedUnitRange{$Ax}}}}) where {T,N,Arr} = LazyArrayStyle{N}()
+    end
+end
+
+BlockArrays._length(::BlockedUnitRange, ::OneToInf) = ∞
+BlockArrays._last(::BlockedUnitRange, ::OneToInf) = ∞
+
+###
+# KronTrav
+###
+
+_krontrav_axes(A::NTuple{N,OneToInf{Int}}, B::NTuple{N,OneToInf{Int}}) where N =
+     @. blockedrange(OneTo(length(A)))