Dl/newblockarrays (#12)

* update for new block arrays

* tests pass!

* Block memory layouts

* CumsumBlockRange -> BlockedUnitRange

* some fixes for tests

* Update InfiniteLinearAlgebra.jl

* sortedunion for block arrays

* LazyBandedMatrices v0.2

* Update Project.toml

* Update .travis.yml

* updates for new LazyArrays

* Increase coverage

Author: dlfivefifty

.gitignore

@@ -3,3 +3,4 @@
 *.jl.mem
 deps/deps.jl
 .DS_Store
+Manifest.toml

.travis.yml

@@ -5,7 +5,6 @@ os:
   - osx
   - windows
 julia:
-  - 1.2
   - 1.3
   - nightly
 matrix:

Project.toml

@@ -1,6 +1,6 @@
 name = "InfiniteLinearAlgebra"
 uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
-version = "0.1.2"
+version = "0.2.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -17,14 +17,14 @@ MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
 [compat]
 ArrayLayouts = "0.1"
 BandedMatrices = "0.14"
-BlockArrays = "0.10"
-BlockBandedMatrices = "0.6"
-FillArrays = "0.8"
-InfiniteArrays = "0.4, 0.5"
-LazyArrays = "0.14.4, 0.15"
-LazyBandedMatrices = "0.0.1, 0.0.2, 0.0.3, 0.1"
+BlockArrays = "0.11"
+BlockBandedMatrices = "0.7.1"
+FillArrays = "0.8.4"
+InfiniteArrays = "0.6.1"
+LazyArrays = "0.15"
+LazyBandedMatrices = "0.2"
 MatrixFactorizations = "0.2"
-julia = "1.2"
+julia = "1.3"
 
 [extras]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

src/InfiniteLinearAlgebra.jl

@@ -1,33 +1,37 @@
 module InfiniteLinearAlgebra
 using BlockArrays, BlockBandedMatrices, BandedMatrices, LazyArrays, LazyBandedMatrices, 
-        FillArrays, InfiniteArrays, MatrixFactorizations, LinearAlgebra
+        FillArrays, InfiniteArrays, MatrixFactorizations, ArrayLayouts, LinearAlgebra
 
 import Base: +, -, *, /, \, ^, OneTo, getindex, promote_op, _unsafe_getindex, print_matrix_row, size, axes,
             AbstractMatrix, AbstractArray, Matrix, Array, Vector, AbstractVector, Slice,
-            show, getproperty
-import Base.Broadcast: BroadcastStyle, Broadcasted
+            show, getproperty, copy, map
+import Base.Broadcast: BroadcastStyle, Broadcasted, broadcasted
 
 import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, sublayout
 import BandedMatrices: BandedMatrix, _BandedMatrix, AbstractBandedMatrix, bandeddata, bandwidths, BandedColumns, bandedcolumns,
                         _default_banded_broadcast
 import FillArrays: AbstractFill, getindex_value      
-import InfiniteArrays: OneToInf, InfUnitRange, Infinity, InfStepRange, AbstractInfUnitRange                  
+import InfiniteArrays: OneToInf, InfUnitRange, Infinity, InfStepRange, AbstractInfUnitRange, InfAxes                  
 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,
-                    CachedMatrix, CachedArray, resizedata!, MemoryLayout, mulapplystyle, LmulStyle, RmulStyle,
-                    factorize, sub_materialize, LazyLayout, LazyArrayStyle,
-                    @lazymul, applylayout, ApplyLayout, PaddedLayout, materialize!, zero!
+                    resizedata!, MemoryLayout, mulapplystyle, LmulStyle, RmulStyle,
+                    factorize, sub_materialize, LazyLayout, LazyArrayStyle, lazy_getindex,
+                    @lazymul, applylayout, ApplyLayout, PaddedLayout, materialize!, zero!, MulAddStyle,
+                    LazyArray, LazyMatrix, LazyVector
 import MatrixFactorizations: ql, ql!, QLPackedQ, getL, getR, reflector!, reflectorApply!, QL, QR, QRPackedQ
 
-import BlockArrays: BlockSizes, cumulsizes, _find_block, AbstractBlockVecOrMat, sizes_from_blocks
+import BlockArrays: AbstractBlockVecOrMat, sizes_from_blocks, _length, BlockedUnitRange
 
 import BandedMatrices: BandedMatrix, bandwidths, AbstractBandedLayout, _banded_qr!, _banded_qr, _BandedMatrix
 
-import BlockBandedMatrices: _BlockSkylineMatrix, _BandedMatrix, AbstractBlockSizes, cumulsizes, _BlockSkylineMatrix, BlockSizes, blockstart, blockstride,
+import LazyBandedMatrices: MulBandedLayout, BroadcastBandedLayout
+
+import BlockBandedMatrices: _BlockSkylineMatrix, _BandedMatrix, _BlockSkylineMatrix, blockstart, blockstride,
         BlockSkylineSizes, BlockSkylineMatrix, BlockBandedMatrix, _BlockBandedMatrix, BlockTridiagonal
 
 
 LazyArrays.@lazymul BandedMatrix{<:Any,<:Any,<:OneToInf}
+*(A::BandedMatrix{<:Any,<:Any,<:OneToInf}, b::CachedVector) = apply(*,A,b)
 
 
 # BroadcastStyle(::Type{<:BandedMatrix{<:Any,<:Any,<:OneToInf}}) = LazyArrayStyle{2}()
@@ -56,7 +60,7 @@ export Vcat, Fill, ql, ql!, ∞, ContinuousSpectrumError, BlockTridiagonal
 include("banded/hessenbergq.jl")
 
 include("banded/infbanded.jl")
-include("blockbanded/infblocktridiagonal.jl")
+include("blockbanded/blockbanded.jl")
 include("banded/infqltoeplitz.jl")
 include("infql.jl")
 include("infqr.jl")
@@ -80,4 +84,28 @@ ApplyStyle(::typeof(*), ::Type{<:Adjoint{<:Any,<:BandedMatrix{<:Any,<:Any,<:OneT
     LazyArrayApplyStyle()
 ApplyStyle(::typeof(*), ::Type{<:Transpose{<:Any,<:BandedMatrix{<:Any,<:Any,<:OneToInf}}}, ::Type{<:AbstractArray}) =
     LazyArrayApplyStyle()
+
+ApplyStyle(::typeof(*), ::Type{<:BandedMatrix{<:Any,<:AbstractFill,<:OneToInf}}, ::Type{<:BandedMatrix{<:Any,<:AbstractFill,<:OneToInf}}) =
+    MulAddStyle()    
+
+#######
+# block broadcasted
+######
+
+const CumsumOneToInf2 = BroadcastArray{Int64,1,typeof(div),Tuple{BroadcastArray{Int64,1,typeof(*),Tuple{InfiniteArrays.OneToInf{Int64},InfiniteArrays.InfUnitRange{Int64}}},Int64}}
+BlockArrays.sortedunion(a::CumsumOneToInf2, ::CumsumOneToInf2) = a
+
+
+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) = ∞
+    
 end # module

src/banded/infbanded.jl

@@ -1,3 +1,9 @@
+###
+# Diagonal
+###
+
+getindex(D::Diagonal, k::InfAxes, j::InfAxes) = lazy_getindex(D, k, j)
+
 const TriToeplitz{T} = Tridiagonal{T,Fill{T,1,Tuple{OneToInf{Int}}}}
 const ConstRowMatrix{T} = ApplyMatrix{T,typeof(*),<:Tuple{<:AbstractVector,<:AbstractFill{<:Any,2,Tuple{OneTo{Int},OneToInf{Int}}}}}
 const PertConstRowMatrix{T} = Hcat{T,<:Tuple{Array{T},<:ConstRowMatrix{T}}}
@@ -322,6 +328,12 @@ _BandedMatrix(::PertToeplitzLayout, A::AbstractMatrix) =
 # end
 
 
+@inline sub_materialize(::MulBandedLayout, V, ::Tuple{InfAxes,InfAxes}) = V
+@inline sub_materialize(::BroadcastBandedLayout, V, ::Tuple{InfAxes,InfAxes}) = V
+@inline sub_materialize(::AbstractBandedLayout, V, ::Tuple{InfAxes,InfAxes}) = BandedMatrix(V)
+@inline sub_materialize(::BandedColumns, V, ::Tuple{InfAxes,InfAxes}) = BandedMatrix(V)
+
+
 ##
 # UniformScaling
 ##
@@ -335,4 +347,25 @@ for op in (:+, :-), Typ in (:(BandedMatrix{<:Any,<:Any,OneToInf{Int}}),
     end
 end
 
-_default_banded_broadcast(bc::Broadcasted, ::Tuple{<:OneToInf,<:Any}) = copy(Broadcasted{LazyArrayStyle{2}}(bc.f, bc.args))
+_default_banded_broadcast(bc::Broadcasted, ::Tuple{<:OneToInf,<:Any}) = copy(Broadcasted{LazyArrayStyle{2}}(bc.f, bc.args))
+
+
+###
+# BandedFill * BandedFill
+###
+
+copy(M::MulAdd{BandedMatrices.BandedColumns{FillLayout},BandedMatrices.BandedColumns{FillLayout},ZerosLayout}) =
+    _bandedfill_mul(M, axes(M.A), axes(M.B))
+
+_bandedfill_mul(M, _, _) = copyto!(similar(M), M)
+function _bandedfill_mul(M::MulAdd, ::Tuple{InfAxes,InfAxes}, ::Tuple{InfAxes,InfAxes})
+    A, B = M.A, M.B
+    Al,Au = bandwidths(A)
+    Bl,Bu = bandwidths(B)
+    l,u = Al+Bl,Au+Bu
+    m = min(Au+Al,Bl+Bu)+1
+    λ = getindex_value(bandeddata(A))*getindex_value(bandeddata(B))
+    ret = _BandedMatrix(Hcat(Array{typeof(λ)}(undef, l+u+1,u), [1:m-1; Fill(m,l+u-2m+3); m-1:-1:1]*Fill(λ,1,∞)), ∞, l, u)
+    mul!(view(ret, 1:l+u,1:u), view(A,1:l+u,1:u+Bl), view(B,1:u+Bl,1:u))
+    ret
+end

src/blockbanded/blockbanded.jl

@@ -0,0 +1,27 @@
+sizes_from_blocks(A::AbstractVector, ::Tuple{OneToInf{Int}}) = (map(length,A),)
+
+
+
+MemoryLayout(::Type{<:PseudoBlockArray{<:Any,N,<:Arr}}) where {N,Arr} = 
+    MemoryLayout(Arr)
+
+struct BlockLayout{LAY} <: MemoryLayout end
+
+MemoryLayout(::Type{<:BlockArray{<:Any,N,<:Arr}}) where {N,Arr} = 
+    blocklayout(MemoryLayout(Arr))
+
+blocklayout(_) = BlockLayout{UnknownLayout}()    
+blocklayout(::LazyLayout) = BlockLayout{LazyLayout}()    
+
+# for LazyLay in (:(BlockLayout{LazyLayout}), :(TriangularLayout{UPLO,UNIT,BlockLayout{LazyLayout}} where {UPLO,UNIT}))
+#     @eval begin
+#         combine_mul_styles(::$LazyLay) = LazyArrayApplyStyle()
+#         mulapplystyle(::QLayout, ::$LazyLay) = LazyArrayApplyStyle()
+#     end
+# end
+
+# BlockArrays.blockbroadcaststyle(::LazyArrayStyle{N}) where N = LazyArrayStyle{N}()
+
+
+include("infblocktridiagonal.jl")
+

src/blockbanded/infblocktridiagonal.jl

@@ -1,5 +1,5 @@
 const BlockTriPertToeplitz{T} = BlockMatrix{T,Tridiagonal{Matrix{T},Vcat{Matrix{T},1,Tuple{Vector{Matrix{T}},Fill{Matrix{T},1,Tuple{OneToInf{Int}}}}}},
-                                        BlockSizes{2,NTuple{2,Vcat{Int,1,Tuple{Int,Vcat{Int,1,Tuple{Vector{Int},InfStepRange{Int,Int}}}}}}}}
+                                        NTuple{2,BlockedUnitRange{Vcat{Int,1,Tuple{Vector{Int},InfStepRange{Int,Int}}}}}}
 
 
 function BlockTridiagonal(adjA::Adjoint{T,BlockTriPertToeplitz{T}}) where T
@@ -30,36 +30,9 @@ end
 
 *(a::AbstractVector, b::AbstractFill{<:Any,2,Tuple{OneTo{Int},OneToInf{Int}}}) = ApplyArray(*,a,b)
 
+sizes_from_blocks(A::Diagonal, ::NTuple{2,OneToInf{Int}}) = size.(A.diag, 1), size.(A.diag,2)
 
-sizes_from_blocks(A::AbstractVector, ::Tuple{OneToInf{Int}}) = BlockSizes((Vcat(1, 1 .+ cumsum(length.(A))),))
-function sizes_from_blocks(A::AbstractMatrix, ::Tuple{OneTo{Int}, OneToInf{Int}})
-    @assert size(A,1) == 1 
-    R,C = vec(size.(A,1)[:,1]), vec(size.(A,2))
-    BlockSizes((vcat(1, 1 .+ cumsum(R)), Vcat(1, 1 .+ cumsum(C))))
-end
-
-function sizes_from_blocks(A::Diagonal, ::NTuple{2,OneToInf{Int}}) 
-    sz = size.(A.diag, 1), size.(A.diag,2)
-    BlockSizes(Vcat.(1,(c -> 1 .+ c).(cumsum.(sz))))
-end
-
-function sizes_from_blocks(A::Tridiagonal, ::NTuple{2,OneToInf{Int}}) 
-    sz = size.(A.d, 1), size.(A.d,2)
-    BlockSizes(Vcat.(1,(c -> 1 .+ c).(cumsum.(sz))))
-end
-
-_find_block(cs::Number, i::Integer) = i ≤ cs ? 1 : 0
-function _find_block(cs::Vcat, i::Integer)
-    n = 0
-    for a in cs.args
-        i < first(a) && return n
-        if i ≤ last(a)
-            return _find_block(a, i) + n
-        end
-        n += length(a)
-    end 
-    return 0
-end
+sizes_from_blocks(A::Tridiagonal, ::NTuple{2,OneToInf{Int}}) = size.(A.d, 1), size.(A.d,2)
 
 print_matrix_row(io::IO,
         X::AbstractBlockVecOrMat, A::Vector,
@@ -98,7 +71,7 @@ function BlockSkylineSizes(A::BlockTriPertToeplitz, (l,u)::NTuple{2,Int})
         bs∞[k] = bs∞[k-1] .+ size(A[Block(N+1-u+k-1,N+1)],1)
     end
 
-    BlockSkylineSizes(blocksizes(A),
+    BlockSkylineSizes(axes(A),
                         _BandedMatrix(Hcat(block_starts.data, Vcat(adjoint.(bs∞)...)), ∞, l, u),
                         Vcat(block_strides, Fill(block_stride∞,∞)),
                         Fill(l,∞),Fill(u,∞))

src/infql.jl

@@ -175,16 +175,15 @@ function lmul!(adjA::Adjoint{<:Any,<:QLPackedQ{<:Any,<:InfBandedMatrix}}, B::Abs
     B
 end
 
-function (*)(A::QLPackedQ{T,<:InfBandedMatrix}, x::AbstractVector{S}) where {T,S}
+function _lmul_cache(A::AbstractMatrix{T}, x::AbstractVector{S}) where {T,S}
     TS = promote_op(matprod, T, S)
     lmul!(A, cache(convert(AbstractVector{TS},x)))
-end
-
-function (*)(A::Adjoint{T,<:QLPackedQ{T,<:InfBandedMatrix}}, x::AbstractVector{S}) where {T,S}
-    TS = promote_op(matprod, T, S)
-    lmul!(A, cache(convert(AbstractVector{TS},x)))
-end
+end    
 
+(*)(A::QLPackedQ{T,<:InfBandedMatrix}, x::AbstractVector) where {T} = _lmul_cache(A, x)
+(*)(A::Adjoint{T,<:QLPackedQ{T,<:InfBandedMatrix}}, x::AbstractVector) where {T} = _lmul_cache(A, x)
+(*)(A::QLPackedQ{T,<:InfBandedMatrix}, x::LazyVector) where {T} = _lmul_cache(A, x)
+(*)(A::Adjoint{T,<:QLPackedQ{T,<:InfBandedMatrix}}, x::LazyVector) where {T} = _lmul_cache(A, x)
 
 
 function blocktailiterate(c,a,b, d=c, e=a)

src/infqr.jl

@@ -158,15 +158,16 @@ function lmul!(adjA::Adjoint{<:Any,<:QRPackedQ{<:Any,<:AdaptiveQRFactors}}, B::C
     B
 end
 
-function (*)(A::QRPackedQ{T,<:AdaptiveQRFactors}, x::AbstractVector{S}) where {T,S}
+function _lmul_copymutable(A::AbstractMatrix{T}, x::AbstractVector{S}) where {T,S}
     TS = promote_op(matprod, T, S)
     lmul!(A, Base.copymutable(convert(AbstractVector{TS},x)))
 end
 
-function (*)(A::Adjoint{T,<:QRPackedQ{T,<:AdaptiveQRFactors}}, x::AbstractVector{S}) where {T,S}
-    TS = promote_op(matprod, T, S)
-    lmul!(A, Base.copymutable(convert(AbstractVector{TS},x)))
-end
+(*)(A::QRPackedQ{T,<:AdaptiveQRFactors}, x::AbstractVector) where {T} = _lmul_copymutable(A, x)
+(*)(A::Adjoint{T,<:QRPackedQ{T,<:AdaptiveQRFactors}}, x::AbstractVector) where {T} = _lmul_copymutable(A, x)
+(*)(A::QRPackedQ{T,<:AdaptiveQRFactors}, x::LazyVector) where {T} = _lmul_copymutable(A, x)
+(*)(A::Adjoint{T,<:QRPackedQ{T,<:AdaptiveQRFactors}}, x::LazyVector) where {T} = _lmul_copymutable(A, x)
+
 
 function ldiv!(R::UpperTriangular{<:Any,<:AdaptiveQRFactors}, B::CachedVector{<:Any,<:Any,<:Zeros{<:Any,1}}) 
     n = B.datasize[1]

test/runtests.jl

@@ -1,4 +1,4 @@
-using InfiniteLinearAlgebra, BlockBandedMatrices, BlockArrays, BandedMatrices, InfiniteArrays, FillArrays, LazyArrays, Test, MatrixFactorizations, LinearAlgebra, Random
+using InfiniteLinearAlgebra, BlockBandedMatrices, BlockArrays, BandedMatrices, InfiniteArrays, FillArrays, LazyArrays, Test, MatrixFactorizations, ArrayLayouts, LinearAlgebra, Random
 import InfiniteLinearAlgebra: qltail, toeptail, tailiterate , tailiterate!, tail_de, ql_X!,
                     InfToeplitz, PertToeplitz, TriToeplitz, InfBandedMatrix, 
                     rightasymptotics, QLHessenberg, ConstRows, PertConstRows, BandedToeplitzLayout, PertToeplitzLayout
@@ -8,6 +8,26 @@ import BlockBandedMatrices: isblockbanded, _BlockBandedMatrix
 import MatrixFactorizations: QLPackedQ
 import BandedMatrices: bandeddata, _BandedMatrix, BandedStyle
 import LazyArrays: colsupport, ApplyStyle, MemoryLayout, ApplyLayout, LazyArrayStyle
+import InfiniteArrays: OneToInf
+
+@testset "∞-banded" begin
+    D = Diagonal(Fill(2,∞))
+    B = D[1:∞,2:∞]
+    @test B isa BandedMatrix
+    @test B[1:10,1:10] == diagm(-1 => Fill(2,9))
+    @test B[1:∞,2:∞] isa BandedMatrix
+end
+
+@testset "∞-block arrays" begin
+    k = Base.OneTo.(Base.OneTo(∞))
+    n = Fill.(Base.OneTo(∞),Base.OneTo(∞))
+    @test broadcast(length,k) == map(length,k) == OneToInf()
+    @test broadcast(length,n) == map(length,n) == OneToInf()
+    b = mortar(Fill([1,2],∞))
+    @test blockaxes(b,1) === Block.(OneToInf())
+    @test b[Block(5)] == [1,2]
+    @test length(axes(b,1)) == last(axes(b,1)) == ∞
+end
 
 @testset "∞-Toeplitz and Pert-Toeplitz" begin
     A = BandedMatrix(1 => Fill(2im,∞), 2 => Fill(-1,∞), 3 => Fill(2,∞), -2 => Fill(-4,∞), -3 => Fill(-2im,∞))
@@ -28,6 +48,15 @@ import LazyArrays: colsupport, ApplyStyle, MemoryLayout, ApplyLayout, LazyArrayS
     @test MemoryLayout(typeof(V)) == PertToeplitzLayout()
     @test BandedMatrix(V) isa PertToeplitz
     @test A[2:∞,2:∞] isa PertToeplitz
+
+    @testset "InfBanded" begin
+        A = _BandedMatrix(Fill(2,4,∞),∞,2,1)
+        B = _BandedMatrix(Fill(3,2,∞),∞,-1,2)
+        @test mul(A,A) isa PertToeplitz
+        @test A*A isa PertToeplitz
+        @test (A*A)[1:20,1:20] == A[1:20,1:23]*A[1:23,1:20]
+        @test (A*B)[1:20,1:20] == A[1:20,1:23]*B[1:23,1:20]
+    end
 end
 
 
@@ -60,7 +89,7 @@ end
         @test At[1:10,1:10] ≈ transpose(A)[1:10,1:10] ≈ transpose(A[1:10,1:10])
 
         A = _BandedMatrix(Fill(1,4,∞),∞,1,2)
-        @test A*A isa ApplyArray
+        @test A*A isa BandedMatrix
         @test (A^2)[1:10,1:10] == (A*A)[1:10,1:10] == (A[1:100,1:100]^2)[1:10,1:10]
         @test (A^3)[1:10,1:10] == (A*A*A)[1:10,1:10] == (A[1:100,1:100]^3)[1:10,1:10]
     end
@@ -82,7 +111,6 @@ end
         @test (I + A)[1:100,1:100] == I+A[1:100,1:100]
         @test (I - A)[1:100,1:100] == I-A[1:100,1:100]
     end
-    
 
     @testset "Fill" begin
         A = _BandedMatrix(Ones(1,∞),∞,-1,1)
@@ -127,7 +155,15 @@ end
     end
 
     @testset "Triangle OP recurrences" begin
-        # mortar((n -> 1:n).(1:∞))
+        k = mortar(Base.OneTo.(1:∞))
+        n = mortar(Fill.(1:∞, 1:∞))
+        @test k[Block.(2:3)] isa BlockArray
+        @test n[Block.(2:3)] isa BlockArray
+        @test k[Block.(2:3)] == [1,2,1,2,3]
+        @test n[Block.(2:3)] == [2,2,3,3,3]
+        @test blocksize(BroadcastVector(exp,k)) == (∞,)
+        @test BroadcastVector(exp,k)[Block.(2:3)] == exp.([1,2,1,2,3])
+        # BroadcastVector(+,k,n)
     end
     # Multivariate OPs Corollary (3)
     # n = 5