work on banded adaptive qr

Author: dlfivefifty

src/InfiniteLinearAlgebra.jl

@@ -7,7 +7,7 @@ import Base: +, -, *, /, \, ^, OneTo, getindex, promote_op, _unsafe_getindex, pr
 import Base.Broadcast: BroadcastStyle
 
 import InfiniteArrays: OneToInf, InfUnitRange, Infinity, InfStepRange
-import FillArrays: AbstractFill
+import FillArrays: AbstractFill, getindex_value
 import BandedMatrices: BandedMatrix, _BandedMatrix, bandeddata, bandwidths
 import LinearAlgebra: lmul!,  ldiv!, matprod, qr, AbstractTriangular, AbstractQ, adjoint, transpose
 import LazyArrays: CachedArray, CachedMatrix, CachedVector, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes, ApplyStyle, LazyArrayApplyStyle, LazyArrayStyle,
@@ -17,7 +17,7 @@ import MatrixFactorizations: ql, ql!, QLPackedQ, getL, getR, reflector!, reflect
 
 import BlockArrays: BlockSizes, cumulsizes, _find_block, AbstractBlockVecOrMat, sizes_from_blocks
 
-import BandedMatrices: BandedMatrix, bandwidths, AbstractBandedLayout, _banded_qr!
+import BandedMatrices: BandedMatrix, bandwidths, AbstractBandedLayout, _banded_qr!, _banded_qr
 
 import BlockBandedMatrices: _BlockSkylineMatrix, _BandedMatrix, AbstractBlockSizes, cumulsizes, _BlockSkylineMatrix, BlockSizes, blockstart, blockstride,
         BlockSkylineSizes, BlockSkylineMatrix, BlockBandedMatrix, _BlockBandedMatrix, BlockTridiagonal

src/banded/infbanded.jl

