Add AdaptiveBlock/BandedLayout (#195)

* Add AdaptiveBlock/BandedLayout * Fix tests * Remove over-templating * start Y + I tests * Y+I fixes * avoid ambiguities * LazyArrays v2.6 * increase coverage
JuliaLinearAlgebra · Feb 12, 2025 · 768594e · 768594e · dlfivefifty · Feb 12, 2025
1 parent 84517d2
commit 768594e
Show file tree

Hide file tree

Showing 7 changed files with 114 additions and 53 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,6 +1,6 @@
 name = "InfiniteLinearAlgebra"
 uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
-version = "0.9.1"
+version = "0.10"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -26,7 +26,7 @@ FillArrays = "1.0"
 InfiniteArrays = "0.15"
 InfiniteRandomArrays = "0.2"
 Infinities = "0.1"
-LazyArrays = "2.5"
+LazyArrays = "2.6"
 LazyBandedMatrices = "0.11"
 LinearAlgebra = "1"
 MatrixFactorizations = "3.0"

diff --git a/src/InfiniteLinearAlgebra.jl b/src/InfiniteLinearAlgebra.jl
@@ -40,16 +40,19 @@ 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, cache_getindex, cache_layout
+                   simplify, islazy, islazy_layout, cache_getindex, cache_layout, AbstractInvLayout
 
-import LazyBandedMatrices: AbstractLazyBandedBlockBandedLayout, AbstractLazyBandedLayout, ApplyBandedLayout, BlockVec,
+import LazyBandedMatrices: AbstractLazyBandedBlockBandedLayout, AbstractLazyBandedLayout, AbstractLazyBlockBandedLayout, ApplyBandedLayout, BlockVec,
                            BroadcastBandedLayout, KronTravBandedBlockBandedLayout, LazyBandedLayout,
                            _block_interlace_axes, _krontrav_axes, krontravargs
 
 const StructuredLayoutTypes{Lay} = Union{SymmetricLayout{Lay}, HermitianLayout{Lay}, TriangularLayout{'L','N',Lay}, TriangularLayout{'U','N',Lay}, TriangularLayout{'L','U',Lay}, TriangularLayout{'U','U',Lay}}
 
-const BandedLayouts = Union{AbstractBandedLayout, StructuredLayoutTypes{<:AbstractBandedLayout}}
-
+LazyArraysBandedMatricesExt = Base.get_extension(LazyArrays, :LazyArraysBandedMatricesExt)
+
+const BandedLazyLayouts = LazyArraysBandedMatricesExt.BandedLazyLayouts
+const BandedLayouts = LazyArraysBandedMatricesExt.BandedLayouts
+const BlockBandedLayouts = Union{AbstractBlockBandedLayout, BlockLayout{<:AbstractBandedLayout}, StructuredLayoutTypes{<:AbstractBlockBandedLayout}} 
 
 import LinearAlgebra: AbstractQ, AdjointQ, AdjOrTrans, factorize, matprod, qr
 
@@ -118,7 +121,7 @@ pad(c::Transpose, ax, bx) = transpose(pad(parent(c), bx, ax))
 pad(c::Adjoint, ax, bx) = adjoint(pad(parent(c), bx, ax))
 pad(c::BlockVec, ax::BlockedOneTo{Int,<:InfStepRange}) = BlockVec(pad(c.args[1], size(c.args[1],1), ∞))
 
-export ∞, ContinuousSpectrumError, BlockTridiagonal
+export ∞, ContinuousSpectrumError, BlockTridiagonal, TridiagonalConjugation, BidiagonalConjugation
 
 include("banded/hessenbergq.jl")
 

diff --git a/src/banded/tridiagonalconjugation.jl b/src/banded/tridiagonalconjugation.jl
@@ -290,14 +290,15 @@ getindex(K::SubArray{<:Any,1,<:TridiagonalConjugationBand}, k::AbstractInfUnitRa
 copy(A::TridiagonalConjugationBand) = A # immutable
 
 
-const TridiagonalConjugation{T} = Tridiagonal{T, TridiagonalConjugationBand{T}}
-const SymTridiagonalConjugation{T} = SymTridiagonal{T, TridiagonalConjugationBand{T}}
+# Use LazyBandedMatrices.Tridiagonal to support algebra
+const TridiagonalConjugation{T} = LazyBandedMatrices.Tridiagonal{T, TridiagonalConjugationBand{T}, TridiagonalConjugationBand{T}, TridiagonalConjugationBand{T}}
+const SymTridiagonalConjugation{T} = LazyBandedMatrices.SymTridiagonal{T, TridiagonalConjugationBand{T}, TridiagonalConjugationBand{T}}
 function TridiagonalConjugation(R, X, Y...)
     data = TridiagonalConjugationData(R, X, Y...)
-    Tridiagonal(TridiagonalConjugationBand(data, :dl), TridiagonalConjugationBand(data, :d), TridiagonalConjugationBand(data, :du))
+    LazyBandedMatrices.Tridiagonal(TridiagonalConjugationBand(data, :dl), TridiagonalConjugationBand(data, :d), TridiagonalConjugationBand(data, :du))
 end
 
 function SymTridiagonalConjugation(R, X, Y...)
     data = TridiagonalConjugationData(R, X, Y...)
-    SymTridiagonal(TridiagonalConjugationBand(data, :d), TridiagonalConjugationBand(data, :du))
+    LazyBandedMatrices.SymTridiagonal(TridiagonalConjugationBand(data, :d), TridiagonalConjugationBand(data, :du))
 end
diff --git a/src/infcholesky.jl b/src/infcholesky.jl
@@ -1,6 +1,6 @@
 
-mutable struct AdaptiveCholeskyFactors{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} <: LazyMatrix{T}
-    data::CachedMatrix{T,DM,M}
+mutable struct AdaptiveCholeskyFactors{T,DM<:AbstractMatrix{T}} <: LazyMatrix{T}
+    data::CachedMatrix{T,DM}
     ncols::Int
 end
 
@@ -16,12 +16,29 @@ function AdaptiveCholeskyFactors(::SymmetricBandedLayouts, S::AbstractMatrix{T})
     AdaptiveCholeskyFactors(CachedArray(data,A), 0)
 end
 AdaptiveCholeskyFactors(A::AbstractMatrix{T}) where T = AdaptiveCholeskyFactors(MemoryLayout(A), A)
-MemoryLayout(::Type{AdaptiveCholeskyFactors{T,DM,M}}) where {T,DM,M} = AdaptiveLayout{typeof(MemoryLayout(DM))}()
+
+struct AdaptiveCholeskyFactorsLayout <: AbstractLazyLayout end
+struct AdaptiveCholeskyFactorsBandedLayout <: AbstractLazyBandedLayout end
+struct AdaptiveCholeskyFactorsBlockBandedLayout <: AbstractLazyBlockBandedLayout end
+
+const AdaptiveCholeskyFactorsLayouts = Union{AdaptiveCholeskyFactorsLayout,AdaptiveCholeskyFactorsBandedLayout,AdaptiveCholeskyFactorsBlockBandedLayout}
+
+# TODO: support other than Banded
+# adaptivecholeskyfactorslayout(_) = AdaptiveCholeskyFactorsLayout()
+adaptivecholeskyfactorslayout(::BandedLayouts) = AdaptiveCholeskyFactorsBandedLayout()
+# adaptivecholeskyfactorslayout(::BlockBandedLayouts) = AdaptiveCholeskyFactorsBlockBandedLayout()
+
+
+MemoryLayout(::Type{AdaptiveCholeskyFactors{T,DM}}) where {T,DM} = adaptivecholeskyfactorslayout(MemoryLayout(DM))
+triangularlayout(::Type{Tri}, ::ML) where {Tri, ML<:AdaptiveCholeskyFactorsLayouts} = Tri{ML}()
+transposelayout(A::AdaptiveCholeskyFactorsLayouts) = A
 
 copy(A::AdaptiveCholeskyFactors) = AdaptiveCholeskyFactors(copy(A.data), copy(A.ncols))
 copy(A::Adjoint{T,<:AdaptiveCholeskyFactors}) where T = copy(parent(A))'
 copy(A::Transpose{T,<:AdaptiveCholeskyFactors}) where T = transpose(copy(parent(A)))
 
+cache_layout(::TriangularLayout{'U', 'N', AdaptiveCholeskyFactorsBandedLayout}, A::AbstractMatrix) = A # already cached
+
 function partialcholesky!(F::AdaptiveCholeskyFactors{T,<:BandedMatrix}, n::Int) where T
     if n > F.ncols
         _,u = bandwidths(F.data.array)
@@ -74,7 +91,7 @@ end
 colsupport(F::AdjOrTrans{<:Any,<:AdaptiveCholeskyFactors}, j) = rowsupport(parent(F), j)
 rowsupport(F::AdjOrTrans{<:Any,<:AdaptiveCholeskyFactors}, j) = colsupport(parent(F), j)
 
-function materialize!(M::MatLdivVec{<:TriangularLayout{'L','N',<:AdaptiveLayout},<:AbstractPaddedLayout})
+function materialize!(M::MatLdivVec{<:TriangularLayout{'L','N',<:AdaptiveCholeskyFactorsLayouts},<:AbstractPaddedLayout})
     A,B = M.A,M.B
     T = eltype(M)
     COLGROWTH = 1000 # rate to grow columns

diff --git a/src/infqr.jl b/src/infqr.jl
@@ -1,11 +1,11 @@
 
-mutable struct AdaptiveQRData{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}}
-    data::CachedMatrix{T,DM,M}
+mutable struct AdaptiveQRData{T,DM<:AbstractMatrix{T}}
+    data::CachedMatrix{T,DM}
     τ::Vector{T}
     ncols::Int
 end
 
-function AdaptiveQRData(::Union{SymmetricLayout{<:AbstractBandedLayout},AbstractBandedLayout}, A::AbstractMatrix{T}) where T
+function AdaptiveQRData(::BandedLayouts, A::AbstractMatrix{T}) where T
     l,u = bandwidths(A)
     FT = float(T)
     data = BandedMatrix{FT}(undef,(2l+u+1,0),(l,l+u)) # pad super
@@ -89,14 +89,24 @@ partialqr!(F::AdaptiveQRData{<:Any,<:BlockSkylineMatrix}, n::Int) =
     partialqr!(F, findblock(axes(F.data,2), n))
 
 
-struct AdaptiveQRFactors{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} <: LayoutMatrix{T}
-    data::AdaptiveQRData{T,DM,M}
+struct AdaptiveQRFactors{T,DM<:AbstractMatrix{T}} <: LayoutMatrix{T}
+    data::AdaptiveQRData{T,DM}
 end
 
-struct AdaptiveLayout{M} <: AbstractLazyLayout end
-MemoryLayout(::Type{AdaptiveQRFactors{T,DM,M}}) where {T,DM,M} = AdaptiveLayout{typeof(MemoryLayout(DM))}()
-triangularlayout(::Type{Tri}, ::ML) where {Tri, ML<:AdaptiveLayout} = Tri{ML}()
-transposelayout(A::AdaptiveLayout{ML}) where ML = AdaptiveLayout{typeof(transposelayout(ML()))}()
+struct AdaptiveQRFactorsLayout <: AbstractLazyLayout end
+struct AdaptiveQRFactorsBandedLayout <: AbstractLazyBandedLayout end
+struct AdaptiveQRFactorsBlockBandedLayout <: AbstractLazyBlockBandedLayout end
+
+const AdaptiveQRFactorsLayouts = Union{AdaptiveQRFactorsLayout,AdaptiveQRFactorsBandedLayout,AdaptiveQRFactorsBlockBandedLayout}
+
+adaptiveqrfactorslayout(_) = AdaptiveQRFactorsLayout()
+adaptiveqrfactorslayout(::BandedLayouts) = AdaptiveQRFactorsBandedLayout()
+adaptiveqrfactorslayout(::BlockBandedLayouts) = AdaptiveQRFactorsBlockBandedLayout()
+
+
+MemoryLayout(::Type{AdaptiveQRFactors{T,DM}}) where {T,DM} =adaptiveqrfactorslayout(MemoryLayout(DM))
+triangularlayout(::Type{Tri}, ::ML) where {Tri, ML<:AdaptiveQRFactorsLayouts} = Tri{ML}()
+# transposelayout(A::AdaptiveQRFactorsLayouts) = A
 
 size(F::AdaptiveQRFactors) = size(F.data.data)
 axes(F::AdaptiveQRFactors) = axes(F.data.data)
@@ -136,8 +146,8 @@ rowsupport(F::QRPackedQ{<:Any,<:AdaptiveQRFactors}, j) = first(rowsupport(F.fact
 blockcolsupport(F::QRPackedQ{<:Any,<:AdaptiveQRFactors}, j) = blockcolsupport(F.factors, j)
 
 
-struct AdaptiveQRTau{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} <: LayoutVector{T}
-    data::AdaptiveQRData{T,DM,M}
+struct AdaptiveQRTau{T,DM<:AbstractMatrix{T}} <: LayoutVector{T}
+    data::AdaptiveQRData{T,DM}
 end
 
 size(F::AdaptiveQRTau) = (size(F.data.data,1),)
@@ -164,7 +174,7 @@ factorize_layout(::BandedLayouts, ::NTuple{2,OneToInf{Int}}, A) = qr(A)
 factorize_layout(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = qr(A)
 
 
-cache_layout(::TriangularLayout{UPLO, UNIT, <:AdaptiveLayout}, A::AbstractMatrix) where {UPLO, UNIT} = A # already cached
+cache_layout(::TriangularLayout{'U', 'N', AdaptiveQRFactorsBandedLayout}, A::AbstractMatrix) = A # already cached
 
 partialqr!(F::QR, n) = partialqr!(F.factors, n)
 partialqr!(F::AdaptiveQRFactors, n) = partialqr!(F.data, n)
@@ -181,7 +191,7 @@ getindex(Q::QRPackedQ{<:Any,<:AdaptiveQRFactors,<:AdaptiveQRTau}, I::AbstractVec
 #########
 
 _view_QRPackedQ(A, kr, jr) = QRPackedQ(view(A.factors.data.data.data,kr,jr), view(A.τ.data.τ,jr))
-function materialize!(M::MatLmulVec{<:QRPackedQLayout{<:AdaptiveLayout},<:AbstractPaddedLayout})
+function materialize!(M::MatLmulVec{<:QRPackedQLayout{<:AdaptiveQRFactorsLayouts},<:AbstractPaddedLayout})
     A,B = M.A,M.B
     sB = size(paddeddata(B),1)
     partialqr!(A.factors.data,sB)
@@ -223,7 +233,7 @@ end
 
 _norm(x::Number) = abs(x)
 
-function materialize!(M::MatLmulVec{<:AdjQRPackedQLayout{<:AdaptiveLayout},<:AbstractPaddedLayout}; tolerance=floatmin(real(eltype(M))))
+function materialize!(M::MatLmulVec{<:AdjQRPackedQLayout{<:AdaptiveQRFactorsLayouts},<:AbstractPaddedLayout}; tolerance=floatmin(real(eltype(M))))
     adjA,B = M.A,M.B
     COLGROWTH = 1000 # rate to grow columns
 
@@ -274,7 +284,7 @@ function _view_QRPackedQ(A, KR::BlockRange, JR::BlockRange)
     QRPackedQ(view(A.factors.data.data.data,KR,JR), view(A.τ.data.τ,jr))
 end
 
-function materialize!(M::MatLmulVec{<:QRPackedQLayout{<:AdaptiveLayout{<:AbstractBlockBandedLayout}},<:AbstractPaddedLayout})
+function materialize!(M::MatLmulVec{QRPackedQLayout{AdaptiveQRFactorsBlockBandedLayout},<:AbstractPaddedLayout})
     A,B_in = M.A,M.B
     sB = length(paddeddata(B_in))
     ax1,ax2 = axes(A.factors.data.data)
@@ -290,7 +300,7 @@ function materialize!(M::MatLmulVec{<:QRPackedQLayout{<:AdaptiveLayout{<:Abstrac
     B
 end
 
-function materialize!(M::MatLmulVec{<:AdjQRPackedQLayout{<:AdaptiveLayout{<:AbstractBlockBandedLayout}},<:AbstractPaddedLayout}; tolerance=1E-30)
+function materialize!(M::MatLmulVec{<:AdjQRPackedQLayout{AdaptiveQRFactorsBlockBandedLayout},<:AbstractPaddedLayout}; tolerance=1E-30)
     adjA,B_in = M.A,M.B
     A = parent(adjA)
     T = eltype(M)
@@ -367,15 +377,19 @@ ldiv!(F::QR{<:Any,<:AdaptiveQRFactors}, b::LayoutVector; kwds...) = ldiv!(F.R, l
 factorize(A::BandedMatrix{<:Any,<:Any,<:OneToInf}) = qr(A)
 qr(A::SymTridiagonal{T,<:AbstractFill{T,1,Tuple{OneToInf{Int}}}}) where T = adaptiveqr(A)
 
-simplifiable(M::Mul{<:QRPackedQLayout{<:AdaptiveLayout}}) = Val(false)
-simplifiable(M::Mul{<:QRPackedQLayout{<:AdaptiveLayout},<:QRPackedQLayout{<:AdaptiveLayout}}) = Val(false)
-simplifiable(M::Mul{<:QRPackedQLayout{<:AdaptiveLayout},<:LazyLayouts}) = Val(false)
-simplifiable(M::Mul{<:Any,<:QRPackedQLayout{<:AdaptiveLayout}}) = Val(false)
-simplifiable(M::Mul{<:LazyLayouts,<:QRPackedQLayout{<:AdaptiveLayout}}) = Val(false)
-
-
-copy(M::Mul{<:QRPackedQLayout{<:AdaptiveLayout},<:QRPackedQLayout{<:AdaptiveLayout}}) = simplify(M)
-copy(M::Mul{<:QRPackedQLayout{<:AdaptiveLayout}}) = simplify(M)
-copy(M::Mul{<:QRPackedQLayout{<:AdaptiveLayout},<:LazyLayouts}) = simplify(M)
-copy(M::Mul{<:Any,<:QRPackedQLayout{<:AdaptiveLayout}}) = simplify(M)
-copy(M::Mul{<:LazyLayouts,<:QRPackedQLayout{<:AdaptiveLayout}}) = simplify(M)
+simplifiable(M::Mul{<:QRPackedQLayout{<:AdaptiveQRFactorsLayouts}}) = Val(false)
+simplifiable(M::Mul{<:QRPackedQLayout{<:AdaptiveQRFactorsLayouts},<:QRPackedQLayout{<:AdaptiveQRFactorsLayouts}}) = Val(false)
+simplifiable(M::Mul{<:QRPackedQLayout{<:AdaptiveQRFactorsLayouts},<:LazyLayouts}) = Val(false)
+simplifiable(M::Mul{<:Any,<:QRPackedQLayout{<:AdaptiveQRFactorsLayouts}}) = Val(false)
+simplifiable(M::Mul{<:LazyLayouts,<:QRPackedQLayout{<:AdaptiveQRFactorsLayouts}}) = Val(false)
+simplifiable(M::Mul{<:AbstractInvLayout,<:QRPackedQLayout{<:AdaptiveQRFactorsLayouts}}) = Val(false)
+
+
+copy(M::Mul{<:QRPackedQLayout{<:AdaptiveQRFactorsLayouts}, <:QRPackedQLayout{<:AdaptiveQRFactorsLayouts}}) = simplify(M)
+copy(M::Mul{<:QRPackedQLayout{<:AdaptiveQRFactorsLayouts}}) = simplify(M)
+copy(M::Mul{<:AdaptiveQRFactorsLayouts, <:BandedLazyLayouts}) = simplify(M)
+copy(M::Mul{<:QRPackedQLayout{<:AdaptiveQRFactorsLayouts}, <:BandedLazyLayouts}) = simplify(M)
+copy(M::Mul{<:QRPackedQLayout{<:AdaptiveQRFactorsLayouts}, <:LazyLayouts}) = simplify(M)
+copy(M::Mul{<:Any, <:QRPackedQLayout{<:AdaptiveQRFactorsLayouts}}) = simplify(M)
+copy(M::Mul{<:LazyLayouts, <:QRPackedQLayout{<:AdaptiveQRFactorsLayouts}}) = simplify(M)
+copy(M::Mul{<:BandedLazyLayouts, <:QRPackedQLayout{<:AdaptiveQRFactorsLayouts}}) = simplify(M)
diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl
@@ -185,12 +185,15 @@ end
         R,X = (_BandedMatrix(Vcat(-Ones(1,∞)/2, Zeros(1,∞), Hcat(Ones(1,1),Ones(1,∞)/2)), ℵ₀, 0,2),
                 LazyBandedMatrices.Tridiagonal(Vcat(1.0, Fill(1/2,∞)), Zeros(∞), Fill(1/2,∞)))
 
-        Y = TridiagonalConjugation(R, X)
-        n = 100_000
-        @test Y[n,n+1] ≈ 1/2
-
-        Y = SymTridiagonalConjugation(R, X)
-        n = 100_000
-        @test Y[n,n+1] ≈ 1/2
+        for Y in (TridiagonalConjugation(R, X), SymTridiagonalConjugation(R, X))
+            n = 100_000
+            @test Y[n,n+1] ≈ 1/2
+
+            @testset "Y+I" begin
+                @test (Y+I)[1:10,1:10] == (I+Y)[1:10,1:10] == I + Y[1:10,1:10]
+                @test (Y-I)[1:10,1:10] == Y[1:10,1:10] - I
+                @test (I-Y)[1:10,1:10] == I - Y[1:10,1:10]
+            end
+        end
     end
 end
diff --git a/test/test_infqr.jl b/test/test_infqr.jl
@@ -1,9 +1,9 @@
 using InfiniteLinearAlgebra, LinearAlgebra, BandedMatrices, InfiniteArrays, MatrixFactorizations, LazyArrays,
         FillArrays, SpecialFunctions, Test, SemiseparableMatrices, LazyBandedMatrices, BlockArrays, BlockBandedMatrices, ArrayLayouts
-import LazyArrays: colsupport, rowsupport, MemoryLayout, DenseColumnMajor, TriangularLayout, resizedata!, arguments
+import LazyArrays: colsupport, rowsupport, MemoryLayout, DenseColumnMajor, TriangularLayout, resizedata!, arguments, simplifiable
 import LazyBandedMatrices: BroadcastBandedLayout, InvDiagTrav, BroadcastBandedBlockBandedLayout
 import BandedMatrices: _BandedMatrix, _banded_qr!, BandedColumns
-import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout, adaptiveqr
+import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveQRFactorsBandedLayout, adaptiveqr
 import SemiseparableMatrices: AlmostBandedLayout, VcatAlmostBandedLayout
 
 
@@ -44,7 +44,7 @@ import SemiseparableMatrices: AlmostBandedLayout, VcatAlmostBandedLayout
         @testset "col/rowsupport" begin
             A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ℵ₀, 1, 1)
             F = qr(A)
-            @test MemoryLayout(typeof(F.factors)) isa AdaptiveLayout{BandedColumns{DenseColumnMajor}}
+            @test MemoryLayout(typeof(F.factors)) isa AdaptiveQRFactorsBandedLayout
             @test bandwidths(F.factors) == (1,2)
             @test colsupport(F.factors,1) ==  1:2
             @test colsupport(F.factors,5) ==  3:6
@@ -165,6 +165,29 @@ import SemiseparableMatrices: AlmostBandedLayout, VcatAlmostBandedLayout
             @test F\b ≈ F\view(b,:)
             @test_broken F\b ≈ ldiv!(F, view(copy(b),:))
         end
+
+        @testset "cache/mul" begin
+            A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ℵ₀, 1, 1)
+            Q,R = qr(A)
+            @test cache(R) ≡ R
+
+            @test simplifiable(*, R, A) == Val(false)
+            @test simplifiable(*, A, R) == Val(false)
+            @test simplifiable(*, parent(R), A) == Val(false)
+            @test simplifiable(*, A, parent(R)) == Val(false)
+            @test simplifiable(*, inv(R), A) == Val(false)
+            @test simplifiable(*, A, inv(R)) == Val(false)
+            @test simplifiable(*, inv(R), Q) == Val(false)
+            @test simplifiable(*, Q, inv(R)) == Val(false)
+            @test simplifiable(*, Symmetric(inv(R)), Q) == Val(false)
+            @test simplifiable(*, Q, Symmetric(inv(R))) == Val(false)
+            @test (R * A)[1:10,1:10] ≈ R[1:10,1:11] * A[1:11,1:10]
+            @test (Q * A)[1:10,1:10] ≈ Q[1:10,1:11] * A[1:11,1:10]
+            @test mul(Q, Symmetric(inv(R))) isa ApplyArray # several bugs avoided heree
+            @test mul(Symmetric(inv(R)), Q) isa ApplyArray # several bugs avoided heree
+            @test (parent(R) * A)[1:10,1:10] ≈ parent(R)[1:10,1:11] * A[1:11,1:10]
+
+        end
     end
 
     @testset "almost-banded" begin