Start almost banded adaptive QR

Author: dlfivefifty

Project.toml

@@ -16,13 +16,13 @@ MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 
 [compat]
-ArrayLayouts = "0.2"
-BandedMatrices = "0.15"
+ArrayLayouts = "0.2.1"
+BandedMatrices = "0.15.1"
 BlockArrays = "0.11"
 BlockBandedMatrices = "0.7.1"
-FillArrays = "0.8.4"
+FillArrays = "0.8.6"
 InfiniteArrays = "0.6.1, 0.7"
-LazyArrays = "0.15"
+LazyArrays = "0.16.1"
 LazyBandedMatrices = "0.2"
 julia = "1.3"
 

src/InfiniteLinearAlgebra.jl

@@ -29,6 +29,7 @@ import LazyBandedMatrices: MulBandedLayout, BroadcastBandedLayout
 import BlockBandedMatrices: _BlockSkylineMatrix, _BandedMatrix, _BlockSkylineMatrix, blockstart, blockstride,
         BlockSkylineSizes, BlockSkylineMatrix, BlockBandedMatrix, _BlockBandedMatrix, BlockTridiagonal
 
+import SemiseparableMatrices: AbstractAlmostBandedLayout, _almostbanded_qr!
 
 LazyArrays.@lazymul BandedMatrix{<:Any,<:Any,<:OneToInf}
 *(A::BandedMatrix{<:Any,<:Any,<:OneToInf}, b::CachedVector) = apply(*,A,b)

src/infqr.jl

@@ -10,6 +10,15 @@ function AdaptiveQRData(::AbstractBandedLayout, A::AbstractMatrix{T}) where T
     data = BandedMatrix{T}(undef,(2l+u+1,0),(l,l+u)) # pad super
     AdaptiveQRData(CachedArray(data,A), Vector{T}(), 0)
 end
+
+function AdaptiveQRData(::AbstractAlmostBandedLayout, A::AbstractMatrix{T}) where T 
+    l,u = almostbandwidths(A)
+    r = almostbandedrank(A)
+    data = AlmostBandedMatrix{T}(undef,(2l+u+1,0),(l,l+u),r) # pad super
+    
+    AdaptiveQRData(CachedArray(data,A), Vector{T}(), 0)
+end
+
 AdaptiveQRData(A::AbstractMatrix{T}) where T = AdaptiveQRData(MemoryLayout(typeof(A)), A)
 
 function partialqr!(F::AdaptiveQRData{<:Any,<:BandedMatrix}, n::Int)
@@ -34,6 +43,28 @@ function partialqr!(F::AdaptiveQRData{<:Any,<:BandedMatrix}, n::Int)
     F
 end
 
+function partialqr!(F::AdaptiveQRData{<:Any,<:AlmostBandedMatrix}, n::Int)
+    if n > F.ncols 
+        l,u = almostbandwidths(F.data.array)
+        resizedata!(F.data,n+l,n+u+l)
+        resize!(F.τ,n)
+        ñ = F.ncols
+        τ = view(F.τ,ñ+1:n)
+        if l ≤ 0 
+            zero!(τ)
+        else
+            factors = view(F.data.data,ñ+1:n+l,ñ+1:n)
+            _almostbanded_qr!(factors, τ)
+            # multiply remaining columns
+            n̄ = max(ñ+1,n-l-u+1) # first column interacting with extra data
+            Q = QRPackedQ(view(F.data.data,n̄:n+l,n̄:n), view(F.τ,n̄:n))
+            lmul!(Q',view(F.data.data,n̄:n+l,n+1:n+u+l))
+        end
+        F.ncols = n
+    end
+    F
+end
+
 struct AdaptiveQRFactors{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} <: AbstractMatrix{T}
     data::AdaptiveQRData{T,DM,M}
 end
@@ -82,11 +113,13 @@ end
 getR(Q::QR, ::Tuple{OneToInf{Int},OneToInf{Int}}) where T = UpperTriangular(Q.factors)
 
 
-function _banded_qr(::NTuple{2,OneToInf{Int}}, A)
+function adaptiveqr(A)
     data = AdaptiveQRData(A)
     QR(AdaptiveQRFactors(data), AdaptiveQRTau(data))
 end
 
+_banded_qr(::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A)
+
 #########
 # lmul!
 #########

test/test_infqr.jl

@@ -1,9 +1,10 @@
-using InfiniteLinearAlgebra, LinearAlgebra, BandedMatrices, InfiniteArrays, MatrixFactorizations, LazyArrays, FillArrays, SpecialFunctions, Test
+using InfiniteLinearAlgebra, LinearAlgebra, BandedMatrices, InfiniteArrays, MatrixFactorizations, LazyArrays,
+        FillArrays, SpecialFunctions, Test, SemiseparableMatrices
 import LazyArrays: colsupport, rowsupport, MemoryLayout, DenseColumnMajor, TriangularLayout, resizedata!, arguments
 import LazyBandedMatrices: BroadcastBandedLayout
 import BandedMatrices: _BandedMatrix, _banded_qr!, BandedColumns
-import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout
-import SemiseparableMatrices: AlmostBandedLayout
+import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout, adaptiveqr
+import SemiseparableMatrices: AlmostBandedLayout, VcatAlmostBandedLayout
 
 
 @testset "Adaptive QR" begin
@@ -136,8 +137,15 @@ import SemiseparableMatrices: AlmostBandedLayout
 
     @testset "almost-banded" begin
         A = Vcat(Ones(1,∞), BandedMatrix(0 => -Ones(∞), 1 => 1:∞))
-        @test MemoryLayout(typeof(A)) == AlmostBandedLayout()
-        B,C = arguments(A)
-        
-    end
+        @test MemoryLayout(typeof(A)) isa VcatAlmostBandedLayout
+        V = view(A,1:10,1:10)
+        @test MemoryLayout(typeof(V)) isa VcatAlmostBandedLayout
+        @test A[1:10,1:10] isa AlmostBandedMatrix
+        @test AlmostBandedMatrix(V) == Matrix(V) == A[1:10,1:10]
+
+        C = cache(A)
+        @test C[1000,1000] ≡ 999.0
+        F = adaptiveqr(A);
+        partialqr!(F.factors.data,10)
+    end 
 end