From e643dfa71e35f40177e269ce783ea7d7aa4653f3 Mon Sep 17 00:00:00 2001
From: DanielVandH <danj.vandenheuvel@gmail.com>
Date: Sun, 30 Jun 2024 22:38:37 +0100
Subject: [PATCH 01/13] Implement BidiagonalConjugationData

---
 Project.toml                            |  4 +-
 src/InfiniteLinearAlgebra.jl            |  3 +-
 src/banded/bidiagonalconjugationdata.jl | 85 +++++++++++++++++++++++++
 test/runtests.jl                        |  1 +
 test/test_bidiagonalconjugationdata.jl  | 33 ++++++++++
 5 files changed, 123 insertions(+), 3 deletions(-)
 create mode 100644 src/banded/bidiagonalconjugationdata.jl
 create mode 100644 test/test_bidiagonalconjugationdata.jl

diff --git a/Project.toml b/Project.toml
index 030395a..3dc14a1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "InfiniteLinearAlgebra"
 uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
-version = "0.8.3"
+version = "0.8.4"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -23,8 +23,8 @@ BandedMatrices = "1.7.2"
 BlockArrays = "1.0"
 BlockBandedMatrices = "0.13"
 FillArrays = "1.0"
-Infinities = "0.1"
 InfiniteArrays = "0.14"
+Infinities = "0.1"
 LazyArrays = "2.0"
 LazyBandedMatrices = "0.10"
 LinearAlgebra = "1"
diff --git a/src/InfiniteLinearAlgebra.jl b/src/InfiniteLinearAlgebra.jl
index 19203f3..f75e59e 100644
--- a/src/InfiniteLinearAlgebra.jl
+++ b/src/InfiniteLinearAlgebra.jl
@@ -21,7 +21,7 @@ import ArrayLayouts: AbstractBandedLayout, AbstractQLayout, AdjQRPackedQLayout,
 
 import BandedMatrices: BandedColumns, BandedMatrix, BandedMatrix, _BandedMatrix, AbstractBandedMatrix,
                        _BandedMatrix, _BandedMatrix, _banded_qr, _banded_qr!, _default_banded_broadcast, banded_chol!,
-                       banded_similar, bandedcolumns, bandeddata, bandwidths, bandwidths
+                       banded_similar, bandedcolumns, bandeddata, bandwidths
 
 import BlockArrays: AbstractBlockLayout, BlockLayout, BlockSlice, BlockSlice1, BlockedOneTo,
                     blockcolsupport, sizes_from_blocks, OneToCumsum, AbstractBlockedUnitRange
@@ -143,5 +143,6 @@ include("infql.jl")
 include("infqr.jl")
 include("inful.jl")
 include("infcholesky.jl")
+include("banded/bidiagonalconjugationdata.jl")
 
 end # module
diff --git a/src/banded/bidiagonalconjugationdata.jl b/src/banded/bidiagonalconjugationdata.jl
new file mode 100644
index 0000000..ed8a3d8
--- /dev/null
+++ b/src/banded/bidiagonalconjugationdata.jl
@@ -0,0 +1,85 @@
+"""
+    BidiagonalConjugationData{T, MU, MC} <: LazyMatrix{T}
+
+Struct for efficiently representing the matrix product `A = inv(U)XV`,
+assuming that 
+
+- `A` is upper bidiagonal, 
+- `U` is upper Hessenberg,
+- `X` is banded,
+- `V` is banded.
+
+None of these properties are checked internally. It is the user's responsibility
+to ensure these properties hold and that the product `inv(U)XV` is indeed bidiagonal.
+
+# Fields 
+- `U`: The upper Hessenberg matrix. 
+- `C`: The matrix product `XV`.
+- `dv`: A vector giving the diagonal of `A`. 
+- `ev`: A vector giving the superdiagonal of `A`.
+
+The vectors `dv` and `ev` grow on demand from `getindex` and should not be 
+used directly. Simply treat 
+
+    A = BidiagonalConjugationData(U, X, V)
+
+as you would an upper bidiagonal matrix.
+"""
+struct BidiagonalConjugationData{T,MU,MC} <: LazyMatrix{T}
+    U::MU
+    C::MC
+    dv::Vector{T}
+    ev::Vector{T}
+end
+function BidiagonalConjugationData(U::MU, X::MX, V::MV) where {MU,MX,MV}
+    C = X * V
+    T = promote_type(typeof(inv(U[begin])), eltype(U), eltype(C)) # include inv so that we can't get Ints
+    dv, ev = T[], T[]
+    return BidiagonalConjugationData{T,MU,typeof(C)}(U, C, dv, ev)
+end
+MemoryLayout(::Type{<:BidiagonalConjugationData}) = BidiagonalLayout{LazyLayout,LazyLayout}()
+bandwidths(A::BidiagonalConjugationData) = (0, 1)
+size(A::BidiagonalConjugationData) = (ℵ₀, ℵ₀)
+axes(A::BidiagonalConjugationData) = (OneToInf(), OneToInf())
+Base.eltype(A::Type{<:BidiagonalConjugationData{T}}) where {T} = T
+
+copy(A::BidiagonalConjugationData) = BidiagonalConjugationData(copy(A.U), copy(A.C), copy(A.dv), copy(A.ev))
+copy(A::Adjoint{T,<:BidiagonalConjugationData}) where {T} = copy(parent(A))'
+
+LazyBandedMatrices.bidiagonaluplo(A::BidiagonalConjugationData) = 'U'
+LazyBandedMatrices.Bidiagonal(A::BidiagonalConjugationData) = LazyBandedMatrices.Bidiagonal(A[band(0)], A[band(1)], 'U')
+
+_colsize(A::BidiagonalConjugationData) = length(A.dv)
+
+function _compute_column!(A::BidiagonalConjugationData, i)
+    # computes A[i, i] and A[i-1, i]
+    i ≤ _colsize(A) && return A
+    dv, ev = A.dv, A.ev
+    U, C = A.U, A.C
+    resize!(dv, i)
+    resize!(ev, i - 1)
+    if i == 1
+        dv[i] = C[1, 1] / U[1, 1]
+    else
+        uᵢ₋₁ᵢ₋₁, uᵢ₋₁ᵢ, uᵢᵢ₋₁, uᵢᵢ = U[i-1, i-1], U[i-1, i], U[i, i-1], U[i, i]
+        cᵢ₋₁ᵢ, cᵢᵢ = C[i-1, i], C[i, i]
+        Uᵢ⁻¹ = inv(uᵢ₋₁ᵢ₋₁ * uᵢᵢ - uᵢ₋₁ᵢ * uᵢᵢ₋₁)
+        dv[i] = Uᵢ⁻¹ * (uᵢ₋₁ᵢ₋₁ * cᵢᵢ - uᵢᵢ₋₁ * cᵢ₋₁ᵢ)
+        ev[i-1] = Uᵢ⁻¹ * (uᵢᵢ * cᵢ₋₁ᵢ - uᵢ₋₁ᵢ * cᵢᵢ)
+    end
+    return A
+end
+
+function getindex(A::BidiagonalConjugationData, i::Int, j::Int)
+    i ≤ 0 || j ≤ 0 && throw(BoundsError(A, (i, j)))
+    T = eltype(A)
+    in_band = i == j || i == j - 1
+    if !in_band
+        return zero(T)
+    elseif j > _colsize(A)
+        _compute_column!(A, j)
+        return i == j ? A.dv[i] : A.ev[i]    
+    else 
+        return i == j ? A.dv[i] : A.ev[i]
+    end
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 2fbaece..8328456 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -476,3 +476,4 @@ include("test_inful.jl")
 include("test_infcholesky.jl")
 include("test_periodic.jl")
 include("test_infreversecholesky.jl")
