_qr

Author: dlfivefifty

Project.toml

@@ -24,6 +24,7 @@ FillArrays = "0.8.6"
 InfiniteArrays = "0.6.1, 0.7"
 LazyArrays = "0.16.1"
 LazyBandedMatrices = "0.2"
+MatrixFactorizations = "0.3.1"
 julia = "1.3"
 
 [extras]

src/InfiniteLinearAlgebra.jl

@@ -7,7 +7,7 @@ import Base: +, -, *, /, \, ^, OneTo, getindex, promote_op, _unsafe_getindex, pr
             show, getproperty, copy, map
 import Base.Broadcast: BroadcastStyle, Broadcasted, broadcasted
 
-import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, sublayout
+import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, sublayout, _qr
 import BandedMatrices: BandedMatrix, _BandedMatrix, AbstractBandedMatrix, bandeddata, bandwidths, BandedColumns, bandedcolumns,
                         _default_banded_broadcast
 import FillArrays: AbstractFill, getindex_value      
@@ -18,7 +18,7 @@ import LazyArrays: applybroadcaststyle, CachedArray, CachedMatrix, CachedVector,
                     factorize, sub_materialize, LazyLayout, LazyArrayStyle, layout_getindex,
                     @lazymul, applylayout, ApplyLayout, PaddedLayout, materialize!, zero!, MulAddStyle,
                     LazyArray, LazyMatrix, LazyVector
-import MatrixFactorizations: ql, ql!, QLPackedQ, getL, getR, reflector!, reflectorApply!, QL, QR, QRPackedQ
+import MatrixFactorizations: ql, ql!, _ql, QLPackedQ, getL, getR, reflector!, reflectorApply!, QL, QR, QRPackedQ
 
 import BlockArrays: AbstractBlockVecOrMat, sizes_from_blocks, _length, BlockedUnitRange
 

src/infql.jl

@@ -208,7 +208,7 @@ end
 # BlockTridiagonal
 ####
 
-function _ql(A::BlockTriPertToeplitz, d, e)
+function _blocktripert_ql(A, d, e)
     N = max(length(A.blocks.du.args[1])+1,length(A.blocks.d.args[1]),length(A.blocks.dl.args[1]))
     c,a,b = A[Block(N+1,N)],A[Block(N,N)],A[Block(N-1,N)]
     P,τ = blocktailiterate(c,a,b,d,e)
@@ -226,7 +226,7 @@ function _ql(A::BlockTriPertToeplitz, d, e)
             Vcat(F.τ,mortar(Fill(τ,∞)))), P[Block(1,1)], P[Block(1,2)]
 end
 
-ql(A::BlockTriPertToeplitz) = _ql(A, A[Block(2,3)], A[Block(3,3)])[1]
+ql(A::BlockTriPertToeplitz) = _blocktripert_ql(A, A[Block(2,3)], A[Block(3,3)])[1]
 
 ql(A::Adjoint{T,BlockTriPertToeplitz{T}}) where T = ql(BlockTridiagonal(A))
 

src/infqr.jl

@@ -16,27 +16,23 @@ function AdaptiveQRData(::AbstractAlmostBandedLayout, A::AbstractMatrix{T}) wher
     r = almostbandedrank(A)
     data = AlmostBandedMatrix{T}(undef,(2l+u+1,0),(l,l+u),r) # pad super
 
-    AdaptiveQRData(CachedArray(data,A), Vector{T}(), 0)
+    AdaptiveQRData(CachedArray(data,A,(0,0)), Vector{T}(), 0)
 end
 
 AdaptiveQRData(A::AbstractMatrix{T}) where T = AdaptiveQRData(MemoryLayout(typeof(A)), A)
 
 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);
+        l,u = bandwidths(F.data.data)
+        resizedata!(F.data,n+l,n+u);
         resize!(F.τ,n);
         ñ = F.ncols
         τ = view(F.τ,ñ+1:n);
         if l ≤ 0 
             zero!(τ)
         else