@@ -8,12 +8,17 @@ const TriPertToeplitz{T} = Tridiagonal{T,Vcat{T,1,Tuple{Vector{T},Fill{T,1,Tuple
 const AdjTriPertToeplitz{T} = Adjoint{T,Tridiagonal{T,Vcat{T,1,Tuple{Vector{T},Fill{T,1,Tuple{OneToInf{Int}}}}}}}
 const InfBandedMatrix{T,V<:AbstractMatrix{T}} = BandedMatrix{T,V,OneToInf{Int}}
 
+_prepad(p, a) = Vcat(Zeros(max(p,0)), a)
+_prepad(p, a::Zeros{T,1}) where T = Zeros{T}(length(a)+p)
+_prepad(p, a::Ones{T,1}) where T = Ones{T}(length(a)+p)
+_prepad(p, a::AbstractFill{T,1}) where T = Fill{T}(getindex_value(a), length(a)+p)
+
 function BandedMatrix{T}(kv::Tuple{Vararg{Pair{<:Integer,<:AbstractVector}}},
                          ::NTuple{2,Infinity},
                          (l,u)::NTuple{2,Integer}) where T
     ks = getproperty.(kv, :first)
-    rws = Vcat(permutedims.(getproperty.(kv, :second))...)
     l,u = -minimum(ks),maximum(ks)
+    rws = Vcat(permutedims.(_prepad.(ks,getproperty.(kv, :second)))...)
     c = zeros(T, l+u+1, length(ks))
     for (k,j) in zip(u .- ks .+ 1,1:length(ks))
         c[k,j] = one(T)

src/infqr.jl

@@ -15,12 +15,12 @@ AdaptiveQRData(A::AbstractMatrix{T}) where T = AdaptiveQRData(MemoryLayout(typeo
 function partialqr!(F::AdaptiveQRData{<:Any,<:BandedMatrix}, n::Int)
     if n > F.ncols 
         l,u = bandwidths(F.data.array)
-        resizedata!(F.data,n+l,n+u+l)
-        resize!(F.τ,n)
+        resizedata!(F.data,n+l,n+u+l);
+        resize!(F.τ,n);
         ñ = F.ncols
-        factors = view(F.data.data,ñ+1:n+l,ñ+1:n)
-        τ = view(F.τ,ñ+1:n)
-        _banded_qr!(factors, τ)
+        factors = view(F.data.data,ñ+1:n+l,ñ+1:n);
+        τ = view(F.τ,ñ+1:n);
+        _banded_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))
@@ -78,7 +78,7 @@ end
 getR(Q::QR, ::Tuple{OneToInf{Int},OneToInf{Int}}) where T = UpperTriangular(Q.factors)
 
 
-function qr(A::BandedMatrix{<:Any,<:Any,<:OneToInf}) 
+function _banded_qr(::NTuple{2,OneToInf{Int}}, A)
     data = AdaptiveQRData(A)
     QR(AdaptiveQRFactors(data), AdaptiveQRTau(data))
 end
@@ -182,4 +182,5 @@ ldiv!(dest::AbstractVector, F::QR{<:Any,<:AdaptiveQRFactors}, b::AbstractVector)
 ldiv!(F::QR{<:Any,<:AdaptiveQRFactors}, b::AbstractVector) = ldiv!(F.R, lmul!(F.Q',b))
 \(F::QR{<:Any,<:AdaptiveQRFactors}, B::AbstractVector) = ldiv!(F.R, F.Q'B)
 
+
 factorize(A::BandedMatrix{<:Any,<:Any,<:OneToInf}) = qr(A)

test/runtests.jl

@@ -9,16 +9,15 @@ import MatrixFactorizations: QLPackedQ
 import BandedMatrices: bandeddata, _BandedMatrix
 import LazyArrays: colsupport, ApplyStyle, MemoryLayout
 
-
-
-
-
 @testset "Algebra" begin 
     @testset "BandedMatrix" begin
         A = BandedMatrix(-3 => Fill(7/10,∞), -2 => 1:∞, 1 => Fill(2im,∞))
         @test A isa BandedMatrix{ComplexF64}
         @test A[1:10,1:10] == diagm(-3 => Fill(7/10,7), -2 => 1:8, 1 => Fill(2im,9))
 
+        A = BandedMatrix(0 => Vcat([1,2,3],Zeros(∞)), 1 => Vcat(1, Zeros(∞)))
+        @test A[1,2] == 1
+
         A = BandedMatrix(-3 => Fill(7/10,∞), -2 => Fill(1,∞), 1 => Fill(2im,∞))
         Ac = BandedMatrix(A')
         At = BandedMatrix(transpose(A))
@@ -112,4 +111,4 @@ end
 
 include("test_hessenbergq.jl")
 include("test_infql.jl")
-
+include("test_infqr.jl")

test/test_infqr.jl

@@ -92,6 +92,28 @@ import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout
         J = A \ Vcat([besselj(1,z)], Zeros(∞))
         @test J[1:20_000] ≈ [besselj(k,z) for k=0:20_000-1]
     end
-end
 
-    
+    @testset "5-band" begin
+        A = BandedMatrix(-2 => Ones(∞), -1 => Vcat(1, Zeros(∞)), 0 => Vcat([1,2,3],Zeros(∞)).+3, 1 => Vcat(1, Zeros(∞)), 2 => Ones(∞))
+        b = Vcat([3,4,5],Zeros(∞))
+        x = qr(A) \ b
+        @test x[1:2000] ≈ (A[1:2000,1:2000]\b[1:2000])
+    end
+
+    @testset "broadcast" begin
+        A = BandedMatrix(0 => -2*(0:∞), 1 => Ones(∞), -1 => Ones(∞))
+        B = BandedMatrix(-2 => Ones(∞), -1 => Vcat(1, Zeros(∞)), 0 => Vcat([1,2,3],Zeros(∞)).+3, 1 => Vcat(1, Zeros(∞)), 2 => Ones(∞))
+
+        AB = BroadcastArray(+,A,B)
+        C = cache(AB);
+        resizedata!(C,103,100); resizedata!(C,203,200);
+        @test C[103,104] ≈ 1.0
+        F = qr(AB);
+        partialqr!(F.factors.data, 100);
+        partialqr!(F.factors.data, 200);
+        @test norm(F.factors.data.data.data) ≤ 4000
+        b = Vcat([3,4,5],Zeros(∞))
+        @time x = qr(AB) \ b;
+        @test x[1:300] ≈ AB[1:300,1:300] \ b[1:300]
+    end
+end