+include("test_bidiagonalconjugationdata.jl")
\ No newline at end of file
diff --git a/test/test_bidiagonalconjugationdata.jl b/test/test_bidiagonalconjugationdata.jl
new file mode 100644
index 0000000..752a822
--- /dev/null
+++ b/test/test_bidiagonalconjugationdata.jl
@@ -0,0 +1,33 @@
+@testset "BidiagonalConjugationData" begin
+    for _ in 1:10
+        ir = InfiniteArrays.InfRandVector
+        r = () -> Random.seed!(rand(1:2^32)) # avoid https://github.com/JuliaArrays/InfiniteArrays.jl/issues/182
+        V = BandedMatrix(-1 => ir(r()), 0 => ir(r()), 1 => ir(r()))
+        A = BandedMatrix(0 => ir(r()), 1 => ir(r()))
+        X = BandedMatrix(0 => ir(r()), 1 => ir(r()), 2 => ir(r()))
+        U = X * V * inv(A)
+        B = InfiniteLinearAlgebra.BidiagonalConjugationData(U, X, V)
+
+        @test MemoryLayout(B) == BidiagonalLayout{LazyArrays.LazyLayout,LazyArrays.LazyLayout}()
+        @test bandwidths(B) == (0, 1)
+        @test size(B) == (ℵ₀, ℵ₀)
+        @test axes(B) == (OneToInf(), OneToInf())
+        @test eltype(B) == Float64
+        @test copy(B)[1:10, 1:10] == B[1:10, 1:10]
+        @test !(copy(B) === B)
+        @test copy(B')[1:10, 1:10] == B[1:10, 1:10]'
+        @test !(copy(B') === B')
+        @test LazyBandedMatrices.bidiagonaluplo(B) == 'U'
+        @test LazyBandedMatrices.Bidiagonal(B)[1:100, 1:100] == LazyBandedMatrices.Bidiagonal(B[band(0)], B[band(1)], 'U')[1:100, 1:100]
+        @test B[1:100, 1:100] ≈ A[1:100, 1:100]
+        @test B[band(0)][1:1000] ≈ A[band(0)][1:1000]
+        @test B[band(1)][1:1000] ≈ A[band(1)][1:1000]
+        @test (B+B)[1:100, 1:100] ≈ 2(A[1:100, 1:100])
+        @test (B*B)[1:100, 1:100] ≈ (A*A)[1:100, 1:100]
+        @test inv(B)[1:100, 1:100] ≈ inv(A)[1:100, 1:100]
+        @test (B*I)[1:100, 1:100] ≈ B[1:100, 1:100]
+        @test (B*Diagonal(1:∞))[1:100, 1:100] ≈ B[1:100, 1:100] * Diagonal(1:100)
+        @test (U*B)[1:100, 1:100] ≈ (X*V)[1:100, 1:100] rtol=1e-2 atol=1e-4
+        @test (B'B)[1:100, 1:100] ≈ B'[1:100, 1:100] * B[1:100, 1:100]
+    end
+end
\ No newline at end of file

From 3e455e9bb2097a7eeac3de68366e45729035dfd8 Mon Sep 17 00:00:00 2001
From: DanielVandH <danj.vandenheuvel@gmail.com>
Date: Tue, 2 Jul 2024 22:58:30 +0100
Subject: [PATCH 02/13] Redesign

---
 Project.toml                            |   1 +
 src/InfiniteLinearAlgebra.jl            |   2 +-
 src/banded/bidiagonalconjugation.jl     | 133 ++++++++++++++++++++++++
 src/banded/bidiagonalconjugationdata.jl |  85 ---------------
 test/runtests.jl                        |   2 +-
 test/test_bidiagonalconjugation.jl      |  60 +++++++++++
 test/test_bidiagonalconjugationdata.jl  |  33 ------
 7 files changed, 196 insertions(+), 120 deletions(-)
 create mode 100644 src/banded/bidiagonalconjugation.jl
 delete mode 100644 src/banded/bidiagonalconjugationdata.jl
 create mode 100644 test/test_bidiagonalconjugation.jl
 delete mode 100644 test/test_bidiagonalconjugationdata.jl

diff --git a/Project.toml b/Project.toml
index 3dc14a1..1d9edd5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
+InfiniteRandomArrays = "2bc77966-89c7-476d-a40f-269028fac4a9"
 Infinities = "e1ba4f0e-776d-440f-acd9-e1d2e9742647"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
diff --git a/src/InfiniteLinearAlgebra.jl b/src/InfiniteLinearAlgebra.jl
index f75e59e..9ebc0a7 100644
--- a/src/InfiniteLinearAlgebra.jl
+++ b/src/InfiniteLinearAlgebra.jl
@@ -143,6 +143,6 @@ include("infql.jl")
 include("infqr.jl")
 include("inful.jl")
 include("infcholesky.jl")
-include("banded/bidiagonalconjugationdata.jl")
+include("banded/bidiagonalconjugation.jl")
 
 end # module
diff --git a/src/banded/bidiagonalconjugation.jl b/src/banded/bidiagonalconjugation.jl
new file mode 100644
index 0000000..98ea8b9
--- /dev/null
+++ b/src/banded/bidiagonalconjugation.jl
@@ -0,0 +1,133 @@
+"""
+    BidiagonalConjugation{T, MU, MC} <: AbstractCachedMatrix{T}
+
+Struct for efficiently projecting the matrix product `inv(U)XV` onto a 
+bidiagonal matrix.
+
+# Fields 
+- `U`: The matrix `U`.
+- `C`: The matrix product `C = XV`.
+- `data::Bidiagonal{T,Vector{T}}`: The cached data of the projection. 
+- `datasize::Int`: The size of the finite section computed thus far.
+- `dv::BandWrapper`: The diagonal part.
+- `ev::BandWrapper`: The offdiagonal part.
+
+# Constructor
+To construct the projection, use the constructor 
+
+    BidiagonalConjugation(U, X, V, uplo)
+
+where `uplo` specifies whether the projection is upper (`U`) or lower (`L`).
+"""
+mutable struct BidiagonalConjugation{T,MU,MC} <: AbstractCachedMatrix{T}
+    const U::MU
+    const C::MC
+    const data::Bidiagonal{T,Vector{T}}
+    datasize::Int # also compute up to a square finite section 
+    function BidiagonalConjugation(U::MU, C::MC, data::Bidiagonal{T}, datasize::Int) where {T,MU,MC}
+        return new{T,MU,MC}(U, C, data, datasize)
+    end # disambiguate with the next constructor
+end
+function BidiagonalConjugation(U::MU, X, V, uplo) where {MU}
+    C = X * V
+    T = promote_type(typeof(inv(U[1, 1])), eltype(U), eltype(C)) # include inv so that we can't get Ints
+    data = Bidiagonal(T[], T[], uplo)
+    return BidiagonalConjugation(U, C, data, 0)
+end
+get_uplo(A::BidiagonalConjugation) = A.data.uplo
+
+MemoryLayout(::Type{<:BidiagonalConjugation}) = BidiagonalLayout{LazyLayout,LazyLayout}()
+diagonaldata(A::BidiagonalConjugation) = A.dv
+supdiagonaldata(A::BidiagonalConjugation) = get_uplo(A) == 'U' ? A.ev : throw(ArgumentError(LazyString(A, " is lower-bidiagonal")))
+subdiagonaldata(A::BidiagonalConjugation) = get_uplo(A) == 'L' ? A.ev : throw(ArgumentError(LazyString(A, " is upper-bidiagonal")))
+
+bandwidths(A::BidiagonalConjugation) = bandwidths(A.data)
+size(A::BidiagonalConjugation) = (ℵ₀, ℵ₀)
+axes(A::BidiagonalConjugation) = (OneToInf(), OneToInf())
+
+copy(A::BidiagonalConjugation) = BidiagonalConjugation(copy(A.U), copy(A.C), copy(A.data), A.datasize)
+copy(A::Adjoint{T,<:BidiagonalConjugation}) where {T} = copy(parent(A))'
+
+LazyBandedMatrices.bidiagonaluplo(A::BidiagonalConjugation) = get_uplo(A)
+LazyBandedMatrices.Bidiagonal(A::BidiagonalConjugation) = LazyBandedMatrices.Bidiagonal(A.dv, A.ev, get_uplo(A))
+
+# We could decouple this from parent since the computation of each vector is 
+# independent of the other, but it would be slower since they get computed 
+# in tandem. Thus, we instead use a wrapper that captures both. 
+# This is needed so that A.dv and A.ev can be defined (and are useful, instead of 
+# just returning A.data.dv and A.data.ev which are finite and have no knowledge 
+# of the parent).
+struct BandWrapper{T,P<:BidiagonalConjugation{T}} <: LazyVector{T}
+    parent::P
+    diag::Bool # true => main diagonal, false => off diagonal 
+end
+size(wrap::BandWrapper) = (size(wrap.parent, 1),)
+function getindex(wrap::BandWrapper, i::Int)
+    parent = wrap.parent
+    uplo = get_uplo(parent)
+    if wrap.diag
+        return parent[i, i]
+    elseif uplo == 'U'
+        return parent[i, i+1]
+    else # uplo == 'L' 
+        return parent[i+1, i]
+    end
+end
+
+function getproperty(A::BidiagonalConjugation, d::Symbol)
+    if d == :dv
+        return BandWrapper(A, true)
+    elseif d == :ev
+        return BandWrapper(A, false)
+    else
+        return getfield(A, d)
+    end
+end
+
+function _compute_column_up!(A::BidiagonalConjugation, i)
+    U, C = A.U, A.C
+    data = A.data
+    dv, ev = data.dv, data.ev
+    if i == 1
+        dv[i] = C[1, 1] / U[1, 1]
+    else
+        uᵢ₋₁ᵢ₋₁, uᵢᵢ₋₁, uᵢ₋₁ᵢ, uᵢᵢ = U[i-1, i-1], U[i, i-1], U[i-1, i], U[i, i]
+        cᵢ₋₁ᵢ, cᵢᵢ = C[i-1, i], C[i, i]
+        Uᵢ⁻¹ = inv(uᵢ₋₁ᵢ₋₁ * uᵢᵢ - uᵢ₋₁ᵢ * uᵢᵢ₋₁)
+        dv[i] = Uᵢ⁻¹ * (uᵢ₋₁ᵢ₋₁ * cᵢᵢ - uᵢᵢ₋₁ * cᵢ₋₁ᵢ)
+        ev[i-1] = Uᵢ⁻¹ * (uᵢᵢ * cᵢ₋₁ᵢ - uᵢ₋₁ᵢ * cᵢᵢ)
+    end
+    return A
+end
+
+function _compute_column_lo!(A::BidiagonalConjugation, i)
+    U, C = A.U, A.C
+    data = A.data
+    dv, ev = data.dv, data.ev
+    uᵢᵢ, uᵢ₊₁ᵢ, uᵢᵢ₊₁, uᵢ₊₁ᵢ₊₁ = U[i, i], U[i+1, i], U[i, i+1], U[i+1, i+1]
+    cᵢᵢ, cᵢ₊₁ᵢ = C[i, i], C[i+1, i]
+    Uᵢ⁻¹ = inv(uᵢᵢ * uᵢ₊₁ᵢ₊₁ - uᵢᵢ₊₁ * uᵢ₊₁ᵢ)
+    dv[i] = Uᵢ⁻¹ * (uᵢ₊₁ᵢ₊₁ * cᵢᵢ - uᵢᵢ₊₁ * cᵢ₊₁ᵢ)
+    ev[i] = Uᵢ⁻¹ * (uᵢᵢ * cᵢ₊₁ᵢ - uᵢ₊₁ᵢ * cᵢᵢ)
+    return A
+end
+
+function _compute_columns!(A::BidiagonalConjugation, i)
+    ds = A.datasize
+    up = get_uplo(A) == 'U'
+    for j in (ds+1):i
+        up ? _compute_column_up!(A, j) : _compute_column_lo!(A, j)
+    end
+    A.datasize = i
+    return A
+end
+
+function resizedata!(A::BidiagonalConjugation, i::Int, j::Int)
+    ds = A.datasize
+    j ≤ ds && return A
+    data = A.data
+    dv, ev = data.dv, data.ev
+    resize!(dv, 2j + 1)
+    resize!(ev, 2j)
+    return _compute_columns!(A, 2j)
+end
\ No newline at end of file
diff --git a/src/banded/bidiagonalconjugationdata.jl b/src/banded/bidiagonalconjugationdata.jl
deleted file mode 100644
index ed8a3d8..0000000
--- a/src/banded/bidiagonalconjugationdata.jl
+++ /dev/null
@@ -1,85 +0,0 @@
-"""
-    BidiagonalConjugationData{T, MU, MC} <: LazyMatrix{T}
-
-Struct for efficiently representing the matrix product `A = inv(U)XV`,
-assuming that 
-
-- `A` is upper bidiagonal, 
-- `U` is upper Hessenberg,
-- `X` is banded,
-- `V` is banded.
-
-None of these properties are checked internally. It is the user's responsibility
-to ensure these properties hold and that the product `inv(U)XV` is indeed bidiagonal.
-
-# Fields 
-- `U`: The upper Hessenberg matrix. 
-- `C`: The matrix product `XV`.
-- `dv`: A vector giving the diagonal of `A`. 
-- `ev`: A vector giving the superdiagonal of `A`.
-
-The vectors `dv` and `ev` grow on demand from `getindex` and should not be 
-used directly. Simply treat 
-
-    A = BidiagonalConjugationData(U, X, V)
-
-as you would an upper bidiagonal matrix.
-"""
-struct BidiagonalConjugationData{T,MU,MC} <: LazyMatrix{T}
-    U::MU
-    C::MC
-    dv::Vector{T}
-    ev::Vector{T}
-end
-function BidiagonalConjugationData(U::MU, X::MX, V::MV) where {MU,MX,MV}
-    C = X * V
-    T = promote_type(typeof(inv(U[begin])), eltype(U), eltype(C)) # include inv so that we can't get Ints
-    dv, ev = T[], T[]
-    return BidiagonalConjugationData{T,MU,typeof(C)}(U, C, dv, ev)
-end
-MemoryLayout(::Type{<:BidiagonalConjugationData}) = BidiagonalLayout{LazyLayout,LazyLayout}()
-bandwidths(A::BidiagonalConjugationData) = (0, 1)
-size(A::BidiagonalConjugationData) = (ℵ₀, ℵ₀)
-axes(A::BidiagonalConjugationData) = (OneToInf(), OneToInf())
-Base.eltype(A::Type{<:BidiagonalConjugationData{T}}) where {T} = T
-
-copy(A::BidiagonalConjugationData) = BidiagonalConjugationData(copy(A.U), copy(A.C), copy(A.dv), copy(A.ev))
-copy(A::Adjoint{T,<:BidiagonalConjugationData}) where {T} = copy(parent(A))'
-
-LazyBandedMatrices.bidiagonaluplo(A::BidiagonalConjugationData) = 'U'
-LazyBandedMatrices.Bidiagonal(A::BidiagonalConjugationData) = LazyBandedMatrices.Bidiagonal(A[band(0)], A[band(1)], 'U')
-
-_colsize(A::BidiagonalConjugationData) = length(A.dv)
-
-function _compute_column!(A::BidiagonalConjugationData, i)
-    # computes A[i, i] and A[i-1, i]
-    i ≤ _colsize(A) && return A
-    dv, ev = A.dv, A.ev
-    U, C = A.U, A.C
-    resize!(dv, i)
-    resize!(ev, i - 1)
-    if i == 1
-        dv[i] = C[1, 1] / U[1, 1]
-    else
-        uᵢ₋₁ᵢ₋₁, uᵢ₋₁ᵢ, uᵢᵢ₋₁, uᵢᵢ = U[i-1, i-1], U[i-1, i], U[i, i-1], U[i, i]
-        cᵢ₋₁ᵢ, cᵢᵢ = C[i-1, i], C[i, i]
-        Uᵢ⁻¹ = inv(uᵢ₋₁ᵢ₋₁ * uᵢᵢ - uᵢ₋₁ᵢ * uᵢᵢ₋₁)
-        dv[i] = Uᵢ⁻¹ * (uᵢ₋₁ᵢ₋₁ * cᵢᵢ - uᵢᵢ₋₁ * cᵢ₋₁ᵢ)
-        ev[i-1] = Uᵢ⁻¹ * (uᵢᵢ * cᵢ₋₁ᵢ - uᵢ₋₁ᵢ * cᵢᵢ)
-    end
-    return A
-end
-
-function getindex(A::BidiagonalConjugationData, i::Int, j::Int)
-    i ≤ 0 || j ≤ 0 && throw(BoundsError(A, (i, j)))
-    T = eltype(A)
-    in_band = i == j || i == j - 1
-    if !in_band
-        return zero(T)
-    elseif j > _colsize(A)
-        _compute_column!(A, j)
-        return i == j ? A.dv[i] : A.ev[i]    
-    else 
-        return i == j ? A.dv[i] : A.ev[i]
-    end
-end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 8328456..cc6ce8c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -476,4 +476,4 @@ include("test_inful.jl")
 include("test_infcholesky.jl")
 include("test_periodic.jl")
 include("test_infreversecholesky.jl")
-include("test_bidiagonalconjugationdata.jl")
\ No newline at end of file
+include("test_bidiagonalconjugation.jl")
\ No newline at end of file
diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl
new file mode 100644
index 0000000..de7553d
--- /dev/null
+++ b/test/test_bidiagonalconjugation.jl
@@ -0,0 +1,60 @@
+@testset "BidiagonalConjugationData" begin
+    for _ in 1:10
+        V1 = InfRandTridiagonal()
+        A1 = InfRandBidiagonal('U')
+        X1 = brand(∞, 0, 2)
+        U1 = X1 * V1 * ApplyArray(inv, A1)
+        B1 = BidiagonalConjugation(U1, X1, V1, 'U')
+
+        V2 = brand(∞, 0, 1)
+        A2 = LazyBandedMatrices.Bidiagonal(Fill(0.2, ∞), 2.0 ./ (1.0 .+ (1:∞)), 'L') # LinearAlgebra.Bidiagonal not playing nice for this case
+        X2 = InfRandBidiagonal('L')
+        U2 = X2 * V2 * ApplyArray(inv, A2) # U2 isn't upper Hessenberg (actually it's lower, oh well), doesn't seem to matter for computing A
+        B2 = BidiagonalConjugation(U2, X2, V2, 'L')
+
+        for (A, B, uplo) in ((A1, B1, 'U'), (A2, B2, 'L'))
+            @test get_uplo(B) == uplo
+            @test MemoryLayout(B) isa BidiagonalLayout{LazyLayout,LazyLayout}
+            @test diagonaldata(B) === B.dv
+            if uplo == 'U'
+                @test supdiagonaldata(B) === B.ev
+                @test_throws ArgumentError subdiagonaldata(B)
+                @test bandwidths(B) == (0, 1)
+            else
+                @test subdiagonaldata(B) === B.ev
+                @test_throws ArgumentError supdiagonaldata(B)
+                @test bandwidths(B) == (1, 0)
+            end
+            @test size(B) == (ℵ₀, ℵ₀)
+            @test axes(B) == (OneToInf(), OneToInf())
+            @test eltype(B) == Float64
+            for _B in (B, B')
+                BB = copy(_B)
+                @test parent(BB).datasize == parent(_B).datasize
+                @test !(BB === B) && !(parent(BB).data === parent(B).data)
+                @test BB[1:10, 1:10] == _B[1:10, 1:10]
+            end
+            @test LazyBandedMatrices.bidiagonaluplo(B) == uplo
+            @test LazyBandedMatrices.Bidiagonal(B) === LazyBandedMatrices.Bidiagonal(B.dv, B.ev, Symbol(uplo))
+            @test B[1:10, 1:10] ≈ A[1:10, 1:10]
+            @test B[230, 230] ≈ A[230, 230]
+            @test B[102, 102] ≈ A[102, 102] # make sure we compute intermediate columns correctly when skipping 
+            @test B[band(0)][1:100] == B.dv[1:100]
+            if uplo == 'U'
+                @test B[band(1)][1:100] == B.ev[1:100]
+                @test B[band(-1)][1:100] == zeros(100)
+            else
+                @test B[band(-1)][1:100] == B.ev[1:100]
+                @test B[band(1)][1:100] == zeros(100)
+            end
+            @test B.dv[500] == B.data.dv[500]
+            @test B.datasize == 1000 # make sure wrapper is working correctly 
+            @test B.ev[1005] == B.data.ev[1005]
+            @test B.datasize == 2010 + 2 * (uplo == 'U')
+            # @test_broken ApplyArray(inv, B)[1:100, 1:100] ≈ ApplyArray(inv, A)[1:100, 1:100] # need to somehow let inv (or even ApplyArray(inv, )) work
+            @test (B+B)[1:100, 1:100] ≈ 2A[1:100, 1:100] ≈ 2B[1:100, 1:100]
+            @test (B*I)[1:100, 1:100] ≈ B[1:100, 1:100]
+            @test (B*Diagonal(1:∞))[1:100, 1:100] ≈ B[1:100, 1:100] * Diagonal(1:100)
+        end
+    end
+end
\ No newline at end of file
diff --git a/test/test_bidiagonalconjugationdata.jl b/test/test_bidiagonalconjugationdata.jl
deleted file mode 100644
index 752a822..0000000
--- a/test/test_bidiagonalconjugationdata.jl
+++ /dev/null
@@ -1,33 +0,0 @@
-@testset "BidiagonalConjugationData" begin
-    for _ in 1:10
-        ir = InfiniteArrays.InfRandVector
-        r = () -> Random.seed!(rand(1:2^32)) # avoid https://github.com/JuliaArrays/InfiniteArrays.jl/issues/182
-        V = BandedMatrix(-1 => ir(r()), 0 => ir(r()), 1 => ir(r()))
-        A = BandedMatrix(0 => ir(r()), 1 => ir(r()))
-        X = BandedMatrix(0 => ir(r()), 1 => ir(r()), 2 => ir(r()))
-        U = X * V * inv(A)
-        B = InfiniteLinearAlgebra.BidiagonalConjugationData(U, X, V)
-
-        @test MemoryLayout(B) == BidiagonalLayout{LazyArrays.LazyLayout,LazyArrays.LazyLayout}()
-        @test bandwidths(B) == (0, 1)
-        @test size(B) == (ℵ₀, ℵ₀)
-        @test axes(B) == (OneToInf(), OneToInf())
-        @test eltype(B) == Float64
-        @test copy(B)[1:10, 1:10] == B[1:10, 1:10]
-        @test !(copy(B) === B)
-        @test copy(B')[1:10, 1:10] == B[1:10, 1:10]'
-        @test !(copy(B') === B')
-        @test LazyBandedMatrices.bidiagonaluplo(B) == 'U'
-        @test LazyBandedMatrices.Bidiagonal(B)[1:100, 1:100] == LazyBandedMatrices.Bidiagonal(B[band(0)], B[band(1)], 'U')[1:100, 1:100]
-        @test B[1:100, 1:100] ≈ A[1:100, 1:100]
-        @test B[band(0)][1:1000] ≈ A[band(0)][1:1000]
-        @test B[band(1)][1:1000] ≈ A[band(1)][1:1000]
-        @test (B+B)[1:100, 1:100] ≈ 2(A[1:100, 1:100])
-        @test (B*B)[1:100, 1:100] ≈ (A*A)[1:100, 1:100]
-        @test inv(B)[1:100, 1:100] ≈ inv(A)[1:100, 1:100]
-        @test (B*I)[1:100, 1:100] ≈ B[1:100, 1:100]
-        @test (B*Diagonal(1:∞))[1:100, 1:100] ≈ B[1:100, 1:100] * Diagonal(1:100)
-        @test (U*B)[1:100, 1:100] ≈ (X*V)[1:100, 1:100] rtol=1e-2 atol=1e-4
-        @test (B'B)[1:100, 1:100] ≈ B'[1:100, 1:100] * B[1:100, 1:100]
-    end
-end
\ No newline at end of file

From e79939a63a7401a81ac4bdd11b6f8ec7f3e786ae Mon Sep 17 00:00:00 2001
From: DanielVandH <danj.vandenheuvel@gmail.com>
Date: Tue, 2 Jul 2024 23:00:18 +0100
Subject: [PATCH 03/13] Fix import and make infiniterandomarrays an extra

---
 Project.toml                 | 5 +++--
 src/InfiniteLinearAlgebra.jl | 2 +-
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/Project.toml b/Project.toml
index 1d9edd5..7d88418 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,7 +9,6 @@ BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
-InfiniteRandomArrays = "2bc77966-89c7-476d-a40f-269028fac4a9"
 Infinities = "e1ba4f0e-776d-440f-acd9-e1d2e9742647"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
@@ -25,6 +24,7 @@ BlockArrays = "1.0"
 BlockBandedMatrices = "0.13"
 FillArrays = "1.0"
 InfiniteArrays = "0.14"
+InfiniteRandomArrays = "0.2"
 Infinities = "0.1"
 LazyArrays = "2.0"
 LazyBandedMatrices = "0.10"
@@ -40,9 +40,10 @@ julia = "1.10"
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+InfiniteRandomArrays = "2bc77966-89c7-476d-a40f-269028fac4a9"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Test", "Random", "SpecialFunctions", "StaticArrays"]
+test = ["Aqua", "Test", "Random", "InfiniteRandomArrays", "SpecialFunctions", "StaticArrays"]
diff --git a/src/InfiniteLinearAlgebra.jl b/src/InfiniteLinearAlgebra.jl
index 9ebc0a7..50bdace 100644
--- a/src/InfiniteLinearAlgebra.jl
+++ b/src/InfiniteLinearAlgebra.jl
@@ -38,7 +38,7 @@ import Infinities: InfiniteCardinal, Infinity
 
 import LazyArrays: AbstractCachedMatrix, AbstractCachedVector, AbstractLazyLayout, ApplyArray, ApplyLayout, ApplyMatrix,
                    CachedArray, CachedLayout, CachedMatrix, CachedVector, LazyArrayStyle, LazyLayout,
-                   LazyLayouts, LazyMatrix, AbstractPaddedLayout, PaddedColumns, _broadcast_sub_arguments,
+                   LazyLayouts, LazyMatrix, LazyVector, AbstractPaddedLayout, PaddedColumns, _broadcast_sub_arguments,
                    applybroadcaststyle, applylayout, arguments, cacheddata, paddeddata, resizedata!, simplifiable,
                    simplify, islazy, islazy_layout
 

From e7c37bde87bcb36b48796b00b20577be7f9aed04 Mon Sep 17 00:00:00 2001
From: DanielVandH <danj.vandenheuvel@gmail.com>
Date: Tue, 2 Jul 2024 23:02:57 +0100
Subject: [PATCH 04/13] typo

---
 src/banded/bidiagonalconjugation.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/banded/bidiagonalconjugation.jl b/src/banded/bidiagonalconjugation.jl
index 98ea8b9..991b107 100644
--- a/src/banded/bidiagonalconjugation.jl
+++ b/src/banded/bidiagonalconjugation.jl
@@ -23,7 +23,7 @@ mutable struct BidiagonalConjugation{T,MU,MC} <: AbstractCachedMatrix{T}
     const U::MU
     const C::MC
     const data::Bidiagonal{T,Vector{T}}
-    datasize::Int # also compute up to a square finite section 
+    datasize::Int # always compute up to a square finite section 
     function BidiagonalConjugation(U::MU, C::MC, data::Bidiagonal{T}, datasize::Int) where {T,MU,MC}
         return new{T,MU,MC}(U, C, data, datasize)
     end # disambiguate with the next constructor

From e2cb23180f51b5787d24560d4f6c59580f60269e Mon Sep 17 00:00:00 2001
From: DanielVandH <danj.vandenheuvel@gmail.com>
Date: Thu, 4 Jul 2024 15:20:21 +0100
Subject: [PATCH 05/13] Tests

---
 test/runtests.jl | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index cc6ce8c..f4695c7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,15 +3,18 @@ using InfiniteLinearAlgebra, BlockBandedMatrices, BlockArrays, BandedMatrices, I
 import InfiniteLinearAlgebra: qltail, toeptail, tailiterate, tailiterate!, tail_de, ql_X!,
     InfToeplitz, PertToeplitz, TriToeplitz, InfBandedMatrix, InfBandCartesianIndices,
     rightasymptotics, QLHessenberg, ConstRows, PertConstRows, chop, chop!, pad,
-    BandedToeplitzLayout, PertToeplitzLayout, TridiagonalToeplitzLayout, BidiagonalToeplitzLayout
+    BandedToeplitzLayout, PertToeplitzLayout, TridiagonalToeplitzLayout, BidiagonalToeplitzLayout,
+    BidiagonalConjugation, get_uplo
 import Base: BroadcastStyle, oneto
 import BlockArrays: _BlockArray, blockcolsupport
 import BlockBandedMatrices: isblockbanded, _BlockBandedMatrix
 import MatrixFactorizations: QLPackedQ
 import BandedMatrices: bandeddata, _BandedMatrix, BandedStyle
-import LazyArrays: colsupport, MemoryLayout, ApplyLayout, LazyArrayStyle, arguments, paddeddata, PaddedColumns
+import LazyArrays: colsupport, MemoryLayout, ApplyLayout, LazyArrayStyle, arguments, paddeddata, PaddedColumns, LazyLayout
 import InfiniteArrays: OneToInf, oneto, RealInfinity
 import LazyBandedMatrices: BroadcastBandedBlockBandedLayout, BroadcastBandedLayout, LazyBandedLayout, BlockVec
+import InfiniteRandomArrays: InfRandTridiagonal, InfRandBidiagonal
+import ArrayLayouts: diagonaldata, supdiagonaldata, subdiagonaldata
 
 using Aqua
 @testset "Project quality" begin

From 563e0b6a21ab5c2368f54bab28c64d1d388833e1 Mon Sep 17 00:00:00 2001
From: DanielVandH <danj.vandenheuvel@gmail.com>
Date: Sat, 6 Jul 2024 13:19:23 +0100
Subject: [PATCH 06/13] Redesign

---
 src/InfiniteLinearAlgebra.jl        |   2 +-
 src/banded/bidiagonalconjugation.jl | 210 +++++++++++++++-------------
 test/runtests.jl                    |   2 +-
 test/test_bidiagonalconjugation.jl  |  40 ++++--
 4 files changed, 138 insertions(+), 116 deletions(-)

diff --git a/src/InfiniteLinearAlgebra.jl b/src/InfiniteLinearAlgebra.jl
index 50bdace..175f1c2 100644
--- a/src/InfiniteLinearAlgebra.jl
+++ b/src/InfiniteLinearAlgebra.jl
@@ -40,7 +40,7 @@ import LazyArrays: AbstractCachedMatrix, AbstractCachedVector, AbstractLazyLayou
                    CachedArray, CachedLayout, CachedMatrix, CachedVector, LazyArrayStyle, LazyLayout,
                    LazyLayouts, LazyMatrix, LazyVector, AbstractPaddedLayout, PaddedColumns, _broadcast_sub_arguments,
                    applybroadcaststyle, applylayout, arguments, cacheddata, paddeddata, resizedata!, simplifiable,
-                   simplify, islazy, islazy_layout
+                   simplify, islazy, islazy_layout, cache_getindex
 
 import LazyBandedMatrices: AbstractLazyBandedBlockBandedLayout, AbstractLazyBandedLayout, ApplyBandedLayout, BlockVec,
                            BroadcastBandedLayout, KronTravBandedBlockBandedLayout, LazyBandedLayout,
diff --git a/src/banded/bidiagonalconjugation.jl b/src/banded/bidiagonalconjugation.jl
index 991b107..d68b51f 100644
--- a/src/banded/bidiagonalconjugation.jl
+++ b/src/banded/bidiagonalconjugation.jl
@@ -1,92 +1,43 @@
-"""
-    BidiagonalConjugation{T, MU, MC} <: AbstractCachedMatrix{T}
-
-Struct for efficiently projecting the matrix product `inv(U)XV` onto a 
-bidiagonal matrix.
-
-# Fields 
-- `U`: The matrix `U`.
-- `C`: The matrix product `C = XV`.
-- `data::Bidiagonal{T,Vector{T}}`: The cached data of the projection. 
-- `datasize::Int`: The size of the finite section computed thus far.
-- `dv::BandWrapper`: The diagonal part.
-- `ev::BandWrapper`: The offdiagonal part.
-
-# Constructor
-To construct the projection, use the constructor 
-
-    BidiagonalConjugation(U, X, V, uplo)
+@inline function _to_uplo(char::Symbol)
+    if char == :U
+        'U'
+    elseif char == :L
+        'L'
+    else
+        _throw_uplo()
+    end
+end
+@inline function _to_uplo(char::Char)
+    if char ∈ ('L', 'U')
+        char
+    else
+        _throw_uplo()
+    end
+end
+@noinline _throw_uplo() = throw(ArgumentError("uplo argument must be either :U (upper) or :L (lower)"))
 
-where `uplo` specifies whether the projection is upper (`U`) or lower (`L`).
-"""
-mutable struct BidiagonalConjugation{T,MU,MC} <: AbstractCachedMatrix{T}
-    const U::MU
-    const C::MC
-    const data::Bidiagonal{T,Vector{T}}
-    datasize::Int # always compute up to a square finite section 
-    function BidiagonalConjugation(U::MU, C::MC, data::Bidiagonal{T}, datasize::Int) where {T,MU,MC}
-        return new{T,MU,MC}(U, C, data, datasize)
-    end # disambiguate with the next constructor
+mutable struct BidiagonalConjugationData{T}
+    const U::AbstractMatrix{T} # Typing these concretely prevents the use of Bidiagonal, unless we want LazyBandedMatrices.Bidiagonal
+    const C::AbstractMatrix{T} # Function barriers help to minimise the penalty from this when resizing anyway.
+    const dv::Vector{T}
+    const ev::Vector{T}
+    const uplo::Char
+    datasize::Int # Number of columns
 end
-function BidiagonalConjugation(U::MU, X, V, uplo) where {MU}
+function BidiagonalConjugationData(U, X, V, uplo::Char)
     C = X * V
     T = promote_type(typeof(inv(U[1, 1])), eltype(U), eltype(C)) # include inv so that we can't get Ints
-    data = Bidiagonal(T[], T[], uplo)
-    return BidiagonalConjugation(U, C, data, 0)
-end
-get_uplo(A::BidiagonalConjugation) = A.data.uplo
-
-MemoryLayout(::Type{<:BidiagonalConjugation}) = BidiagonalLayout{LazyLayout,LazyLayout}()
-diagonaldata(A::BidiagonalConjugation) = A.dv
-supdiagonaldata(A::BidiagonalConjugation) = get_uplo(A) == 'U' ? A.ev : throw(ArgumentError(LazyString(A, " is lower-bidiagonal")))
-subdiagonaldata(A::BidiagonalConjugation) = get_uplo(A) == 'L' ? A.ev : throw(ArgumentError(LazyString(A, " is upper-bidiagonal")))
-
-bandwidths(A::BidiagonalConjugation) = bandwidths(A.data)
-size(A::BidiagonalConjugation) = (ℵ₀, ℵ₀)
-axes(A::BidiagonalConjugation) = (OneToInf(), OneToInf())
-
-copy(A::BidiagonalConjugation) = BidiagonalConjugation(copy(A.U), copy(A.C), copy(A.data), A.datasize)
-copy(A::Adjoint{T,<:BidiagonalConjugation}) where {T} = copy(parent(A))'
-
-LazyBandedMatrices.bidiagonaluplo(A::BidiagonalConjugation) = get_uplo(A)
-LazyBandedMatrices.Bidiagonal(A::BidiagonalConjugation) = LazyBandedMatrices.Bidiagonal(A.dv, A.ev, get_uplo(A))
-
-# We could decouple this from parent since the computation of each vector is 
-# independent of the other, but it would be slower since they get computed 
-# in tandem. Thus, we instead use a wrapper that captures both. 
-# This is needed so that A.dv and A.ev can be defined (and are useful, instead of 
-# just returning A.data.dv and A.data.ev which are finite and have no knowledge 
-# of the parent).
-struct BandWrapper{T,P<:BidiagonalConjugation{T}} <: LazyVector{T}
-    parent::P
-    diag::Bool # true => main diagonal, false => off diagonal 
-end
-size(wrap::BandWrapper) = (size(wrap.parent, 1),)
-function getindex(wrap::BandWrapper, i::Int)
-    parent = wrap.parent
-    uplo = get_uplo(parent)
-    if wrap.diag
-        return parent[i, i]
-    elseif uplo == 'U'
-        return parent[i, i+1]
-    else # uplo == 'L' 
-        return parent[i+1, i]
-    end
+    dv, ev = T[], T[]
+    return BidiagonalConjugationData(U, C, dv, ev, uplo, 0)
 end
 
-function getproperty(A::BidiagonalConjugation, d::Symbol)
-    if d == :dv
-        return BandWrapper(A, true)
-    elseif d == :ev
-        return BandWrapper(A, false)
-    else
-        return getfield(A, d)
-    end
+function copy(data::BidiagonalConjugationData)
+    U, C, dv, ev, uplo, datasize = data.U, data.C, data.dv, data.ev, data.uplo, data.datasize
+    return BidiagonalConjugationData(copy(U), copy(C), copy(dv), copy(ev), uplo, datasize)
 end
 
-function _compute_column_up!(A::BidiagonalConjugation, i)
-    U, C = A.U, A.C
-    data = A.data
+function _compute_column_up!(data::BidiagonalConjugationData, i)
+    U, C = data.U, data.C
     dv, ev = data.dv, data.ev
     if i == 1
         dv[i] = C[1, 1] / U[1, 1]
@@ -97,37 +48,96 @@ function _compute_column_up!(A::BidiagonalConjugation, i)
         dv[i] = Uᵢ⁻¹ * (uᵢ₋₁ᵢ₋₁ * cᵢᵢ - uᵢᵢ₋₁ * cᵢ₋₁ᵢ)
         ev[i-1] = Uᵢ⁻¹ * (uᵢᵢ * cᵢ₋₁ᵢ - uᵢ₋₁ᵢ * cᵢᵢ)
     end
-    return A
+    return data
 end
 
-function _compute_column_lo!(A::BidiagonalConjugation, i)
-    U, C = A.U, A.C
-    data = A.data
+function _compute_column_lo!(data::BidiagonalConjugationData, i)
+    U, C = data.U, data.C
     dv, ev = data.dv, data.ev
     uᵢᵢ, uᵢ₊₁ᵢ, uᵢᵢ₊₁, uᵢ₊₁ᵢ₊₁ = U[i, i], U[i+1, i], U[i, i+1], U[i+1, i+1]
     cᵢᵢ, cᵢ₊₁ᵢ = C[i, i], C[i+1, i]
     Uᵢ⁻¹ = inv(uᵢᵢ * uᵢ₊₁ᵢ₊₁ - uᵢᵢ₊₁ * uᵢ₊₁ᵢ)
     dv[i] = Uᵢ⁻¹ * (uᵢ₊₁ᵢ₊₁ * cᵢᵢ - uᵢᵢ₊₁ * cᵢ₊₁ᵢ)
     ev[i] = Uᵢ⁻¹ * (uᵢᵢ * cᵢ₊₁ᵢ - uᵢ₊₁ᵢ * cᵢᵢ)
-    return A
+    return data
 end
 
-function _compute_columns!(A::BidiagonalConjugation, i)
-    ds = A.datasize
-    up = get_uplo(A) == 'U'
+function _compute_columns!(data::BidiagonalConjugationData, i)
+    ds = data.datasize
+    up = data.uplo == 'U'
     for j in (ds+1):i
-        up ? _compute_column_up!(A, j) : _compute_column_lo!(A, j)
+        up ? _compute_column_up!(data, j) : _compute_column_lo!(data, j)
     end
-    A.datasize = i
-    return A
+    data.datasize = i
+    return data
 end
 
-function resizedata!(A::BidiagonalConjugation, i::Int, j::Int)
-    ds = A.datasize
-    j ≤ ds && return A
-    data = A.data
+function resizedata!(data::BidiagonalConjugationData, n)
+    ds = data.datasize
+    n ≤ ds && return data
     dv, ev = data.dv, data.ev
-    resize!(dv, 2j + 1)
-    resize!(ev, 2j)
-    return _compute_columns!(A, 2j)
-end
\ No newline at end of file
+    resize!(dv, 2n + 1)
+    resize!(ev, 2n)
+    return _compute_columns!(data, 2n)
+end
+
+struct BidiagonalConjugationBand{T} <: AbstractCachedVector{T}
+    data::BidiagonalConjugationData{T}
+    diag::Bool # true => diagonal, false => offdiagonal 
+end
+@inline size(A::BidiagonalConjugationBand) = (ℵ₀,)
+@inline resizedata!(A::BidiagonalConjugationBand, n) = resizedata!(A.data, n)
+
+function _bcb_getindex(band::BidiagonalConjugationBand, I)
+    resizedata!(band, maximum(I) + 1)
+    if band.diag 
+        return band.data.dv[I]
+    else 
+        return band.data.ev[I]
+    end
+end
+
+@inline cache_getindex(band::BidiagonalConjugationBand, I::Integer) = _bcb_getindex(band, I) # needed for show
+@inline getindex(band::BidiagonalConjugationBand, I::Integer) = cache_getindex(band, I) # also needed for show because of isassigned
+@inline getindex(band::BidiagonalConjugationBand, I::AbstractVector) = _bcb_getindex(band, I)
+#@inline getindex(band::BidiagonalConjugationBand, I::AbstractInfUnitRange{<:Integer}) = view(band, I)
+#@inline getindex(band::SubArray{<:Any, 1, <:BidiagonalConjugationBand}, I::AbstractInfUnitRange{<:Integer}) = view(band, I)
+
+copy(band::BidiagonalConjugationBand) = band
+
+function convert(::Type{BidiagonalConjugationBand{T}}, ::Zeros{V, 1, Tuple{OneToInf{Int}}}) where {T, V}
+    # without this method, we can't do e.g. B[Band(-1)] for upper bidiagonal
+    
+end
+
+const BidiagonalConjugation{T} = Bidiagonal{T, BidiagonalConjugationBand{T}}
+
+"""
+    BidiagonalConjugation(U, X, V, uplo)
+
+Efficiently compute the projection of the matrix product
+`inv(U)XV` onto a bidiagonal matrix. The `uplo` argument 
+specifies whether the projection is upper (`uplo = 'U'`)
+or lower (`uplo = 'L'`) bidiagonal. 
+
+The computation is returned as a `Bidiagonal` matrix whose 
+diagonal and off-diagonal vectors are computed lazily.
+"""
+function BidiagonalConjugation(U, X, V, uplo)
+    _uplo = _to_uplo(uplo)
+    data = BidiagonalConjugationData(U, X, V, _uplo)
+    return _BidiagonalConjugation(data, _uplo)
+end
+
+function _BidiagonalConjugation(data, uplo) # need uplo argument so that we can take transposes
+    dv = BidiagonalConjugationBand(data, true)
+    ev = BidiagonalConjugationBand(data, false)
+    return Bidiagonal(dv, ev, uplo)
+end
+
+function copy(A::BidiagonalConjugation) # need a new method so that the data remains aliased between dv and ev 
+    data = copy(A.dv.data)
+    return _BidiagonalConjugation(data, A.uplo)
+end
+
+LazyBandedMatrices.Bidiagonal(A::BidiagonalConjugation) = LazyBandedMatrices.Bidiagonal(A.dv, A.ev, A.uplo)
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index f4695c7..4481369 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,7 +4,7 @@ import InfiniteLinearAlgebra: qltail, toeptail, tailiterate, tailiterate!, tail_
     InfToeplitz, PertToeplitz, TriToeplitz, InfBandedMatrix, InfBandCartesianIndices,
     rightasymptotics, QLHessenberg, ConstRows, PertConstRows, chop, chop!, pad,
     BandedToeplitzLayout, PertToeplitzLayout, TridiagonalToeplitzLayout, BidiagonalToeplitzLayout,
-    BidiagonalConjugation, get_uplo
+    BidiagonalConjugation
 import Base: BroadcastStyle, oneto
 import BlockArrays: _BlockArray, blockcolsupport
 import BlockBandedMatrices: isblockbanded, _BlockBandedMatrix
diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl
index de7553d..79db456 100644
--- a/test/test_bidiagonalconjugation.jl
+++ b/test/test_bidiagonalconjugation.jl
@@ -1,4 +1,11 @@
 @testset "BidiagonalConjugationData" begin
+    @test InfiniteLinearAlgebra._to_uplo('U') == 'U'
+    @test InfiniteLinearAlgebra._to_uplo('L') == 'L'
+    @test_throws ArgumentError InfiniteLinearAlgebra._to_uplo('a')
+    @test InfiniteLinearAlgebra._to_uplo(:U) == 'U'
+    @test InfiniteLinearAlgebra._to_uplo(:L) == 'L'
+    @test_throws ArgumentError InfiniteLinearAlgebra._to_uplo(:a)
+
     for _ in 1:10
         V1 = InfRandTridiagonal()
         A1 = InfRandBidiagonal('U')
@@ -10,10 +17,10 @@
         A2 = LazyBandedMatrices.Bidiagonal(Fill(0.2, ∞), 2.0 ./ (1.0 .+ (1:∞)), 'L') # LinearAlgebra.Bidiagonal not playing nice for this case
         X2 = InfRandBidiagonal('L')
         U2 = X2 * V2 * ApplyArray(inv, A2) # U2 isn't upper Hessenberg (actually it's lower, oh well), doesn't seem to matter for computing A
-        B2 = BidiagonalConjugation(U2, X2, V2, 'L')
+        B2 = BidiagonalConjugation(U2, X2, V2, :L)
 
         for (A, B, uplo) in ((A1, B1, 'U'), (A2, B2, 'L'))
-            @test get_uplo(B) == uplo
+            @test B.dv.data === B.ev.data
             @test MemoryLayout(B) isa BidiagonalLayout{LazyLayout,LazyLayout}
             @test diagonaldata(B) === B.dv
             if uplo == 'U'
@@ -30,9 +37,12 @@
             @test eltype(B) == Float64
             for _B in (B, B')
                 BB = copy(_B)
-                @test parent(BB).datasize == parent(_B).datasize
-                @test !(BB === B) && !(parent(BB).data === parent(B).data)
-                @test BB[1:10, 1:10] == _B[1:10, 1:10]
+                @test BB.dv.data === BB.ev.data
+                @test parent(BB).dv.data.datasize == parent(_B).dv.data.datasize
+                @test !(BB === B) && !(parent(BB).dv.data === parent(B).dv.data)
+                @test BB[1:100, 1:100] == _B[1:100, 1:100]
+                @test BB[1:2:50, 1:3:40] == _B[1:2:50, 1:3:40]
+                @test view(BB, [1, 3, 7, 10], 1:10) == _B[[1, 3, 7, 10], 1:10]
             end
             @test LazyBandedMatrices.bidiagonaluplo(B) == uplo
             @test LazyBandedMatrices.Bidiagonal(B) === LazyBandedMatrices.Bidiagonal(B.dv, B.ev, Symbol(uplo))
@@ -42,19 +52,21 @@
             @test B[band(0)][1:100] == B.dv[1:100]
             if uplo == 'U'
                 @test B[band(1)][1:100] == B.ev[1:100]
-                @test B[band(-1)][1:100] == zeros(100)
+                # @test B[band(-1)][1:100] == zeros(100) # This test requires that we define a 
+                                                         # convert(::Type{BidiagonalConjugationBand{T}}, ::Zeros{V, 1, Tuple{OneToInf{Int}}}) where {T, V} method, 
+                                                         # which we probably don't need beyond this test
             else
                 @test B[band(-1)][1:100] == B.ev[1:100]
-                @test B[band(1)][1:100] == zeros(100)
+                # @test B[band(1)][1:100] == zeros(100)
             end
-            @test B.dv[500] == B.data.dv[500]
-            @test B.datasize == 1000 # make sure wrapper is working correctly 
-            @test B.ev[1005] == B.data.ev[1005]
-            @test B.datasize == 2010 + 2 * (uplo == 'U')
+            @test B.dv[500] == B.dv.data.dv[500]
+            @test B.dv.data.datasize == 1002 
+            @test B.ev[1005] == B.ev.data.ev[1005]
+            @test B.ev.data.datasize == 2012 # 2010 + 2 * (uplo == 'L')
             # @test_broken ApplyArray(inv, B)[1:100, 1:100] ≈ ApplyArray(inv, A)[1:100, 1:100] # need to somehow let inv (or even ApplyArray(inv, )) work
-            @test (B+B)[1:100, 1:100] ≈ 2A[1:100, 1:100] ≈ 2B[1:100, 1:100]
-            @test (B*I)[1:100, 1:100] ≈ B[1:100, 1:100]
-            @test (B*Diagonal(1:∞))[1:100, 1:100] ≈ B[1:100, 1:100] * Diagonal(1:100)
+            # @test (B+B)[1:100, 1:100] ≈ 2A[1:100, 1:100] ≈ 2B[1:100, 1:100]
+            # @test (B*I)[1:100, 1:100] ≈ B[1:100, 1:100]
+            # @test (B*Diagonal(1:∞))[1:100, 1:100] ≈ B[1:100, 1:100] * Diagonal(1:100)
         end
     end
 end
\ No newline at end of file

From 550a03a36ea66bc5f28c3f8188c9eb64622bbc22 Mon Sep 17 00:00:00 2001
From: DanielVandH <danj.vandenheuvel@gmail.com>
Date: Sat, 6 Jul 2024 13:27:14 +0100
Subject: [PATCH 07/13] Function barriers

---
 src/banded/bidiagonalconjugation.jl | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/src/banded/bidiagonalconjugation.jl b/src/banded/bidiagonalconjugation.jl
index d68b51f..f98936b 100644
--- a/src/banded/bidiagonalconjugation.jl
+++ b/src/banded/bidiagonalconjugation.jl
@@ -36,8 +36,7 @@ function copy(data::BidiagonalConjugationData)
     return BidiagonalConjugationData(copy(U), copy(C), copy(dv), copy(ev), uplo, datasize)
 end
 
-function _compute_column_up!(data::BidiagonalConjugationData, i)
-    U, C = data.U, data.C
+function _compute_column_up!(data::BidiagonalConjugationData, U, C, i)
     dv, ev = data.dv, data.ev
     if i == 1
         dv[i] = C[1, 1] / U[1, 1]
@@ -51,8 +50,7 @@ function _compute_column_up!(data::BidiagonalConjugationData, i)
     return data
 end
 
-function _compute_column_lo!(data::BidiagonalConjugationData, i)
-    U, C = data.U, data.C
+function _compute_column_lo!(data::BidiagonalConjugationData, U, C, i)
     dv, ev = data.dv, data.ev
     uᵢᵢ, uᵢ₊₁ᵢ, uᵢᵢ₊₁, uᵢ₊₁ᵢ₊₁ = U[i, i], U[i+1, i], U[i, i+1], U[i+1, i+1]
     cᵢᵢ, cᵢ₊₁ᵢ = C[i, i], C[i+1, i]
@@ -65,8 +63,14 @@ end
 function _compute_columns!(data::BidiagonalConjugationData, i)
     ds = data.datasize
     up = data.uplo == 'U'
+    U, C = data.U, data.C # Treat _compute_column_(up/lo) as function barriers and take these out early
+    return __compute_columns!(data, U, C, i)
+end
+function __compute_columns!(data::BidiagonalConjugationData, U, C, i)
+    ds = data.datasize 
+    up = data.uplo == 'U'
     for j in (ds+1):i
-        up ? _compute_column_up!(data, j) : _compute_column_lo!(data, j)
+        up ? _compute_column_up!(data, U, C, j) : _compute_column_lo!(data, U, C, j)
     end
     data.datasize = i
     return data

From c45472f82863e2c77de826271c69c4825001f145 Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Sat, 6 Jul 2024 13:39:07 +0100
Subject: [PATCH 08/13] BidiagonalConjugationBand should not be <:
 AbstractCachedVector

---
 src/banded/bidiagonalconjugation.jl | 2 +-
 test/test_bidiagonalconjugation.jl  | 5 +++++
 2 files changed, 6 insertions(+), 1 deletion(-)

diff --git a/src/banded/bidiagonalconjugation.jl b/src/banded/bidiagonalconjugation.jl
index f98936b..5d541e6 100644
--- a/src/banded/bidiagonalconjugation.jl
+++ b/src/banded/bidiagonalconjugation.jl
@@ -85,7 +85,7 @@ function resizedata!(data::BidiagonalConjugationData, n)
     return _compute_columns!(data, 2n)
 end
 
-struct BidiagonalConjugationBand{T} <: AbstractCachedVector{T}
+struct BidiagonalConjugationBand{T} <: LazyVector{T}
     data::BidiagonalConjugationData{T}
     diag::Bool # true => diagonal, false => offdiagonal 
 end
diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl
index 79db456..bad9396 100644
--- a/test/test_bidiagonalconjugation.jl
+++ b/test/test_bidiagonalconjugation.jl
@@ -1,3 +1,8 @@
+using InfiniteLinearAlgebra, InfiniteRandomArrays, BandedMatrices, LazyArrays, LazyBandedMatrices, InfiniteArrays, ArrayLayouts, Test
+using InfiniteLinearAlgebra: BidiagonalConjugation, OneToInf
+using ArrayLayouts: supdiagonaldata, subdiagonaldata, diagonaldata
+using LazyArrays: LazyLayout
+
 @testset "BidiagonalConjugationData" begin
     @test InfiniteLinearAlgebra._to_uplo('U') == 'U'
     @test InfiniteLinearAlgebra._to_uplo('L') == 'L'

From f9ded918c7c7a71af9e773dc608c900156393d70 Mon Sep 17 00:00:00 2001
From: DanielVandH <danj.vandenheuvel@gmail.com>
Date: Sat, 6 Jul 2024 15:01:32 +0100
Subject: [PATCH 09/13] oops

---
 src/banded/bidiagonalconjugation.jl | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/src/banded/bidiagonalconjugation.jl b/src/banded/bidiagonalconjugation.jl
index 5d541e6..e5c980f 100644
--- a/src/banded/bidiagonalconjugation.jl
+++ b/src/banded/bidiagonalconjugation.jl
@@ -109,11 +109,6 @@ end
 
 copy(band::BidiagonalConjugationBand) = band
 
-function convert(::Type{BidiagonalConjugationBand{T}}, ::Zeros{V, 1, Tuple{OneToInf{Int}}}) where {T, V}
-    # without this method, we can't do e.g. B[Band(-1)] for upper bidiagonal
-    
-end
-
 const BidiagonalConjugation{T} = Bidiagonal{T, BidiagonalConjugationBand{T}}
 
 """

From 833a6a9c6289b68c64413d53b83473a765cccd55 Mon Sep 17 00:00:00 2001
From: DanielVandH <danj.vandenheuvel@gmail.com>
Date: Sat, 6 Jul 2024 16:24:59 +0100
Subject: [PATCH 10/13] Changes

---
 src/banded/bidiagonalconjugation.jl |  3 +--
 test/test_bidiagonalconjugation.jl  | 31 ++++++++++++++++++++---------
 2 files changed, 23 insertions(+), 11 deletions(-)

diff --git a/src/banded/bidiagonalconjugation.jl b/src/banded/bidiagonalconjugation.jl
index e5c980f..165cdb3 100644
--- a/src/banded/bidiagonalconjugation.jl
+++ b/src/banded/bidiagonalconjugation.jl
@@ -101,8 +101,7 @@ function _bcb_getindex(band::BidiagonalConjugationBand, I)
     end
 end
 
-@inline cache_getindex(band::BidiagonalConjugationBand, I::Integer) = _bcb_getindex(band, I) # needed for show
-@inline getindex(band::BidiagonalConjugationBand, I::Integer) = cache_getindex(band, I) # also needed for show because of isassigned
+@inline getindex(band::BidiagonalConjugationBand, I::Integer) = _bcb_getindex(band, I)
 @inline getindex(band::BidiagonalConjugationBand, I::AbstractVector) = _bcb_getindex(band, I)
 #@inline getindex(band::BidiagonalConjugationBand, I::AbstractInfUnitRange{<:Integer}) = view(band, I)
 #@inline getindex(band::SubArray{<:Any, 1, <:BidiagonalConjugationBand}, I::AbstractInfUnitRange{<:Integer}) = view(band, I)
diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl
index bad9396..956b7b2 100644
--- a/test/test_bidiagonalconjugation.jl
+++ b/test/test_bidiagonalconjugation.jl
@@ -1,6 +1,7 @@
 using InfiniteLinearAlgebra, InfiniteRandomArrays, BandedMatrices, LazyArrays, LazyBandedMatrices, InfiniteArrays, ArrayLayouts, Test
 using InfiniteLinearAlgebra: BidiagonalConjugation, OneToInf
 using ArrayLayouts: supdiagonaldata, subdiagonaldata, diagonaldata
+using LinearAlgebra
 using LazyArrays: LazyLayout
 
 @testset "BidiagonalConjugationData" begin
@@ -11,7 +12,7 @@ using LazyArrays: LazyLayout
     @test InfiniteLinearAlgebra._to_uplo(:L) == 'L'
     @test_throws ArgumentError InfiniteLinearAlgebra._to_uplo(:a)
 
-    for _ in 1:10
+    for _ in 1:3
         V1 = InfRandTridiagonal()
         A1 = InfRandBidiagonal('U')
         X1 = brand(∞, 0, 2)
@@ -58,20 +59,32 @@ using LazyArrays: LazyLayout
             if uplo == 'U'
                 @test B[band(1)][1:100] == B.ev[1:100]
                 # @test B[band(-1)][1:100] == zeros(100) # This test requires that we define a 
-                                                         # convert(::Type{BidiagonalConjugationBand{T}}, ::Zeros{V, 1, Tuple{OneToInf{Int}}}) where {T, V} method, 
-                                                         # which we probably don't need beyond this test
+                # convert(::Type{BidiagonalConjugationBand{T}}, ::Zeros{V, 1, Tuple{OneToInf{Int}}}) where {T, V} method, 
+                # which we probably don't need beyond this test
             else
                 @test B[band(-1)][1:100] == B.ev[1:100]
                 # @test B[band(1)][1:100] == zeros(100)
             end
             @test B.dv[500] == B.dv.data.dv[500]
-            @test B.dv.data.datasize == 1002 
+            @test B.dv.data.datasize == 1002
             @test B.ev[1005] == B.ev.data.ev[1005]
-            @test B.ev.data.datasize == 2012 # 2010 + 2 * (uplo == 'L')
-            # @test_broken ApplyArray(inv, B)[1:100, 1:100] ≈ ApplyArray(inv, A)[1:100, 1:100] # need to somehow let inv (or even ApplyArray(inv, )) work
-            # @test (B+B)[1:100, 1:100] ≈ 2A[1:100, 1:100] ≈ 2B[1:100, 1:100]
-            # @test (B*I)[1:100, 1:100] ≈ B[1:100, 1:100]
-            # @test (B*Diagonal(1:∞))[1:100, 1:100] ≈ B[1:100, 1:100] * Diagonal(1:100)
+            @test B.ev.data.datasize == 2012
+            @test ApplyArray(inv, B)[1:100, 1:100] ≈ ApplyArray(inv, A)[1:100, 1:100] # need to somehow let inv (or even ApplyArray(inv, )) work
+            @test (B+B)[1:100, 1:100] ≈ 2A[1:100, 1:100] ≈ 2B[1:100, 1:100]
+            @test (B*I)[1:100, 1:100] ≈ B[1:100, 1:100]
+            # @test (B*Diagonal(1:∞))[1:100, 1:100] ≈ B[1:100, 1:100] * Diagonal(1:100) # Uncomment once https://github.com/JuliaLinearAlgebra/ArrayLayouts.jl/pull/241 is registered
+
+            # Pointwise tests 
+            for i in 1:10
+                for j in 1:10
+                    @test B[i, j] ≈ A[i, j]
+                    @test B'[i, j] ≈ A[j, i]
+                end
+            end
+            @inferred B[5, 5]
+
+            # Make sure that, when indexing the transpose, B expands correctly
+            @test B'[3000:3005, 2993:3006] ≈ A[2993:3006, 3000:3005]'
         end
     end
 end
\ No newline at end of file

From e9f33d96f2947ac7042151375bd5face25bc8178 Mon Sep 17 00:00:00 2001
From: DanielVandH <danj.vandenheuvel@gmail.com>
Date: Sun, 7 Jul 2024 14:17:47 +0100
Subject: [PATCH 11/13] Make copy a no-op

---
 src/banded/bidiagonalconjugation.jl | 5 +----
 test/test_bidiagonalconjugation.jl  | 2 +-
 2 files changed, 2 insertions(+), 5 deletions(-)

diff --git a/src/banded/bidiagonalconjugation.jl b/src/banded/bidiagonalconjugation.jl
index 165cdb3..b6eef5f 100644
--- a/src/banded/bidiagonalconjugation.jl
+++ b/src/banded/bidiagonalconjugation.jl
@@ -133,9 +133,6 @@ function _BidiagonalConjugation(data, uplo) # need uplo argument so that we can
     return Bidiagonal(dv, ev, uplo)
 end
 
-function copy(A::BidiagonalConjugation) # need a new method so that the data remains aliased between dv and ev 
-    data = copy(A.dv.data)
-    return _BidiagonalConjugation(data, A.uplo)
-end
+copy(A::BidiagonalConjugation)  = A # no-op
 
 LazyBandedMatrices.Bidiagonal(A::BidiagonalConjugation) = LazyBandedMatrices.Bidiagonal(A.dv, A.ev, A.uplo)
\ No newline at end of file
diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl
index 956b7b2..b8e902f 100644
--- a/test/test_bidiagonalconjugation.jl
+++ b/test/test_bidiagonalconjugation.jl
@@ -45,7 +45,7 @@ using LazyArrays: LazyLayout
                 BB = copy(_B)
                 @test BB.dv.data === BB.ev.data
                 @test parent(BB).dv.data.datasize == parent(_B).dv.data.datasize
-                @test !(BB === B) && !(parent(BB).dv.data === parent(B).dv.data)
+                # @test !(BB === B) && !(parent(BB).dv.data === parent(B).dv.data) # copy is a no-op
                 @test BB[1:100, 1:100] == _B[1:100, 1:100]
                 @test BB[1:2:50, 1:3:40] == _B[1:2:50, 1:3:40]
                 @test view(BB, [1, 3, 7, 10], 1:10) == _B[[1, 3, 7, 10], 1:10]

From a74a47ef29f074fa57d28ed214d8f98c8d730ae3 Mon Sep 17 00:00:00 2001
From: DanielVandH <danj.vandenheuvel@gmail.com>
Date: Mon, 8 Jul 2024 14:50:55 +0100
Subject: [PATCH 12/13] Cleanup

---
 src/banded/bidiagonalconjugation.jl | 16 ++++++----------
 1 file changed, 6 insertions(+), 10 deletions(-)

diff --git a/src/banded/bidiagonalconjugation.jl b/src/banded/bidiagonalconjugation.jl
index b6eef5f..7c5de69 100644
--- a/src/banded/bidiagonalconjugation.jl
+++ b/src/banded/bidiagonalconjugation.jl
@@ -61,13 +61,11 @@ function _compute_column_lo!(data::BidiagonalConjugationData, U, C, i)
 end
 
 function _compute_columns!(data::BidiagonalConjugationData, i)
-    ds = data.datasize
-    up = data.uplo == 'U'
     U, C = data.U, data.C # Treat _compute_column_(up/lo) as function barriers and take these out early
     return __compute_columns!(data, U, C, i)
 end
 function __compute_columns!(data::BidiagonalConjugationData, U, C, i)
-    ds = data.datasize 
+    ds = data.datasize
     up = data.uplo == 'U'
     for j in (ds+1):i
         up ? _compute_column_up!(data, U, C, j) : _compute_column_lo!(data, U, C, j)
@@ -89,26 +87,24 @@ struct BidiagonalConjugationBand{T} <: LazyVector{T}
     data::BidiagonalConjugationData{T}
     diag::Bool # true => diagonal, false => offdiagonal 
 end
-@inline size(A::BidiagonalConjugationBand) = (ℵ₀,)
+@inline size(::BidiagonalConjugationBand) = (ℵ₀,)
 @inline resizedata!(A::BidiagonalConjugationBand, n) = resizedata!(A.data, n)
 
 function _bcb_getindex(band::BidiagonalConjugationBand, I)
     resizedata!(band, maximum(I) + 1)
-    if band.diag 
+    if band.diag
         return band.data.dv[I]
-    else 
+    else
         return band.data.ev[I]
     end
 end
 
 @inline getindex(band::BidiagonalConjugationBand, I::Integer) = _bcb_getindex(band, I)
 @inline getindex(band::BidiagonalConjugationBand, I::AbstractVector) = _bcb_getindex(band, I)
-#@inline getindex(band::BidiagonalConjugationBand, I::AbstractInfUnitRange{<:Integer}) = view(band, I)
-#@inline getindex(band::SubArray{<:Any, 1, <:BidiagonalConjugationBand}, I::AbstractInfUnitRange{<:Integer}) = view(band, I)
 
 copy(band::BidiagonalConjugationBand) = band
 
-const BidiagonalConjugation{T} = Bidiagonal{T, BidiagonalConjugationBand{T}}
+const BidiagonalConjugation{T} = Bidiagonal{T,BidiagonalConjugationBand{T}}
 
 """
     BidiagonalConjugation(U, X, V, uplo)
@@ -133,6 +129,6 @@ function _BidiagonalConjugation(data, uplo) # need uplo argument so that we can
     return Bidiagonal(dv, ev, uplo)
 end
 
-copy(A::BidiagonalConjugation)  = A # no-op
+copy(A::BidiagonalConjugation) = A # no-op
 
 LazyBandedMatrices.Bidiagonal(A::BidiagonalConjugation) = LazyBandedMatrices.Bidiagonal(A.dv, A.ev, A.uplo)
\ No newline at end of file

From ebc1d2f512463f3c0251b2580c1756f7a1f7e35b Mon Sep 17 00:00:00 2001
From: DanielVandH <danj.vandenheuvel@gmail.com>
Date: Fri, 12 Jul 2024 12:24:01 +0100
Subject: [PATCH 13/13] Fix resize

---
 src/banded/bidiagonalconjugation.jl | 14 +++++++++-----
 test/test_bidiagonalconjugation.jl  | 12 ++++++------
 2 files changed, 15 insertions(+), 11 deletions(-)

diff --git a/src/banded/bidiagonalconjugation.jl b/src/banded/bidiagonalconjugation.jl
index 7c5de69..57eaa06 100644
--- a/src/banded/bidiagonalconjugation.jl
+++ b/src/banded/bidiagonalconjugation.jl
@@ -75,12 +75,16 @@ function __compute_columns!(data::BidiagonalConjugationData, U, C, i)
 end
 
 function resizedata!(data::BidiagonalConjugationData, n)
-    ds = data.datasize
-    n ≤ ds && return data
+    n ≤ 0 && return data
+    v = data.datasize
+    n = max(v, n)
     dv, ev = data.dv, data.ev
-    resize!(dv, 2n + 1)
-    resize!(ev, 2n)
-    return _compute_columns!(data, 2n)
+    if n > length(ev) # Avoid O(n²) growing. Note min(length(dv), length(ev)) == length(ev)
+        resize!(dv, 2n + 1)
+        resize!(ev, 2n)
+    end
+    n > v && _compute_columns!(data, n)
+    return data
 end
 
 struct BidiagonalConjugationBand{T} <: LazyVector{T}
diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl
index b8e902f..510d180 100644
--- a/test/test_bidiagonalconjugation.jl
+++ b/test/test_bidiagonalconjugation.jl
@@ -17,13 +17,13 @@ using LazyArrays: LazyLayout
         A1 = InfRandBidiagonal('U')
         X1 = brand(∞, 0, 2)
         U1 = X1 * V1 * ApplyArray(inv, A1)
-        B1 = BidiagonalConjugation(U1, X1, V1, 'U')
+        B1 = BidiagonalConjugation(U1, X1, V1, 'U');
 
         V2 = brand(∞, 0, 1)
         A2 = LazyBandedMatrices.Bidiagonal(Fill(0.2, ∞), 2.0 ./ (1.0 .+ (1:∞)), 'L') # LinearAlgebra.Bidiagonal not playing nice for this case
         X2 = InfRandBidiagonal('L')
-        U2 = X2 * V2 * ApplyArray(inv, A2) # U2 isn't upper Hessenberg (actually it's lower, oh well), doesn't seem to matter for computing A
-        B2 = BidiagonalConjugation(U2, X2, V2, :L)
+        U2 = X2 * V2 * ApplyArray(inv, A2) 
+        B2 = BidiagonalConjugation(U2, X2, V2, :L);
 
         for (A, B, uplo) in ((A1, B1, 'U'), (A2, B2, 'L'))
             @test B.dv.data === B.ev.data
@@ -66,9 +66,9 @@ using LazyArrays: LazyLayout
                 # @test B[band(1)][1:100] == zeros(100)
             end
             @test B.dv[500] == B.dv.data.dv[500]
-            @test B.dv.data.datasize == 1002
+            @test B.dv.data.datasize == 501
             @test B.ev[1005] == B.ev.data.ev[1005]
-            @test B.ev.data.datasize == 2012
+            @test B.ev.data.datasize == 1006
             @test ApplyArray(inv, B)[1:100, 1:100] ≈ ApplyArray(inv, A)[1:100, 1:100] # need to somehow let inv (or even ApplyArray(inv, )) work
             @test (B+B)[1:100, 1:100] ≈ 2A[1:100, 1:100] ≈ 2B[1:100, 1:100]
             @test (B*I)[1:100, 1:100] ≈ B[1:100, 1:100]
@@ -87,4 +87,4 @@ using LazyArrays: LazyLayout
             @test B'[3000:3005, 2993:3006] ≈ A[2993:3006, 3000:3005]'
         end
     end
-end
\ No newline at end of file
+end