-            factors = view(F.data.data,ñ+1:n+l,ñ+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))
-            lmul!(Q',view(F.data.data,n̄:n+l,n+1:n+u+l))
+            factors = view(F.data.data,ñ+1:n+l,ñ+1:n+u);
+            _banded_qr!(factors, τ, n-ñ)
         end
         F.ncols = n
     end
@@ -45,20 +41,16 @@ 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)
+        l,u = almostbandwidths(F.data.data)
+        resizedata!(F.data,n+l,n+l+u);
+        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))
+            factors = view(F.data.data,ñ+1:n+l,ñ+1:n+l+u)
+            _almostbanded_qr!(factors, τ, n-ñ)
         end
         F.ncols = n
     end
@@ -118,7 +110,8 @@ function adaptiveqr(A)
     QR(AdaptiveQRFactors(data), AdaptiveQRTau(data))
 end
 
-_banded_qr(::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A)
+_qr(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A)
+_qr(::AbstractAlmostBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A)
 
 #########
 # lmul!
@@ -175,11 +168,12 @@ function lmul!(adjA::Adjoint{<:Any,<:QRPackedQ{<:Any,<:AdaptiveQRFactors}}, B::C
 
     @inbounds begin
         while first(jr) < nA
+            j = first(jr)
             cs = colsupport(A.factors, last(jr))
             cs_max = maximum(cs)
-            kr = first(jr):cs_max
+            kr = j:cs_max
             resizedata!(B, min(cs_max,mB))
-            if (first(jr) > sB && maximum(abs,view(B.data,colsupport(A.factors,first(jr)))) ≤ tol)
+            if (j > sB && maximum(abs,view(B.data,j:last(colsupport(A.factors,j)))) ≤ tol)
                 break
             end
             partialqr!(A.factors.data, min(cs_max,nA))

test/test_infqr.jl

@@ -143,9 +143,26 @@ import SemiseparableMatrices: AlmostBandedLayout, VcatAlmostBandedLayout
         @test A[1:10,1:10] isa AlmostBandedMatrix
         @test AlmostBandedMatrix(V) == Matrix(V) == A[1:10,1:10]
 
-        C = cache(A)
+        C = cache(A);
         @test C[1000,1000] ≡ 999.0
         F = adaptiveqr(A);
-        partialqr!(F.factors.data,10)
+        partialqr!(F.factors.data,2);
+        @test F.factors.data.data[1:3,1:5] ≈ qr(A[1:3,1:5]).factors
+        partialqr!(F.factors.data,3);
+        @test F.factors.data.data[1:4,1:6] ≈ qr(A[1:4,1:6]).factors
+
+        F = adaptiveqr(A);
+        partialqr!(F.factors.data,10);
+        @test F.factors[1:11,1:10] ≈ qr(A[1:11,1:10]).factors        
+        @test F.τ[1:10] ≈ qr(A[1:11,1:10]).τ        
+        partialqr!(F.factors.data,20);
+        @test F.factors[1:21,1:20] ≈ qr(A[1:21,1:20]).factors      
+
+        @test adaptiveqr(A).R[1:10,1:10] ≈ qr(A[1:11,1:10]).R
+        
+        @test qr(A) isa MatrixFactorizations.QR{Float64,<:InfiniteLinearAlgebra.AdaptiveQRFactors}
+        @test factorize(A) isa MatrixFactorizations.QR{Float64,<:InfiniteLinearAlgebra.AdaptiveQRFactors}
+
+        @test (adaptiveqr(A) \ [ℯ; zeros(∞)])[1:1000] ≈ (qr(A) \ [ℯ; zeros(∞)])[1:1000] ≈ (A \ [ℯ; zeros(∞)])[1:1000] ≈ [1/factorial(1.0k) for k=0:999]
     end 
 end