From 66f0301e0eb275a21b088f5d7244117abbb59c5a Mon Sep 17 00:00:00 2001
From: MikaelSlevinsky <richard.slevinsky@umanitoba.ca>
Date: Mon, 31 Mar 2025 11:28:10 -0500
Subject: [PATCH 1/5] start adding bivariate gram

---
 Project.toml                     |   6 +
 src/BivariateGramMatrix.jl       | 565 +++++++++++++++++++++++++++++++
 src/FastTransforms.jl            |  15 +-
 test/bivariategrammatrixtests.jl |  25 ++
 test/runtests.jl                 |   1 +
 5 files changed, 609 insertions(+), 3 deletions(-)
 create mode 100644 src/BivariateGramMatrix.jl
 create mode 100644 test/bivariategrammatrixtests.jl

diff --git a/Project.toml b/Project.toml
index 2959561e..95a88b48 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,8 @@ version = "0.17"
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FastTransforms_jll = "34b6f7d7-08f9-5794-9e10-3819e4c7e49a"
@@ -17,12 +19,15 @@ Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 RecurrenceRelationships = "807425ed-42ea-44d6-a357-6771516d7b2c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
 [compat]
 AbstractFFTs = "1.0"
 ArrayLayouts = "1.10"
 BandedMatrices = "1.5"
+BlockArrays = "1"
+BlockBandedMatrices = "0.13"
 FFTW = "1.7"
 FastGaussQuadrature = "0.4, 0.5, 1"
 FastTransforms_jll = "0.6.2"
@@ -31,6 +36,7 @@ GenericFFT = "0.1"
 LazyArrays = "2.2"
 RecurrenceRelationships = "0.2"
 SpecialFunctions = "0.10, 1, 2"
+StaticArrays = "1.9"
 ToeplitzMatrices = "0.7.1, 0.8"
 julia = "1.7"
 
