cleanup blockbanded constructors

Author: MikaelSlevinsky

src/BivariateGramMatrix.jl

@@ -17,6 +17,20 @@ function bivariatemoments(μ1::AbstractVector{T}, μ2::AbstractVector{T}) where
     return μ
 end
 
+function bivariatemoments(μ1::PaddedVector{T}, μ2::PaddedVector{T}) where T
+    @assert length(μ1) == length(μ2)
+    N = length(μ1)
+    b = length(μ1.args[2])+length(μ2.args[2])-1
+    v = Vector{T}(undef, b*(b+1)÷2)
+    for n in 0:b-1
+        for k in 0:n
+            v[n*(n+1)÷2+k+1] = μ1[n-k+1]*μ2[k+1]
+        end
+    end
+    μ = BlockedVector(PaddedVector(v, N*(N+1)÷2), 1:N)
+    return μ
+end
+
 # These should live in BlockBandedMatrices.jl after PR #223
 
 @inline function inbands_viewblock(A::BandedBlockBandedMatrix, KJ::Block{2})
@@ -130,7 +144,8 @@ end
 
 function BivariateGramMatrix(μ::BlockedVector{T, <: PaddedVector{T}}, X::XT, Y::YT, p0::T) where {T, XT <: AbstractMatrix{T}, YT <: AbstractMatrix{T}}
     N = blocklength(μ)
-    b = blocklength(μ.args[2])-1
+    bb = length(μ.blocks.args[2])
+    b = ceil(Int, (-1+sqrt(1+8bb))/2) - 1
     n = (N+1)÷2
     @assert N == blocksize(X, 1) == blocksize(X, 2) == blocksize(Y, 1) == blocksize(Y, 2)
     @assert blockbandwidths(X) == blockbandwidths(Y) == (1, 1)

test/bivariategrammatrixtests.jl

@@ -1,4 +1,4 @@
-using BlockArrays, FastTransforms, LazyArrays, LinearAlgebra, Test
+using BlockArrays, BlockBandedMatrices, FastTransforms, LazyArrays, LinearAlgebra, Test
 
 import FastTransforms: chebyshevmoments1, chebyshevabsmoments1, bivariatemoments
 
@@ -18,11 +18,6 @@ import FastTransforms: chebyshevmoments1, chebyshevabsmoments1, bivariatemoments
     WC = BivariateChebyshevGramMatrix(μ)
     @test W ≈ WC
 
-    R = cholesky(W).U
-    RC = cholesky(WC).U
-
-    @test R ≈ RC
-
     Gx = FastTransforms.compute_skew_generators(Val(1), W)
     GxC = FastTransforms.compute_skew_generators(Val(1), WC)
     @test Gx ≈ GxC
@@ -34,4 +29,20 @@ import FastTransforms: chebyshevmoments1, chebyshevabsmoments1, bivariatemoments
     J = [zeros(n, n) Matrix{Float64}(I, n, n); Matrix{Float64}(-I, n, n) zeros(n, n)]
     @test W.X'W-W*W.X ≈ Gx*J*Gx'
     @test W.Y'W-W*W.Y ≈ Gy*J*Gy'
+
+    R = cholesky(W).U
+    RC = cholesky(WC).U
+
+    @test R ≈ RC
+
+    μ1 = PaddedVector(1 ./ [1,2,3], 2n-1)
+    μ2 = PaddedVector(1 ./ [1,2,3,4,5,6], 2n-1)
+    μ = bivariatemoments(μ1, μ2)
+    μ̂ = bivariatemoments(Vector(μ1), Vector(μ2))
+    @test μ ≈ μ̂
+
+    W = BivariateGramMatrix(μ, X, Y)
+    WC = BivariateChebyshevGramMatrix(μ)
+    @test blockbandwidths(W) == blockbandwidths(WC) == subblockbandwidths(W) == subblockbandwidths(WC) == (7, 7)
+    @test W ≈ WC
 end