diff --git a/src/BivariateGramMatrix.jl b/src/BivariateGramMatrix.jl
new file mode 100644
index 00000000..f15f3cf9
--- /dev/null
+++ b/src/BivariateGramMatrix.jl
@@ -0,0 +1,565 @@
+"""
+Bivariate modified moments of ``P_{n,k}(x,y) = p_{n-k}(x) q_k(y)`` with respect to the separable measure ``{\\rm d}\\mu(x, y) = {\\rm d}\\mu_1(x) {\\rm\\,d}\\mu_2(y)``:
+
+```math
+    \\iint_{\\mathbb{R}^2} P_{n,k}(x, y) {\\rm\\,d}\\mu(x, y) = \\int_{\\mathbb{R}} p_{n-k}(x) {\\rm\\,d}\\mu_1(x) \\int_{\\mathbb{R}} q_{k}(y) {\\rm\\,d}\\mu_2(y).
+```
+"""
+function bivariatemoments(μ1::AbstractVector{T}, μ2::AbstractVector{T}) where T
+    @assert length(μ1) == length(μ2)
+    N = length(μ1)
+    μ = BlockVector{T}(undef, 1:N)
+    for n in 0:N-1
+        for k in 0:n
+            μ.blocks[n+1][k+1] = μ1[n-k+1]*μ2[k+1]
+        end
+    end
+    return μ
+end
+
+# These should live in BlockBandedMatrices.jl after PR #223
+
+@inline function inbands_viewblock(A::BandedBlockBandedMatrix, KJ::Block{2})
+    l,u = blockbandwidths(A)
+    K,J = KJ.n
+    BandedMatrices._BandedMatrix(view(A.data, Block(K-J+u+1, J)), length(axes(A,1)[Block(K)]), subblockbandwidths(A)...)
+end
+
+@inline function viewblock(A::BandedBlockBandedMatrix, KJ::Block{2})
+    @boundscheck checkbounds(A, KJ)
+    K,J = KJ.n
+    if -A.u ≤ K-J ≤ A.l
+        inbands_viewblock(A, KJ)
+    else
+        BandedMatrices._BandedMatrix(view(A.data, Block(1,1)), blocklengths(A,1)[Block(K)], (-40320,-40320))
+    end
+end
+
+@inline function viewblock(A::AbstractMatrix, KJ::Block{2})
+    view(A, KJ)
+end
+
+
+
+abstract type AbstractBivariateGramMatrix{T} <: AbstractBlockMatrix{T} end
+
+@inline issymmetric(G::AbstractBivariateGramMatrix) = true
+@inline isposdef(G::AbstractBivariateGramMatrix) = true
+
+
+struct BivariateGramMatrix{T, WT <: AbstractMatrix{T}, XT <: AbstractMatrix{T}, YT <: AbstractMatrix{T}} <: AbstractBivariateGramMatrix{T}
+    W::WT
+    X::XT
+    Y::YT
+    function BivariateGramMatrix{T, WT, XT, YT}(W::WT, X::XT, Y::YT) where {T, WT, XT, YT}
+        if size(W) ≠ size(X) ≠ size(Y)
+            throw(ArgumentError("Cannot construct a BivariateGramMatrix with W, X, and Y of different sizes."))
+        end
+        if !issymmetric(W)
+            throw(ArgumentError("Cannot construct a BivariateGramMatrix with a nonsymmetric W."))
+        end
+        if blockbandwidths(X) ≠ (1, 1)
+            throw(ArgumentError("Cannot construct a BivariateGramMatrix with a nonblocktridiagonal X."))
+        end
+        if blockbandwidths(Y) ≠ (1, 1)
+            throw(ArgumentError("Cannot construct a BivariateGramMatrix with a nonblocktridiagonal Y."))
+        end
+        new{T, WT, XT, YT}(W, X, Y)
+    end
+end
+
+BivariateGramMatrix(W::WT, X::XT, Y::YT) where {T, WT <: AbstractMatrix{T}, XT <: AbstractMatrix{T}, YT <: AbstractMatrix{T}} = BivariateGramMatrix{T, WT, XT, YT}(W, X, Y)
+
+@inline axes(G::BivariateGramMatrix) = axes(G.W)
+@inline getindex(G::BivariateGramMatrix, i::Integer, j::Integer) = G.W[i, j]
+@inline _blockindex_getindex(G::BivariateGramMatrix, bi::BlockIndex{2}) = _blockindex_getindex(G.W, bi)
+@inline blockbandwidths(G::BivariateGramMatrix) = blockbandwidths(G.W)
+@inline subblockbandwidths(G::BivariateGramMatrix) = subblockbandwidths(G.W)
+@inline MemoryLayout(G::BivariateGramMatrix) = MemoryLayout(G.W)
+@inline blockrowsupport(G::BivariateGramMatrix, j) = blockrowsupport(MemoryLayout(G), G.W, j)
+@inline blockcolsupport(G::BivariateGramMatrix, j) = blockcolsupport(MemoryLayout(G), G.W, j)
+
+BivariateGramMatrix(μ::AbstractVector{T}, X::XT, Y::YT) where {T, XT <: AbstractMatrix{T}, YT <: AbstractMatrix{T}} = BivariateGramMatrix(μ, X, Y, one(T))
+function BivariateGramMatrix(μ::AbstractVector{T}, X::XT, Y::YT, p0::T) where {T, XT <: AbstractMatrix{T}, YT <: AbstractMatrix{T}}
+    N = blocklength(μ)
+    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)
+    @assert subblockbandwidths(X) == (0, 0)
+    @assert subblockbandwidths(Y) == (1, 1)
+    W = BlockMatrix{T}(undef, 1:N, 1:N)
+    if n > 0
+        for m in 1:N
+            setblock!(W, p0*μ[Block(m, 1)], m, 1)
+        end
+    end
+    if n > 1
+        for m in 2:N-1
+            recurse!(W, X, Y, m)
+        end
+        symmetrize_block!(view(W, Block(2, 2)))
+    end
+    for n in 3:n
+        for m in n:N-n+1
+            recurse!(W, X, Y, m, n)
+        end
+        symmetrize_block!(view(W, Block(n, n)))
+    end
+    WN = BlockMatrix{T}(undef, 1:n, 1:n)
+    for j in 1:n
+        for k in j:n
+            setblock!(WN, viewblock(W, Block(k, j)), k, j)
+        end
+    end
+    WS = Symmetric(WN, :L)
+    XN = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (1, 1), (0, 0))
+    for j in 1:n
+        for k in max(1, j-1):min(n, j+1)
+            setblock!(XN, viewblock(X, Block(k, j)), k, j)
+        end
+    end
+    YN = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (1, 1), (1, 1))
+    for j in 1:n
+        for k in max(1, j-1):min(n, j+1)
+            setblock!(YN, viewblock(Y, Block(k, j)), k, j)
+        end
+    end
+    return BivariateGramMatrix(WS, XN, YN)
+    #return BivariateGramMatrix(Symmetric(W[Block.(1:n), Block.(1:n)], :L), X[Block.(1:n), Block.(1:n)], Y[Block.(1:n), Block.(1:n)])
+end
+
+function BivariateGramMatrix(μ::PaddedVector{T}, X::XT, Y::YT, p0::T) where {T, XT <: AbstractMatrix{T}, YT <: AbstractMatrix{T}}
+    N = blocklength(μ)
+    b = blocklength(μ.args[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)
+    @assert subblockbandwidths(X) == (0, 0)
+    @assert subblockbandwidths(Y) == (1, 1)
+    W = BandedBlockBandedMatrix{T}(undef, 1:N, 1:N, (b, 0), (b, b))
+    if n > 0
+        for m in 1:min(N, b+1)
+            setblock!(W, p0*μ[Block(m, 1)], m, 1)
+        end
+    end
+    if n > 1
+        for m in 2:min(N-1, b+2)
+            recurse!(W, X, Y, m)
+        end
+        symmetrize_block!(view(W, Block(2, 2)))
+    end
+    for n in 3:n
+        for m in n:min(N-n+1, b+n)
+            recurse!(W, X, Y, m, n)
+        end
+        symmetrize_block!(view(W, Block(n, n)))
+    end
+    WN = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (b, 0), (b, b))
+    for j in 1:n
+        for k in j:min(n, j+b)
+            setblock!(WN, viewblock(W, Block(k, j)), k, j)
+        end
+    end
+    WS = Symmetric(WN, :L)
+    XN = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (1, 1), (0, 0))
+    for j in 1:n
+        for k in max(1, j-1):min(n, j+1)
+            setblock!(XN, viewblock(X, Block(k, j)), k, j)
+        end
+    end
+    YN = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (1, 1), (1, 1))
+    for j in 1:n
+        for k in max(1, j-1):min(n, j+1)
+            setblock!(YN, viewblock(Y, Block(k, j)), k, j)
+        end
+    end
+    return BivariateGramMatrix(WS, XN, YN)
+    #return BivariateGramMatrix(Symmetric(W[Block.(1:n), Block.(1:n)], :L), X[Block.(1:n), Block.(1:n)], Y[Block.(1:n), Block.(1:n)])
+end
+
+function symmetrize_block!(W)
+    @assert size(W, 1) == size(W, 2)
+    @inbounds for j in 2:size(W, 2)
+        r = colrange(W, j)
+        ir = first(r):j-1
+        for i in ir
+            W[i, j] = W[j, i]
+        end
+    end
+end
+
+# n == 2
+function recurse!(W, X, Y, m)
+    # XW = X[Block(m-1, m)]'W[Block(m-1, 1)] + X[Block(m, m)]'W[Block(m, 1)] + X[Block(m+1, m)]'W[Block(m+1, 1)] - W[Block(m, 1)]*X[Block(1, 1)]
+    # Xn = Diagonal(X[Block(2, 1)])
+    # Wmn = zeros(T, m, 2)
+    # Wmn[:, 1] .= XW/Xn
+    Wa = view(W, Block(m-1, 1))
+    Wb = view(W, Block(m, 1))
+    Wc = view(W, Block(m+1, 1))
+    We = view(viewblock(W, Block(m, 2)), 1:m, 1:1)
+
+    Xa = view(X.data, Block(1, m))
+    Xb = view(X.data, Block(2, m))
+    Xc = view(X.data, Block(3, m))
+    Xd = view(X.data, Block(2, 1))
+    Xn = view(X.data, Block(3, 1))
+
+    r = colrange(Wa, 1)
+    jr = first(r):min(last(r), m-1)
+    @inbounds for j in jr
+        We[j, 1] = Xa[j]*Wa[j, 1]
+    end
+    We[m, 1] = 0
+    @inbounds for j in colrange(Wb, 1)
+        We[j, 1] += Xb[j]*Wb[j, 1]
+    end
+    r = colrange(Wc, 1)
+    jr = first(r):min(last(r), m)
+    @inbounds for j in jr
+        We[j, 1] += Xc[j]*Wc[j, 1]
+    end
+    @inbounds for j in colrange(Wb, 1)
+        We[j, 1] -= Wb[j, 1]*Xd[1]
+    end
+    @inbounds for j in colrange(We, 1)
+        We[j, 1] /= Xn[1]
+    end
+
+    # YWnm1 = Y[Block(m-1, m)]'view(Wa, :, 1) + Y[Block(m, m)]'view(Wb, :, 1) + Y[Block(m+1, m)]'view(Wc, :, 1) - Wb*view(Y[Block(1, 1)], :, 1)
+    # Yn = Y[Block(2, 1)]
+    # Wmn[:, 2] .= (YWnm1 .- Yn[1, 1].*Wmn[:, 1])./Yn[2, 1]
+    Ya = viewblock(Y, Block(m-1, m))
+    Yb = viewblock(Y, Block(m, m))
+    Yc = viewblock(Y, Block(m+1, m))
+    Yd = viewblock(Y, Block(1, 1))
+
+    We = view(viewblock(W, Block(m, 2)), 1:m, 2)
+    mul!(We, Ya', view(Wa, :, 1))
+    mul!(We, Yb', view(Wb, :, 1), 1, 1)
+    mul!(We, Yc', view(Wc, :, 1), 1, 1)
+    mul!(We, Wb, view(Yd, :, 1), -1, 1)
+    Yn = viewblock(Y, Block(2, 1))
+    We .-= Yn[1, 1].*view(viewblock(W, Block(m, 2)), 1:m, 1)
+    rdiv!(We, Yn[2, 1])
+end
+
+function recurse!(W, X, Y, m, n)
+    # XW = X[Block(m-1, m)]'W[Block(m-1, n-1)] + X[Block(m, m)]'W[Block(m, n-1)] + X[Block(m+1, m)]'W[Block(m+1, n-1)] - W[Block(m, n-1)]*X[Block(n-1, n-1)] - W[Block(m, n-2)]*X[Block(n-2, n-1)]
+    # Xn = Diagonal(X[Block(n, n-1)])
+    # Wmn = zeros(T, m, n)
+    # Wmn[:, 1:n-1] .= XW/Xn
+    Wa = view(W, Block(m-1, n-1))
+    Wb = view(W, Block(m, n-1))
+    Wc = view(W, Block(m+1, n-1))
+    Wd = view(W, Block(m, n-2))
+    We = view(viewblock(W, Block(m, n)), 1:m, 1:n-1)
+
+    Xa = view(X.data, Block(1, m))
+    Xb = view(X.data, Block(2, m))
+    Xc = view(X.data, Block(3, m))
+    Xd = view(X.data, Block(2, n-1))
+    Xe = view(X.data, Block(1, n-1))
+    Xn = view(X.data, Block(3, n-1))
+
+    @inbounds for k in 1:n-1
+        r = colrange(Wa, k)
+        jr = first(r):min(last(r), m-1)
+        for j in jr
+            We[j, k] = Xa[j]*Wa[j, k]
+        end
+        We[m, k] = 0
+    end
+    @inbounds for k in 1:n-1
+        for j in colrange(Wb, k)
+            We[j, k] += Xb[j]*Wb[j, k]
+        end
+    end
+    @inbounds for k in 1:n-1
+        r = colrange(Wc, k)
+        jr = first(r):min(last(r), m)
+        for j in jr
+            We[j, k] += Xc[j]*Wc[j, k]
+        end
+    end
+    @inbounds for k in 1:n-1
+        for j in colrange(Wb, k)
+            We[j, k] -= Wb[j, k]*Xd[k]
+        end
+    end
+    @inbounds for k in 1:n-2
+        for j in colrange(Wd, k)
+            We[j, k] -= Wd[j, k]*Xe[k]
+        end
+    end
+    @inbounds for k in 1:n-1
+        for j in colrange(We, k)
+            We[j, k] /= Xn[k]
+        end
+    end
+
+    # YWnm1 = Y[Block(m-1, m)]'view(Wa, :, n-1) + Y[Block(m, m)]'view(Wb, :, n-1) + Y[Block(m+1, m)]'view(Wc, :, n-1) - Wb*view(Y[Block(n-1, n-1)], :, n-1) - Wd*view(Y[Block(n-2, n-1)], :, n-1)
+    # Yn = Y[Block(n, n-1)]
+    # Wmn[:, n] .= (YWnm1 .- Yn[n-2, n-1].*Wmn[:, n-2] .- Yn[n-1, n-1].*Wmn[:, n-1])./Yn[n, n-1]
+    Ya = viewblock(Y, Block(m-1, m))
+    Yb = viewblock(Y, Block(m, m))
+    Yc = viewblock(Y, Block(m+1, m))
+    Yd = viewblock(Y, Block(n-1, n-1))
+    Ye = viewblock(Y, Block(n-2, n-1))
+
+    We = view(viewblock(W, Block(m, n)), 1:m, n)
+    mul!(We, Ya', view(Wa, :, n-1))
+    mul!(We, Yb', view(Wb, :, n-1), 1, 1)
+    mul!(We, Yc', view(Wc, :, n-1), 1, 1)
+    mul!(We, Wb, view(Yd, :, n-1), -1, 1)
+    mul!(We, Wd, view(Ye, :, n-1), -1, 1)
+    Yn = viewblock(Y, Block(n, n-1))
+    We .-= Yn[n-2, n-1].*view(viewblock(W, Block(m, n)), 1:m, n-2) .+ Yn[n-1, n-1].*view(viewblock(W, Block(m, n)), 1:m, n-1)
+    rdiv!(We, Yn[n, n-1])
+end
+
+
+
+compute_skew_generators(::Val{1}, G::BivariateGramMatrix) = compute_skew_generators(G.X, G)
+compute_skew_generators(::Val{2}, G::BivariateGramMatrix) = compute_skew_generators(G.Y, G)
+#
+# The computation of the skew-symmetric generators for the Gram matrix equation:
+# Z'W - WZ = GJGᵀ,
+# where Z = X or Y is a bivariate Jacobi matrix that is block tridiagonal.
+# We know that the left-hand side has zeros in the first (n-1) x (n-1) blocks.
+# Thus, determining G can be done with the matrix-matrix product in the last column.
+# This determines the first (n-1) blocks of the generators. To find the last,
+# we observe that the skew-symmetry in the last block leaves it underdetermined,
+# but it can be chosen to be strictly upper triangular.
+#
+function compute_skew_generators(Z::AbstractMatrix{T}, W::BivariateGramMatrix{T}) where T
+    n = blocksize(W, 1)
+    G = BlockMatrix{T}(undef_blocks, 1:n, SVector(n, n))
+    for j in 1:n-1
+        G.blocks[j, 1] = zeros(T, j, n)
+    end
+    G.blocks[n, 1] = Matrix{T}(I, n, n)
+    v = W[Block.(1:n), Block(n-1)]*Z[Block(n-1), Block(n)] + W[Block.(1:n), Block(n)]*Z[Block(n), Block(n)] - Z'W[Block.(1:n), Block(n)]
+    for j in 1:n-1
+        G.blocks[j, 2] = v[Block(j)]
+    end
+    G.blocks[n, 2] = triu!(v[Block(n)])
+    return G
+end
+
+
+
+
+
+struct BivariateChebyshevGramMatrix{T, BV <: AbstractBlockVector{T}, BS <: NTuple{2, AbstractUnitRange{Int}}, XT <: AbstractMatrix{T}, YT <: AbstractMatrix{T}} <: AbstractBivariateGramMatrix{T}
+    μ::BV
+    axes::BS
+    X::XT
+    Y::YT
+    function BivariateChebyshevGramMatrix{T, BV, BS, XT, YT}(μ::BV, axes::BS, X::XT, Y::YT) where {T, BV, BS, XT, YT}
+        if size(X) ≠ size(Y)
+            throw(ArgumentError("Cannot construct a BivariateChebyshevGramMatrix with X and Y of different sizes."))
+        end
+        if blockbandwidths(X) ≠ (1, 1)
+            throw(ArgumentError("Cannot construct a BivariateChebyshevGramMatrix with a nonblocktridiagonal X."))
+        end
+        if blockbandwidths(Y) ≠ (1, 1)
+            throw(ArgumentError("Cannot construct a BivariateChebyshevGramMatrix with a nonblocktridiagonal Y."))
+        end
+        new{T, BV, BS, XT, YT}(μ, axes, X, Y)
+    end
+end
+
+BivariateChebyshevGramMatrix(μ::BV, axes::BS, X::XT, Y::YT) where {T, BV <: AbstractBlockVector{T}, BS <: NTuple{2, AbstractUnitRange{Int}}, XT <: AbstractMatrix{T}, YT <: AbstractMatrix{T}} = BivariateChebyshevGramMatrix{T, BV, BS, XT, YT}(μ, axes, X, Y)
+
+function BivariateChebyshevGramMatrix(μ::AbstractBlockVector{T}) where T
+    n = (length(μ.blocks)+1)÷2
+    BivariateChebyshevGramMatrix(μ, (blockedrange(1:n), blockedrange(1:n)), _chebyshev_x(T, n), _chebyshev_y(T, n))
+end
+
+@inline axes(W::BivariateChebyshevGramMatrix) = W.axes
+
+Base.@propagate_inbounds getindex(W::BivariateChebyshevGramMatrix{T}, blockindex::BlockIndex{2}) where T = _blockindex_getindex(W, blockindex)
+
+@inline function _blockindex_getindex(W::BivariateChebyshevGramMatrix{T}, bi::BlockIndex{2}) where T
+    μ = W.μ
+    @boundscheck blockcheckbounds(W, Block(bi.I))
+    m, n = bi.I
+    j, k = bi.α
+    @boundscheck (1 ≤ j ≤ m) && (1 ≤ k ≤ n)
+    v = (μ.blocks[m+n-1][j+k-1]+μ.blocks[abs(m-j-n+k)+j+k-1][j+k-1]+μ.blocks[m-j+n-k+abs(j-k)+1][abs(j-k)+1]+μ.blocks[abs(m-j-n+k)+abs(j-k)+1][abs(j-k)+1])/4
+    return v
+end
+
+@inline function getindex(W::BivariateChebyshevGramMatrix{T}, i::Vararg{Integer, 2}) where T
+    @boundscheck checkbounds(W, i...)
+    @inbounds v = W[findblockindex.(axes(W), i)...]
+    return v
+end
+
+function _chebyshev_x(::Type{T}, n::Integer) where T
+    X = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (1, 1), (0, 0))
+    dat = X.data.blocks
+    @inbounds for j in 1:n*(n+1)÷2
+        dat[2, j] = zero(T)
+    end
+    @inbounds for j in 1:n-1
+        for k in 1:j-1
+            dat[3, k+j*(j-1)÷2] = one(T)/2
+        end
+        dat[3, j+j*(j-1)÷2] = one(T)
+    end
+    @inbounds for j in 2:n
+        for k in 1:j
+            dat[1, k+j*(j-1)÷2] = one(T)/2
+        end
+    end
+
+    return X
+end
+
+function _chebyshev_y(::Type{T}, n::Integer) where T
+    Y = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (1, 1), (1, 1))
+    dat = Y.data.blocks
+    @inbounds for k in 2:8
+        for j in 1:n*(n+1)÷2
+            dat[k, j] = zero(T)
+        end
+    end
+    @inbounds for j in 1:n-1
+        dat[9, 1+j*(j-1)÷2] = one(T)
+        for k in 2:j
+            dat[9, k+j*(j-1)÷2] = one(T)/2
+        end
+    end
+    @inbounds for j in 2:n
+        for k in 2:j
+            dat[1, k+j*(j-1)÷2] = one(T)/2
+        end
+    end
+
+    return Y
+end
+
+#
+# The computation of the skew-symmetric generators for the bivariate Chebyshev--Gram matrix equation:
+# Z'W-WZ = GJGᵀ,
+# where Z = X or Y is a bivariate Chebyshev Jacobi matrix that is block tridiagonal.
+# Since the bivariate Chebyhsev-Gram matrix is special, these formulas compute G
+# without the use of X or Y (thus without matrix multiplication).
+#
+function compute_skew_generators(::Val{1}, W::BivariateChebyshevGramMatrix{T}) where T
+    N = blocksize(W, 1)
+    G = BlockMatrix{T}(undef_blocks, 1:N, SVector(N, N))
+    for j in 1:N-1
+        G.blocks[j, 1] = zeros(T, j, N)
+    end
+    G.blocks[N, 1] = Matrix{T}(I, N, N)
+    μ = W.μ
+    @inbounds for m in 1:N-1
+        GB = zeros(T, m, N)
+        for k in 1:N-1
+            for j in 1:m
+                GB[j, k] = -(μ.blocks[m+N][j+k-1]+μ.blocks[abs(m-j-N-1+k)+j+k-1][j+k-1]+μ.blocks[m-j+N+1-k+abs(j-k)+1][abs(j-k)+1]+μ.blocks[abs(m-j-N-1+k)+abs(j-k)+1][abs(j-k)+1])/8
+            end
+        end
+        for j in 1:m
+            GB[j, N] = -(μ.blocks[m+N][j+N-1]+μ.blocks[abs(m-j-1)+j+N-1][j+N-1]+μ.blocks[m-j+1+abs(j-N)+1][abs(j-N)+1]+μ.blocks[abs(m-j-1)+abs(j-N)+1][abs(j-N)+1])/4
+        end
+        G.blocks[m, 2] = GB
+    end
+    GB = zeros(T, N, N)
+    @inbounds for k in 1:N-1
+        for j in 1:N
+            if abs(k-j-1)+j+k-1 < 2N
+                GB[j, k] -= μ.blocks[abs(k-j-1)+j+k-1][j+k-1]
+            end
+            if 2N-j+1-k+abs(j-k)+1 < 2N
+                GB[j, k] -= μ.blocks[2N-j+1-k+abs(j-k)+1][abs(j-k)+1]
+            end
+            if abs(k-j-1)+abs(j-k)+1 < 2N
+                GB[j, k] -= μ.blocks[abs(k-j-1)+abs(j-k)+1][abs(j-k)+1]
+            end
+            GB[j, k] /= 8
+        end
+    end
+    @inbounds for j in 1:N
+        if abs(N-j-1)+j+N-1 < 2N
+            GB[j, N] -= μ.blocks[abs(N-j-1)+j+N-1][j+N-1]
+        end
+        if N-j+1+abs(j-N)+1 < 2N
+            GB[j, N] -= μ.blocks[N-j+1+abs(j-N)+1][abs(j-N)+1]
+        end
+        if abs(N-j-1)+abs(j-N)+1 < 2N
+            GB[j, N] -= μ.blocks[abs(N-j-1)+abs(j-N)+1][abs(j-N)+1]
+        end
+        GB[j, N] /= 4
+    end
+    @inbounds for k in 1:N
+        for j in 1:k
+            GB[j, k] -= GB[k, j]
+            GB[k, j] = zero(T)
+        end
+    end
+    G.blocks[N, 2] = GB
+
+    G
+end
+
+function compute_skew_generators(::Val{2}, W::BivariateChebyshevGramMatrix{T}) where T
+    N = blocksize(W, 1)
+    G = BlockMatrix{T}(undef_blocks, 1:N, SVector(N, N))
+    for j in 1:N-1
+        G.blocks[j, 1] = zeros(T, j, N)
+    end
+    G.blocks[N, 1] = Matrix{T}(I, N, N)
+    μ = W.μ
+    @inbounds for m in 1:N-1
+        GB = zeros(T, m, N)
+        for j in 1:m
+            GB[j, 1] = -(μ.blocks[m+N][j+1]+μ.blocks[abs(m-j-N+1)+j+1][j+1]+μ.blocks[m-j+N-1+abs(j-2)+1][abs(j-2)+1]+μ.blocks[abs(m-j-N+1)+abs(j-2)+1][abs(j-2)+1])/4
+        end
+        for k in 2:N
+            for j in 1:m
+                GB[j, k] = -(μ.blocks[m+N][j+k]+μ.blocks[abs(m-j-N+k)+j+k][j+k]+μ.blocks[m-j+N-k+abs(j-k-1)+1][abs(j-k-1)+1]+μ.blocks[abs(m-j-N+k)+abs(j-k-1)+1][abs(j-k-1)+1])/8
+            end
+        end
+        G.blocks[m, 2] = GB
+    end
+    GB = zeros(T, N, N)
+    @inbounds for j in 1:N
+        if abs(1-j)+j+1 < 2N
+            GB[j, 1] -= μ.blocks[abs(1-j)+j+1][j+1]
+        end
+        if 2N-j-1+abs(j-2)+1 < 2N
+            GB[j, 1] -= μ.blocks[2N-j-1+abs(j-2)+1][abs(j-2)+1]
+        end
+        if abs(1-j)+abs(j-2)+1 < 2N
+            GB[j, 1] -= μ.blocks[abs(1-j)+abs(j-2)+1][abs(j-2)+1]
+        end
+        GB[j, 1] /= 4
+    end
+    @inbounds for k in 2:N
+        for j in 1:N
+            if abs(k-j)+j+k < 2N
+                GB[j, k] -= μ.blocks[abs(k-j)+j+k][j+k]
+            end
+            if 2N-j-k+abs(j-k-1)+1 < 2N
+                GB[j, k] -= μ.blocks[2N-j-k+abs(j-k-1)+1][abs(j-k-1)+1]
+            end
+            if abs(k-j)+abs(j-k-1)+1 < 2N
+                GB[j, k] -= μ.blocks[abs(k-j)+abs(j-k-1)+1][abs(j-k-1)+1]
+            end
+            GB[j, k] /= 8
+        end
+    end
+    @inbounds for k in 1:N
+        for j in 1:k
+            GB[j, k] -= GB[k, j]
+            GB[k, j] = zero(T)
+        end
+    end
+    G.blocks[N, 2] = GB
+
+    G
+end
diff --git a/src/FastTransforms.jl b/src/FastTransforms.jl
index a9324735..9526feee 100644
--- a/src/FastTransforms.jl
+++ b/src/FastTransforms.jl
@@ -1,13 +1,14 @@
 module FastTransforms
 
-using ArrayLayouts, BandedMatrices, FastGaussQuadrature, FillArrays, LazyArrays, LinearAlgebra,
-      SpecialFunctions, ToeplitzMatrices, RecurrenceRelationships
+using ArrayLayouts, BandedMatrices, BlockArrays, BlockBandedMatrices,
+      FastGaussQuadrature, FillArrays, LazyArrays, LinearAlgebra,
+      SpecialFunctions, StaticArrays, ToeplitzMatrices, RecurrenceRelationships
 
 using AbstractFFTs
 using FFTW
 using GenericFFT
 
-import Base: convert, unsafe_convert, eltype, ndims, adjoint, transpose, show,
+import Base: axes, convert, unsafe_convert, eltype, ndims, adjoint, transpose, show, showerror,
              *, \, inv, length, size, view, getindex, tail, OneTo
 
 import Base.GMP: Limb
@@ -23,6 +24,10 @@ import ArrayLayouts: rowsupport, colsupport, LayoutMatrix, MemoryLayout, Abstrac
 
 import BandedMatrices: bandwidths, BandedLayout
 
+import BlockArrays: _blockindex_getindex, blockrowsupport, blockcolsupport
+
+import BlockBandedMatrices: blockbandwidths, subblockbandwidths, blockcolrange
+
 import FFTW: dct, dct!, idct, idct!, plan_dct!, plan_idct!,
              plan_dct, plan_idct, fftwNumber
 
@@ -105,6 +110,10 @@ export GramMatrix, ChebyshevGramMatrix
 
 include("GramMatrix.jl")
 
+export BivariateGramMatrix, BivariateChebyshevGramMatrix
+
+include("BivariateGramMatrix.jl")
+
 export weightedhermitetransform, iweightedhermitetransform
 
 include("hermite.jl")
diff --git a/test/bivariategrammatrixtests.jl b/test/bivariategrammatrixtests.jl
new file mode 100644
index 00000000..d1d18e30
--- /dev/null
+++ b/test/bivariategrammatrixtests.jl
@@ -0,0 +1,25 @@
+using BlockArrays, FastTransforms, LazyArrays, LinearAlgebra, Test
+
+import FastTransforms: chebyshevmoments1, chebyshevabsmoments1, bivariatemoments
+
+@testset "BivariateGramMatrix" begin
+    n = 20
+    # w(x, y) = 1 + |x| + |y| - |xy|
+    μ1 = chebyshevmoments1(Float64, 2n-1)
+    μ2 = chebyshevabsmoments1(Float64, 2n-1)
+    μ = bivariatemoments(μ1, μ1) + bivariatemoments(μ2, μ1) + bivariatemoments(μ1, μ2) - bivariatemoments(μ2, μ2)
+
+    X = FastTransforms._chebyshev_x(Float64, 2n-1)
+    Y = FastTransforms._chebyshev_y(Float64, 2n-1)
+
+    W = BivariateGramMatrix(μ, X, Y)
+    @test issymmetric(W)
+    @test isposdef(W)
+    WC = BivariateChebyshevGramMatrix(μ)
+    @test W ≈ WC
+
+    R = cholesky(W).U
+    RC = cholesky(WC).U
+
+    @test R ≈ RC
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 36c95de8..f3305b1b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -12,4 +12,5 @@ include("toeplitzplanstests.jl")
 include("toeplitzhankeltests.jl")
 include("toeplitzplushankeltests.jl")
 include("grammatrixtests.jl")
+include("bivariategrammatrixtests.jl")
 include("arraystests.jl")

From c90bfd3123e6e3adaba5958ff7b7b319f2d4cbb3 Mon Sep 17 00:00:00 2001
From: MikaelSlevinsky <richard.slevinsky@umanitoba.ca>
Date: Tue, 1 Apr 2025 13:01:19 -0500
Subject: [PATCH 2/5] refactor for AbstractBlockVector moments
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

- remove most instances of μ.blocks in favour of μ[BlockIndex(n, k)] which should exist for any AbstractBlockVector

- have the bivariatemoments constructor return a BlockedVector

- the rationale is that the block-analogue of a zero-padded vector isa BlockedVector{T, <: PaddedVector{T}} rather than a BlockVector

- remove X and Y from BivariateChebyshevGramMatrix. the whole point for it is to be a light type
---
 src/BivariateGramMatrix.jl       | 196 +++++++++++++++----------------
 src/FastTransforms.jl            |   4 +-
 src/GramMatrix.jl                |   2 +-
 test/bivariategrammatrixtests.jl |  12 ++
 4 files changed, 113 insertions(+), 101 deletions(-)

diff --git a/src/BivariateGramMatrix.jl b/src/BivariateGramMatrix.jl
index f15f3cf9..8a346106 100644
--- a/src/BivariateGramMatrix.jl
+++ b/src/BivariateGramMatrix.jl
@@ -8,10 +8,10 @@ Bivariate modified moments of ``P_{n,k}(x,y) = p_{n-k}(x) q_k(y)`` with respect
 function bivariatemoments(μ1::AbstractVector{T}, μ2::AbstractVector{T}) where T
     @assert length(μ1) == length(μ2)
     N = length(μ1)
-    μ = BlockVector{T}(undef, 1:N)
+    μ = BlockedVector{T}(undef, 1:N)
     for n in 0:N-1
         for k in 0:n
-            μ.blocks[n+1][k+1] = μ1[n-k+1]*μ2[k+1]
+            μ[BlockIndex(n+1, k+1)] = μ1[n-k+1]*μ2[k+1]
         end
     end
     return μ
@@ -79,8 +79,8 @@ BivariateGramMatrix(W::WT, X::XT, Y::YT) where {T, WT <: AbstractMatrix{T}, XT <
 @inline blockrowsupport(G::BivariateGramMatrix, j) = blockrowsupport(MemoryLayout(G), G.W, j)
 @inline blockcolsupport(G::BivariateGramMatrix, j) = blockcolsupport(MemoryLayout(G), G.W, j)
 
-BivariateGramMatrix(μ::AbstractVector{T}, X::XT, Y::YT) where {T, XT <: AbstractMatrix{T}, YT <: AbstractMatrix{T}} = BivariateGramMatrix(μ, X, Y, one(T))
-function BivariateGramMatrix(μ::AbstractVector{T}, X::XT, Y::YT, p0::T) where {T, XT <: AbstractMatrix{T}, YT <: AbstractMatrix{T}}
+BivariateGramMatrix(μ::AbstractBlockVector{T}, X::XT, Y::YT) where {T, XT <: AbstractMatrix{T}, YT <: AbstractMatrix{T}} = BivariateGramMatrix(μ, X, Y, one(T))
+function BivariateGramMatrix(μ::AbstractBlockVector{T}, X::XT, Y::YT, p0::T) where {T, XT <: AbstractMatrix{T}, YT <: AbstractMatrix{T}}
     N = blocklength(μ)
     n = (N+1)÷2
     @assert N == blocksize(X, 1) == blocksize(X, 2) == blocksize(Y, 1) == blocksize(Y, 2)
@@ -90,7 +90,7 @@ function BivariateGramMatrix(μ::AbstractVector{T}, X::XT, Y::YT, p0::T) where {
     W = BlockMatrix{T}(undef, 1:N, 1:N)
     if n > 0
         for m in 1:N
-            setblock!(W, p0*μ[Block(m, 1)], m, 1)
+            W[Block(m, 1)] = p0*μ[Block(m, 1)]
         end
     end
     if n > 1
@@ -108,27 +108,27 @@ function BivariateGramMatrix(μ::AbstractVector{T}, X::XT, Y::YT, p0::T) where {
     WN = BlockMatrix{T}(undef, 1:n, 1:n)
     for j in 1:n
         for k in j:n
-            setblock!(WN, viewblock(W, Block(k, j)), k, j)
+            WN[Block(k, j)] = viewblock(W, Block(k, j))
         end
     end
     WS = Symmetric(WN, :L)
     XN = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (1, 1), (0, 0))
     for j in 1:n
         for k in max(1, j-1):min(n, j+1)
-            setblock!(XN, viewblock(X, Block(k, j)), k, j)
+            XN[Block(k, j)] = viewblock(X, Block(k, j))
         end
     end
     YN = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (1, 1), (1, 1))
     for j in 1:n
         for k in max(1, j-1):min(n, j+1)
-            setblock!(YN, viewblock(Y, Block(k, j)), k, j)
+            YN[Block(k, j)] = viewblock(Y, Block(k, j))
         end
     end
     return BivariateGramMatrix(WS, XN, YN)
     #return BivariateGramMatrix(Symmetric(W[Block.(1:n), Block.(1:n)], :L), X[Block.(1:n), Block.(1:n)], Y[Block.(1:n), Block.(1:n)])
 end
 
-function BivariateGramMatrix(μ::PaddedVector{T}, X::XT, Y::YT, p0::T) where {T, XT <: AbstractMatrix{T}, YT <: AbstractMatrix{T}}
+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
     n = (N+1)÷2
@@ -139,7 +139,7 @@ function BivariateGramMatrix(μ::PaddedVector{T}, X::XT, Y::YT, p0::T) where {T,
     W = BandedBlockBandedMatrix{T}(undef, 1:N, 1:N, (b, 0), (b, b))
     if n > 0
         for m in 1:min(N, b+1)
-            setblock!(W, p0*μ[Block(m, 1)], m, 1)
+            W[Block(m, 1)] = p0*μ[Block(m, 1)]
         end
     end
     if n > 1
@@ -157,20 +157,20 @@ function BivariateGramMatrix(μ::PaddedVector{T}, X::XT, Y::YT, p0::T) where {T,
     WN = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (b, 0), (b, b))
     for j in 1:n
         for k in j:min(n, j+b)
-            setblock!(WN, viewblock(W, Block(k, j)), k, j)
+            WN[Block(k, j)] = viewblock(W, Block(k, j))
         end
     end
     WS = Symmetric(WN, :L)
     XN = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (1, 1), (0, 0))
     for j in 1:n
         for k in max(1, j-1):min(n, j+1)
-            setblock!(XN, viewblock(X, Block(k, j)), k, j)
+            XN[Block(k, j)] = viewblock(X, Block(k, j))
         end
     end
     YN = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (1, 1), (1, 1))
     for j in 1:n
         for k in max(1, j-1):min(n, j+1)
-            setblock!(YN, viewblock(Y, Block(k, j)), k, j)
+            YN[Block(k, j)] = viewblock(Y, Block(k, j))
         end
     end
     return BivariateGramMatrix(WS, XN, YN)
@@ -351,95 +351,47 @@ end
 
 
 
-struct BivariateChebyshevGramMatrix{T, BV <: AbstractBlockVector{T}, BS <: NTuple{2, AbstractUnitRange{Int}}, XT <: AbstractMatrix{T}, YT <: AbstractMatrix{T}} <: AbstractBivariateGramMatrix{T}
+struct BivariateChebyshevGramMatrix{T, BV <: AbstractBlockVector{T}, BS <: NTuple{2, AbstractUnitRange{Int}}} <: AbstractBivariateGramMatrix{T}
     μ::BV
     axes::BS
-    X::XT
-    Y::YT
-    function BivariateChebyshevGramMatrix{T, BV, BS, XT, YT}(μ::BV, axes::BS, X::XT, Y::YT) where {T, BV, BS, XT, YT}
-        if size(X) ≠ size(Y)
-            throw(ArgumentError("Cannot construct a BivariateChebyshevGramMatrix with X and Y of different sizes."))
-        end
-        if blockbandwidths(X) ≠ (1, 1)
-            throw(ArgumentError("Cannot construct a BivariateChebyshevGramMatrix with a nonblocktridiagonal X."))
-        end
-        if blockbandwidths(Y) ≠ (1, 1)
-            throw(ArgumentError("Cannot construct a BivariateChebyshevGramMatrix with a nonblocktridiagonal Y."))
-        end
-        new{T, BV, BS, XT, YT}(μ, axes, X, Y)
-    end
 end
 
-BivariateChebyshevGramMatrix(μ::BV, axes::BS, X::XT, Y::YT) where {T, BV <: AbstractBlockVector{T}, BS <: NTuple{2, AbstractUnitRange{Int}}, XT <: AbstractMatrix{T}, YT <: AbstractMatrix{T}} = BivariateChebyshevGramMatrix{T, BV, BS, XT, YT}(μ, axes, X, Y)
-
 function BivariateChebyshevGramMatrix(μ::AbstractBlockVector{T}) where T
-    n = (length(μ.blocks)+1)÷2
-    BivariateChebyshevGramMatrix(μ, (blockedrange(1:n), blockedrange(1:n)), _chebyshev_x(T, n), _chebyshev_y(T, n))
+    n = (blocklength(μ)+1)÷2
+    BivariateChebyshevGramMatrix(μ, (blockedrange(1:n), blockedrange(1:n)))
 end
 
-@inline axes(W::BivariateChebyshevGramMatrix) = W.axes
+@inline axes(G::BivariateChebyshevGramMatrix) = G.axes
 
-Base.@propagate_inbounds getindex(W::BivariateChebyshevGramMatrix{T}, blockindex::BlockIndex{2}) where T = _blockindex_getindex(W, blockindex)
+Base.@propagate_inbounds getindex(G::BivariateChebyshevGramMatrix{T}, blockindex::BlockIndex{2}) where T = _blockindex_getindex(G, blockindex)
 
-@inline function _blockindex_getindex(W::BivariateChebyshevGramMatrix{T}, bi::BlockIndex{2}) where T
-    μ = W.μ
-    @boundscheck blockcheckbounds(W, Block(bi.I))
+@inline function _blockindex_getindex(G::BivariateChebyshevGramMatrix{T}, bi::BlockIndex{2}) where T
+    @boundscheck blockcheckbounds(G, Block(bi.I))
     m, n = bi.I
     j, k = bi.α
     @boundscheck (1 ≤ j ≤ m) && (1 ≤ k ≤ n)
-    v = (μ.blocks[m+n-1][j+k-1]+μ.blocks[abs(m-j-n+k)+j+k-1][j+k-1]+μ.blocks[m-j+n-k+abs(j-k)+1][abs(j-k)+1]+μ.blocks[abs(m-j-n+k)+abs(j-k)+1][abs(j-k)+1])/4
+    μ = G.μ
+    v = (μ[BlockIndex(m+n-1, j+k-1)]+μ[BlockIndex(abs(m-j-n+k)+j+k-1, j+k-1)]+μ[BlockIndex(m-j+n-k+abs(j-k)+1, abs(j-k)+1)]+μ[BlockIndex(abs(m-j-n+k)+abs(j-k)+1, abs(j-k)+1)])/4
     return v
 end
 
-@inline function getindex(W::BivariateChebyshevGramMatrix{T}, i::Vararg{Integer, 2}) where T
-    @boundscheck checkbounds(W, i...)
-    @inbounds v = W[findblockindex.(axes(W), i)...]
+@inline function getindex(G::BivariateChebyshevGramMatrix{T}, i::Vararg{Integer, 2}) where T
+    @boundscheck checkbounds(G, i...)
+    @inbounds v = G[findblockindex.(axes(G), i)...]
     return v
 end
 
-function _chebyshev_x(::Type{T}, n::Integer) where T
-    X = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (1, 1), (0, 0))
-    dat = X.data.blocks
-    @inbounds for j in 1:n*(n+1)÷2
-        dat[2, j] = zero(T)
-    end
-    @inbounds for j in 1:n-1
-        for k in 1:j-1
-            dat[3, k+j*(j-1)÷2] = one(T)/2
-        end
-        dat[3, j+j*(j-1)÷2] = one(T)
-    end
-    @inbounds for j in 2:n
-        for k in 1:j
-            dat[1, k+j*(j-1)÷2] = one(T)/2
-        end
-    end
-
-    return X
+@inline function blockbandwidths(G::BivariateChebyshevGramMatrix{T, <: BlockedVector{T, <: PaddedVector{T}}}) where T
+    N = length(G.μ.blocks.args[2])
+    b = ceil(Int, (-1+sqrt(1+8N))/2) - 1
+    return (b, b)
 end
-
-function _chebyshev_y(::Type{T}, n::Integer) where T
-    Y = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (1, 1), (1, 1))
-    dat = Y.data.blocks
-    @inbounds for k in 2:8
-        for j in 1:n*(n+1)÷2
-            dat[k, j] = zero(T)
-        end
-    end
-    @inbounds for j in 1:n-1
-        dat[9, 1+j*(j-1)÷2] = one(T)
-        for k in 2:j
-            dat[9, k+j*(j-1)÷2] = one(T)/2
-        end
-    end
-    @inbounds for j in 2:n
-        for k in 2:j
-            dat[1, k+j*(j-1)÷2] = one(T)/2
-        end
-    end
-
-    return Y
+@inline function subblockbandwidths(G::BivariateChebyshevGramMatrix{T, <: BlockedVector{T, <: PaddedVector{T}}}) where T
+    N = length(G.μ.blocks.args[2])
+    b = ceil(Int, (-1+sqrt(1+8N))/2) - 1
+    return (b, b)
 end
+@inline MemoryLayout(G::BivariateChebyshevGramMatrix{T, <: BlockedVector{T, <: PaddedVector{T}}}) where T = BandedBlockBandedLayout()
 
 #
 # The computation of the skew-symmetric generators for the bivariate Chebyshev--Gram matrix equation:
@@ -460,11 +412,11 @@ function compute_skew_generators(::Val{1}, W::BivariateChebyshevGramMatrix{T}) w
         GB = zeros(T, m, N)
         for k in 1:N-1
             for j in 1:m
-                GB[j, k] = -(μ.blocks[m+N][j+k-1]+μ.blocks[abs(m-j-N-1+k)+j+k-1][j+k-1]+μ.blocks[m-j+N+1-k+abs(j-k)+1][abs(j-k)+1]+μ.blocks[abs(m-j-N-1+k)+abs(j-k)+1][abs(j-k)+1])/8
+                GB[j, k] = -(μ[BlockIndex(m+N, j+k-1)]+μ[BlockIndex(abs(m-j-N-1+k)+j+k-1, j+k-1)]+μ[BlockIndex(m-j+N+1-k+abs(j-k)+1, abs(j-k)+1)]+μ[BlockIndex(abs(m-j-N-1+k)+abs(j-k)+1, abs(j-k)+1)])/8
             end
         end
         for j in 1:m
-            GB[j, N] = -(μ.blocks[m+N][j+N-1]+μ.blocks[abs(m-j-1)+j+N-1][j+N-1]+μ.blocks[m-j+1+abs(j-N)+1][abs(j-N)+1]+μ.blocks[abs(m-j-1)+abs(j-N)+1][abs(j-N)+1])/4
+            GB[j, N] = -(μ[BlockIndex(m+N, j+N-1)]+μ[BlockIndex(abs(m-j-1)+j+N-1, j+N-1)]+μ[BlockIndex(m-j+1+abs(j-N)+1, abs(j-N)+1)]+μ[BlockIndex(abs(m-j-1)+abs(j-N)+1, abs(j-N)+1)])/4
         end
         G.blocks[m, 2] = GB
     end
@@ -472,26 +424,26 @@ function compute_skew_generators(::Val{1}, W::BivariateChebyshevGramMatrix{T}) w
     @inbounds for k in 1:N-1
         for j in 1:N
             if abs(k-j-1)+j+k-1 < 2N
-                GB[j, k] -= μ.blocks[abs(k-j-1)+j+k-1][j+k-1]
+                GB[j, k] -= μ[BlockIndex(abs(k-j-1)+j+k-1, j+k-1)]
             end
             if 2N-j+1-k+abs(j-k)+1 < 2N
-                GB[j, k] -= μ.blocks[2N-j+1-k+abs(j-k)+1][abs(j-k)+1]
+                GB[j, k] -= μ[BlockIndex(2N-j+1-k+abs(j-k)+1, abs(j-k)+1)]
             end
             if abs(k-j-1)+abs(j-k)+1 < 2N
-                GB[j, k] -= μ.blocks[abs(k-j-1)+abs(j-k)+1][abs(j-k)+1]
+                GB[j, k] -= μ[BlockIndex(abs(k-j-1)+abs(j-k)+1, abs(j-k)+1)]
             end
             GB[j, k] /= 8
         end
     end
     @inbounds for j in 1:N
         if abs(N-j-1)+j+N-1 < 2N
-            GB[j, N] -= μ.blocks[abs(N-j-1)+j+N-1][j+N-1]
+            GB[j, N] -= μ[BlockIndex(abs(N-j-1)+j+N-1, j+N-1)]
         end
         if N-j+1+abs(j-N)+1 < 2N
-            GB[j, N] -= μ.blocks[N-j+1+abs(j-N)+1][abs(j-N)+1]
+            GB[j, N] -= μ[BlockIndex(N-j+1+abs(j-N)+1, abs(j-N)+1)]
         end
         if abs(N-j-1)+abs(j-N)+1 < 2N
-            GB[j, N] -= μ.blocks[abs(N-j-1)+abs(j-N)+1][abs(j-N)+1]
+            GB[j, N] -= μ[BlockIndex(abs(N-j-1)+abs(j-N)+1, abs(j-N)+1)]
         end
         GB[j, N] /= 4
     end
@@ -517,11 +469,11 @@ function compute_skew_generators(::Val{2}, W::BivariateChebyshevGramMatrix{T}) w
     @inbounds for m in 1:N-1
         GB = zeros(T, m, N)
         for j in 1:m
-            GB[j, 1] = -(μ.blocks[m+N][j+1]+μ.blocks[abs(m-j-N+1)+j+1][j+1]+μ.blocks[m-j+N-1+abs(j-2)+1][abs(j-2)+1]+μ.blocks[abs(m-j-N+1)+abs(j-2)+1][abs(j-2)+1])/4
+            GB[j, 1] = -(μ[BlockIndex(m+N, j+1)]+μ[BlockIndex(abs(m-j-N+1)+j+1, j+1)]+μ[BlockIndex(m-j+N-1+abs(j-2)+1, abs(j-2)+1)]+μ[BlockIndex(abs(m-j-N+1)+abs(j-2)+1, abs(j-2)+1)])/4
         end
         for k in 2:N
             for j in 1:m
-                GB[j, k] = -(μ.blocks[m+N][j+k]+μ.blocks[abs(m-j-N+k)+j+k][j+k]+μ.blocks[m-j+N-k+abs(j-k-1)+1][abs(j-k-1)+1]+μ.blocks[abs(m-j-N+k)+abs(j-k-1)+1][abs(j-k-1)+1])/8
+                GB[j, k] = -(μ[BlockIndex(m+N, j+k)]+μ[BlockIndex(abs(m-j-N+k)+j+k, j+k)]+μ[BlockIndex(m-j+N-k+abs(j-k-1)+1, abs(j-k-1)+1)]+μ[BlockIndex(abs(m-j-N+k)+abs(j-k-1)+1, abs(j-k-1)+1)])/8
             end
         end
         G.blocks[m, 2] = GB
@@ -529,26 +481,26 @@ function compute_skew_generators(::Val{2}, W::BivariateChebyshevGramMatrix{T}) w
     GB = zeros(T, N, N)
     @inbounds for j in 1:N
         if abs(1-j)+j+1 < 2N
-            GB[j, 1] -= μ.blocks[abs(1-j)+j+1][j+1]
+            GB[j, 1] -= μ[BlockIndex(abs(1-j)+j+1, j+1)]
         end
         if 2N-j-1+abs(j-2)+1 < 2N
-            GB[j, 1] -= μ.blocks[2N-j-1+abs(j-2)+1][abs(j-2)+1]
+            GB[j, 1] -= μ[BlockIndex(2N-j-1+abs(j-2)+1, abs(j-2)+1)]
         end
         if abs(1-j)+abs(j-2)+1 < 2N
-            GB[j, 1] -= μ.blocks[abs(1-j)+abs(j-2)+1][abs(j-2)+1]
+            GB[j, 1] -= μ[BlockIndex(abs(1-j)+abs(j-2)+1, abs(j-2)+1)]
         end
         GB[j, 1] /= 4
     end
     @inbounds for k in 2:N
         for j in 1:N
             if abs(k-j)+j+k < 2N
-                GB[j, k] -= μ.blocks[abs(k-j)+j+k][j+k]
+                GB[j, k] -= μ[BlockIndex(abs(k-j)+j+k, j+k)]
             end
             if 2N-j-k+abs(j-k-1)+1 < 2N
-                GB[j, k] -= μ.blocks[2N-j-k+abs(j-k-1)+1][abs(j-k-1)+1]
+                GB[j, k] -= μ[BlockIndex(2N-j-k+abs(j-k-1)+1, abs(j-k-1)+1)]
             end
             if abs(k-j)+abs(j-k-1)+1 < 2N
-                GB[j, k] -= μ.blocks[abs(k-j)+abs(j-k-1)+1][abs(j-k-1)+1]
+                GB[j, k] -= μ[BlockIndex(abs(k-j)+abs(j-k-1)+1, abs(j-k-1)+1)]
             end
             GB[j, k] /= 8
         end
@@ -563,3 +515,51 @@ function compute_skew_generators(::Val{2}, W::BivariateChebyshevGramMatrix{T}) w
 
     G
 end
+
+
+
+## Move to tests?
+
+function _chebyshev_x(::Type{T}, n::Integer) where T
+    X = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (1, 1), (0, 0))
+    dat = X.data.blocks
+    @inbounds for j in 1:n*(n+1)÷2
+        dat[2, j] = zero(T)
+    end
+    @inbounds for j in 1:n-1
+        for k in 1:j-1
+            dat[3, k+j*(j-1)÷2] = one(T)/2
+        end
+        dat[3, j+j*(j-1)÷2] = one(T)
+    end
+    @inbounds for j in 2:n
+        for k in 1:j
+            dat[1, k+j*(j-1)÷2] = one(T)/2
+        end
+    end
+
+    return X
+end
+
+function _chebyshev_y(::Type{T}, n::Integer) where T
+    Y = BandedBlockBandedMatrix{T}(undef, 1:n, 1:n, (1, 1), (1, 1))
+    dat = Y.data.blocks
+    @inbounds for k in 2:8
+        for j in 1:n*(n+1)÷2
+            dat[k, j] = zero(T)
+        end
+    end
+    @inbounds for j in 1:n-1
+        dat[9, 1+j*(j-1)÷2] = one(T)
+        for k in 2:j
+            dat[9, k+j*(j-1)÷2] = one(T)/2
+        end
+    end
+    @inbounds for j in 2:n
+        for k in 2:j
+            dat[1, k+j*(j-1)÷2] = one(T)/2
+        end
+    end
+
+    return Y
+end
diff --git a/src/FastTransforms.jl b/src/FastTransforms.jl
index 9526feee..f0b93a86 100644
--- a/src/FastTransforms.jl
+++ b/src/FastTransforms.jl
@@ -8,7 +8,7 @@ using AbstractFFTs
 using FFTW
 using GenericFFT
 
-import Base: axes, convert, unsafe_convert, eltype, ndims, adjoint, transpose, show, showerror,
+import Base: axes, convert, unsafe_convert, eltype, ndims, adjoint, transpose, show,
              *, \, inv, length, size, view, getindex, tail, OneTo
 
 import Base.GMP: Limb
@@ -26,7 +26,7 @@ import BandedMatrices: bandwidths, BandedLayout
 
 import BlockArrays: _blockindex_getindex, blockrowsupport, blockcolsupport
 
-import BlockBandedMatrices: blockbandwidths, subblockbandwidths, blockcolrange
+import BlockBandedMatrices: blockbandwidths, subblockbandwidths, blockcolrange, BandedBlockBandedLayout
 
 import FFTW: dct, dct!, idct, idct!, plan_dct!, plan_idct!,
              plan_dct, plan_idct, fftwNumber
diff --git a/src/GramMatrix.jl b/src/GramMatrix.jl
index 75f0cc22..f5d54284 100644
--- a/src/GramMatrix.jl
+++ b/src/GramMatrix.jl
@@ -277,7 +277,7 @@ See also [`GramMatrix`](@ref) for the general case.
 """
 function ChebyshevGramMatrix(μ::V) where {T, V <: AbstractVector{T}}
     n = (length(μ)+1)÷2
-    ChebyshevGramMatrix{T, V}(μ, n)
+    ChebyshevGramMatrix(μ, n)
 end
 
 @inline size(G::ChebyshevGramMatrix) = (G.n, G.n)
diff --git a/test/bivariategrammatrixtests.jl b/test/bivariategrammatrixtests.jl
index d1d18e30..4af5e38c 100644
--- a/test/bivariategrammatrixtests.jl
+++ b/test/bivariategrammatrixtests.jl
@@ -22,4 +22,16 @@ import FastTransforms: chebyshevmoments1, chebyshevabsmoments1, bivariatemoments
     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
+
+    Gy = FastTransforms.compute_skew_generators(Val(2), W)
+    GyC = FastTransforms.compute_skew_generators(Val(2), WC)
+    @test Gy ≈ GyC
+
+    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'
 end

From 688994cf3952c99e2b8ac8b848418b3a7d17d97d Mon Sep 17 00:00:00 2001
From: MikaelSlevinsky <richard.slevinsky@umanitoba.ca>
Date: Tue, 1 Apr 2025 14:41:33 -0500
Subject: [PATCH 3/5] cleanup blockbanded constructors

---
 src/BivariateGramMatrix.jl       | 17 ++++++++++++++++-
 test/bivariategrammatrixtests.jl | 23 +++++++++++++++++------
 2 files changed, 33 insertions(+), 7 deletions(-)

diff --git a/src/BivariateGramMatrix.jl b/src/BivariateGramMatrix.jl
index 8a346106..641dffe6 100644
--- a/src/BivariateGramMatrix.jl
+++ b/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)
diff --git a/test/bivariategrammatrixtests.jl b/test/bivariategrammatrixtests.jl
index 4af5e38c..1475df64 100644
--- a/test/bivariategrammatrixtests.jl
+++ b/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

From 04135c56c2851787393c938d2a7ab7711c6dd893 Mon Sep 17 00:00:00 2001
From: MikaelSlevinsky <richard.slevinsky@umanitoba.ca>
Date: Mon, 12 May 2025 11:42:13 -0500
Subject: [PATCH 4/5] convert W's data to a BlockedMatrix

---
 src/BivariateGramMatrix.jl | 6 ++++--
 src/FastTransforms.jl      | 2 +-
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/BivariateGramMatrix.jl b/src/BivariateGramMatrix.jl
index 641dffe6..d8665bfc 100644
--- a/src/BivariateGramMatrix.jl
+++ b/src/BivariateGramMatrix.jl
@@ -90,6 +90,8 @@ BivariateGramMatrix(W::WT, X::XT, Y::YT) where {T, WT <: AbstractMatrix{T}, XT <
 @inline blockbandwidths(G::BivariateGramMatrix) = blockbandwidths(G.W)
 @inline subblockbandwidths(G::BivariateGramMatrix) = subblockbandwidths(G.W)
 @inline MemoryLayout(G::BivariateGramMatrix) = MemoryLayout(G.W)
+@inline symmetricdata(G::BivariateGramMatrix) = symmetricdata(G.W)
+@inline symmetricuplo(G::BivariateGramMatrix) = symmetricuplo(G.W)
 @inline blockrowsupport(G::BivariateGramMatrix, j) = blockrowsupport(MemoryLayout(G), G.W, j)
 @inline blockcolsupport(G::BivariateGramMatrix, j) = blockcolsupport(MemoryLayout(G), G.W, j)
 
@@ -101,7 +103,7 @@ function BivariateGramMatrix(μ::AbstractBlockVector{T}, X::XT, Y::YT, p0::T) wh
     @assert blockbandwidths(X) == blockbandwidths(Y) == (1, 1)
     @assert subblockbandwidths(X) == (0, 0)
     @assert subblockbandwidths(Y) == (1, 1)
-    W = BlockMatrix{T}(undef, 1:N, 1:N)
+    W = BlockedMatrix{T}(undef, 1:N, 1:N)
     if n > 0
         for m in 1:N
             W[Block(m, 1)] = p0*μ[Block(m, 1)]
@@ -119,7 +121,7 @@ function BivariateGramMatrix(μ::AbstractBlockVector{T}, X::XT, Y::YT, p0::T) wh
         end
         symmetrize_block!(view(W, Block(n, n)))
     end
-    WN = BlockMatrix{T}(undef, 1:n, 1:n)
+    WN = BlockedMatrix{T}(undef, 1:n, 1:n)
     for j in 1:n
         for k in j:n
             WN[Block(k, j)] = viewblock(W, Block(k, j))
diff --git a/src/FastTransforms.jl b/src/FastTransforms.jl
index f0b93a86..d90f8569 100644
--- a/src/FastTransforms.jl
+++ b/src/FastTransforms.jl
@@ -20,7 +20,7 @@ import AbstractFFTs: Plan, ScaledPlan,
                      fftshift, ifftshift, rfft_output_size, brfft_output_size,
                      normalization
 
-import ArrayLayouts: rowsupport, colsupport, LayoutMatrix, MemoryLayout, AbstractBandedLayout
+import ArrayLayouts: rowsupport, colsupport, LayoutMatrix, MemoryLayout, AbstractBandedLayout, symmetricdata, symmetricuplo
 
 import BandedMatrices: bandwidths, BandedLayout
 

From da1dc1783c3a250bc5c3c4b3ab11bf522ef2132d Mon Sep 17 00:00:00 2001
From: MikaelSlevinsky <richard.slevinsky@umanitoba.ca>
Date: Thu, 15 May 2025 14:36:53 -0500
Subject: [PATCH 5/5] fixup recurrence for tridiag-block-tridiag Y

---
 src/BivariateGramMatrix.jl | 80 +++++++++++++++++++++++++++++---------
 1 file changed, 61 insertions(+), 19 deletions(-)

diff --git a/src/BivariateGramMatrix.jl b/src/BivariateGramMatrix.jl
index d8665bfc..cb2026ad 100644
--- a/src/BivariateGramMatrix.jl
+++ b/src/BivariateGramMatrix.jl
@@ -167,6 +167,7 @@ function BivariateGramMatrix(μ::BlockedVector{T, <: PaddedVector{T}}, X::XT, Y:
     end
     for n in 3:n
         for m in n:min(N-n+1, b+n)
+            #println("m = $m, n = $n")
             recurse!(W, X, Y, m, n)
         end
         symmetrize_block!(view(W, Block(n, n)))
@@ -214,7 +215,7 @@ function recurse!(W, X, Y, m)
     Wa = view(W, Block(m-1, 1))
     Wb = view(W, Block(m, 1))
     Wc = view(W, Block(m+1, 1))
-    We = view(viewblock(W, Block(m, 2)), 1:m, 1:1)
+    We = viewblock(W, Block(m, 2))
 
     Xa = view(X.data, Block(1, m))
     Xb = view(X.data, Block(2, m))
@@ -250,15 +251,37 @@ function recurse!(W, X, Y, m)
     Yb = viewblock(Y, Block(m, m))
     Yc = viewblock(Y, Block(m+1, m))
     Yd = viewblock(Y, Block(1, 1))
-
-    We = view(viewblock(W, Block(m, 2)), 1:m, 2)
-    mul!(We, Ya', view(Wa, :, 1))
-    mul!(We, Yb', view(Wb, :, 1), 1, 1)
-    mul!(We, Yc', view(Wc, :, 1), 1, 1)
-    mul!(We, Wb, view(Yd, :, 1), -1, 1)
     Yn = viewblock(Y, Block(2, 1))
-    We .-= Yn[1, 1].*view(viewblock(W, Block(m, 2)), 1:m, 1)
-    rdiv!(We, Yn[2, 1])
+    for j in colrange(We, 2)
+        if j == 1
+            if j == size(Ya, 2)-1
+                We[j, 2] = Ya[j, j]*Wa[j, 1]
+            else
+                We[j, 2] = Ya[j, j]*Wa[j, 1] + Ya[j+1, j]*Wa[j+1, 1]
+            end
+        elseif j == size(Ya, 2)
+            We[j, 2] = Ya[j-1, j]*Wa[j-1, 1]
+        elseif j == size(Ya, 2)-1
+            We[j, 2] = Ya[j-1, j]*Wa[j-1, 1] + Ya[j, j]*Wa[j, 1]
+        else
+            We[j, 2] = Ya[j-1, j]*Wa[j-1, 1] + Ya[j, j]*Wa[j, 1] + Ya[j+1, j]*Wa[j+1, 1]
+        end
+        if j == 1
+            We[j, 2] += Yb[j, j]*Wb[j, 1] + Yb[j+1, j]*Wb[j+1, 1]
+        elseif j == size(Yb, 2)
+            We[j, 2] += Yb[j-1, j]*Wb[j-1, 1] + Yb[j, j]*Wb[j, 1]
+        else
+            We[j, 2] += Yb[j-1, j]*Wb[j-1, 1] + Yb[j, j]*Wb[j, 1] + Yb[j+1, j]*Wb[j+1, 1]
+        end
+        if j == 1
+            We[j, 2] += Yc[j, j]*Wc[j, 1] + Yc[j+1, j]*Wc[j+1, 1]
+        else
+            We[j, 2] += Yc[j-1, j]*Wc[j-1, 1] + Yc[j, j]*Wc[j, 1] + Yc[j+1, j]*Wc[j+1, 1]
+        end
+        We[j, 2] -= Wb[j, 1]*Yd[1, 1]
+        We[j, 2] -= Yn[1, 1]*We[j, 1]
+        We[j, 2] /= Yn[2, 1]
+    end
 end
 
 function recurse!(W, X, Y, m, n)
@@ -270,7 +293,7 @@ function recurse!(W, X, Y, m, n)
     Wb = view(W, Block(m, n-1))
     Wc = view(W, Block(m+1, n-1))
     Wd = view(W, Block(m, n-2))
-    We = view(viewblock(W, Block(m, n)), 1:m, 1:n-1)
+    We = viewblock(W, Block(m, n))
 
     Xa = view(X.data, Block(1, m))
     Xb = view(X.data, Block(2, m))
@@ -323,16 +346,35 @@ function recurse!(W, X, Y, m, n)
     Yc = viewblock(Y, Block(m+1, m))
     Yd = viewblock(Y, Block(n-1, n-1))
     Ye = viewblock(Y, Block(n-2, n-1))
-
-    We = view(viewblock(W, Block(m, n)), 1:m, n)
-    mul!(We, Ya', view(Wa, :, n-1))
-    mul!(We, Yb', view(Wb, :, n-1), 1, 1)
-    mul!(We, Yc', view(Wc, :, n-1), 1, 1)
-    mul!(We, Wb, view(Yd, :, n-1), -1, 1)
-    mul!(We, Wd, view(Ye, :, n-1), -1, 1)
     Yn = viewblock(Y, Block(n, n-1))
-    We .-= Yn[n-2, n-1].*view(viewblock(W, Block(m, n)), 1:m, n-2) .+ Yn[n-1, n-1].*view(viewblock(W, Block(m, n)), 1:m, n-1)
-    rdiv!(We, Yn[n, n-1])
+
+    for j in colrange(We, n)
+        if j == 1
+            We[j, n] = Ya[j, j]*Wa[j, n-1] + Ya[j+1, j]*Wa[j+1, n-1]
+        elseif j == size(Ya, 2)
+            We[j, n] = Ya[j-1, j]*Wa[j-1, n-1]
+        elseif j == size(Ya, 2)-1
+            We[j, n] = Ya[j-1, j]*Wa[j-1, n-1] + Ya[j, j]*Wa[j, n-1]
+        else
+            We[j, n] = Ya[j-1, j]*Wa[j-1, n-1] + Ya[j, j]*Wa[j, n-1] + Ya[j+1, j]*Wa[j+1, n-1]
+        end
+        if j == 1
+            We[j, n] += Yb[j, j]*Wb[j, n-1] + Yb[j+1, j]*Wb[j+1, n-1]
+        elseif j == size(Yb, 2)
+            We[j, n] += Yb[j-1, j]*Wb[j-1, n-1] + Yb[j, j]*Wb[j, n-1]
+        else
+            We[j, n] += Yb[j-1, j]*Wb[j-1, n-1] + Yb[j, j]*Wb[j, n-1] + Yb[j+1, j]*Wb[j+1, n-1]
+        end
+        if j == 1
+            We[j, n] += Yc[j, j]*Wc[j, n-1] + Yc[j+1, j]*Wc[j+1, n-1]
+        else
+            We[j, n] += Yc[j-1, j]*Wc[j-1, n-1] + Yc[j, j]*Wc[j, n-1] + Yc[j+1, j]*Wc[j+1, n-1]
+        end
+        We[j, n] -= Wb[j, n-2]*Yd[n-2, n-1] + Wb[j, n-1]*Yd[n-1, n-1]
+        We[j, n] -= Wd[j, n-2]*Ye[n-2, n-1]
+        We[j, n] -= Yn[n-2, n-1]*We[j, n-2] + Yn[n-1, n-1]*We[j, n-1]
+        We[j, n] /= Yn[n, n-1]
+    end
 end