Special case diagonal Q in adaptive QR (#67)

* Special case diagonal Q in adaptive QR

* fix routine, add test

* tests pass on Julia v1.6

* import unitrange

* Update runtests.jl

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "InfiniteLinearAlgebra"
 uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
-version = "0.4.7"
+version = "0.4.8"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -18,7 +18,7 @@ SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 [compat]
 ArrayLayouts = "0.5.1"
 BandedMatrices = "0.16"
-BlockArrays = "0.14.2"
+BlockArrays = "0.14.3"
 BlockBandedMatrices = "0.10"
 FillArrays = "0.11"
 InfiniteArrays = "0.9.3"

src/InfiniteLinearAlgebra.jl

@@ -36,6 +36,12 @@ import BlockBandedMatrices: _BlockSkylineMatrix, _BandedMatrix, _BlockSkylineMat
 import SemiseparableMatrices: AbstractAlmostBandedLayout, _almostbanded_qr!
 
 
+if VERSION < v"1.6-"
+    oneto(n) = Base.OneTo(n)
+else
+    import Base: oneto, unitrange
+end
+
 # BroadcastStyle(::Type{<:BandedMatrix{<:Any,<:Any,<:OneToInf}}) = LazyArrayStyle{2}()
 
 function ArrayLayouts._power_by_squaring(_, ::NTuple{2,Infinity}, A::AbstractMatrix{T}, p::Integer) where T

src/banded/hessenbergq.jl

@@ -84,7 +84,7 @@ for Typ in (:UpperHessenbergQ, :LowerHessenbergQ)
 end
 
 size(Q::AbstractHessenbergQ, k::Integer) = length(Q.q)+1
-axes(Q::AbstractHessenbergQ, k::Integer) = Base.OneTo(length(Q.q)+1)
+axes(Q::AbstractHessenbergQ, k::Integer) = oneto(length(Q.q)+1)
 size(F::QLHessenberg, dim::Integer) = size(getfield(F, :factors), dim)
 size(F::QLHessenberg) = size(getfield(F, :factors))
 

src/blockbanded/blockbanded.jl

@@ -33,6 +33,8 @@ const OneToBlocks = BlockedUnitRange{OneToCumsum}
 axes(a::OneToInfBlocks) = (a,)
 axes(a::OneToBlocks) = (a,)
 
+unitrange(b::OneToInfBlocks) = first(b):∞
+
 
 function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{Ones{T,1,Tuple{OneToInfBlocks}},AbstractArray{V,N}}}) where {N,T,V}
     a,b = bc.args
@@ -78,4 +80,4 @@ BlockArrays._last(::BlockedUnitRange, ::OneToInf) = ∞
 ###
 
 _krontrav_axes(A::NTuple{N,OneToInf{Int}}, B::NTuple{N,OneToInf{Int}}) where N =
-     @. blockedrange(OneTo(length(A)))
+     @. blockedrange(oneto(length(A)))

src/infqr.jl

@@ -164,9 +164,9 @@ function materialize!(M::MatLmulVec{<:QRPackedQLayout{<:AdaptiveLayout},<:Padded
     A,B = M.A,M.B
     sB = B.datasize[1]
     partialqr!(A.factors.data,sB)
-    jr = Base.OneTo(sB)
+    jr = oneto(sB)
     m = maximum(colsupport(A,jr))
-    kr = Base.OneTo(m)
+    kr = oneto(m)
     resizedata!(B, m)
     b = paddeddata(B)
     lmul!(_view_QRPackedQ(A,kr,jr), b)
@@ -186,7 +186,12 @@ function materialize!(M::MatLmulVec{<:AdjQRPackedQLayout{<:AdaptiveLayout},<:Pad
     if mA != mB
         throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
     end
-    sB = length(B.data)
+    sB = B.datasize[1]
+    l,u = bandwidths(A.factors)
+    if l == 0 # diagonal special case
+        return B
+    end
+
     jr = 1:min(COLGROWTH,nA)
 
     @inbounds begin
@@ -287,7 +292,7 @@ end
 function ldiv!(R::UpperTriangular{<:Any,<:AdaptiveQRFactors}, B::CachedVector{<:Any,<:Any,<:Zeros{<:Any,1}})
     n = B.datasize[1]
     partialqr!(parent(R).data, n)
-    materialize!(Ldiv(UpperTriangular(view(parent(R).data.data.data,OneTo(n),OneTo(n))), view(B.data,OneTo(n))))
+    materialize!(Ldiv(UpperTriangular(view(parent(R).data.data.data,oneto(n),oneto(n))), view(B.data,oneto(n))))
     B
 end
 

test/runtests.jl

@@ -9,7 +9,7 @@ import BlockBandedMatrices: isblockbanded, _BlockBandedMatrix
 import MatrixFactorizations: QLPackedQ
 import BandedMatrices: bandeddata, _BandedMatrix, BandedStyle
 import LazyArrays: colsupport, ApplyStyle, MemoryLayout, ApplyLayout, LazyArrayStyle, arguments
-import InfiniteArrays: OneToInf
+import InfiniteArrays: OneToInf, oneto
 import LazyBandedMatrices: BroadcastBandedBlockBandedLayout, BroadcastBandedLayout
 
 
@@ -40,8 +40,8 @@ end
 
 @testset "∞-block arrays" begin
     @testset "fixed block size" begin
-        k = Base.OneTo.(Base.OneTo(∞))
-        n = Fill.(Base.OneTo(∞),Base.OneTo(∞))
+        k = Base.OneTo.(oneto(∞))
+        n = Fill.(oneto(∞),oneto(∞))
         @test broadcast(length,k) ≡ map(length,k) ≡ OneToInf()
         @test broadcast(length,n) ≡ map(length,n) ≡ OneToInf()
         b = mortar(Fill([1,2],∞))
@@ -52,7 +52,7 @@ end
     end
 
     @testset "1:∞ blocks" begin
-        a = blockedrange(Base.OneTo(∞))
+        a = blockedrange(oneto(∞))
         @test axes(a,1) == a
         o = Ones((a,))
         @test Base.BroadcastStyle(typeof(a)) isa LazyArrayStyle{1}
@@ -102,8 +102,8 @@ end
 
     @testset "triangle recurrences" begin
         @testset "n and k" begin
-            n = mortar(Fill.(Base.OneTo(∞),Base.OneTo(∞)))
-            k = mortar(Base.OneTo.(Base.OneTo(∞)))
+            n = mortar(Fill.(oneto(∞),oneto(∞)))
+            k = mortar(Base.OneTo.(oneto(∞)))
 
             @test n[Block(5)] ≡ layout_getindex(n, Block(5)) ≡ view(n,Block(5)) ≡ Fill(5,5)
             @test k[Block(5)] ≡ layout_getindex(k, Block(5)) ≡ view(k,Block(5)) ≡ Base.OneTo(5)
@@ -135,8 +135,8 @@ end
         end
 
         @testset "BlockHcat copyto!" begin
-            n = mortar(Fill.(Base.OneTo(∞),Base.OneTo(∞)))
-            k = mortar(Base.OneTo.(Base.OneTo(∞)))
+            n = mortar(Fill.(oneto(∞),oneto(∞)))
+            k = mortar(Base.OneTo.(oneto(∞)))
 
             a = b = c = 0.0
             dat = BlockHcat(
@@ -160,8 +160,8 @@ end
 
         @testset "BlockBanded" begin
             a = b = c = 0.0
-            n = mortar(Fill.(Base.OneTo(∞),Base.OneTo(∞)))
-            k = mortar(Base.OneTo.(Base.OneTo(∞)))
+            n = mortar(Fill.(oneto(∞),oneto(∞)))
+            k = mortar(Base.OneTo.(oneto(∞)))
             Dy = BlockBandedMatrices._BandedBlockBandedMatrix((k .+ (b+c))', axes(k,1), (-1,1), (-1,1))
             N = 100; 
             @test Dy[Block.(1:N), Block.(1:N)] == BlockBandedMatrices._BandedBlockBandedMatrix((k .+ (b+c))[Block.(1:N)]', axes(k,1)[Block.(1:N)], (-1,1), (-1,1))

test/test_infqr.jl

@@ -136,6 +136,12 @@ import SemiseparableMatrices: AlmostBandedLayout, VcatAlmostBandedLayout
             @test F.Q[1:10,1:10] == Eye(10)
             @test F.R[1:10,1:10] == A[1:10,1:10]
         end
+
+        @testset "diag special case" begin
+            A = _BandedMatrix((1:∞)', ∞, 0, 0)
+            b = [[1,2,3]; zeros(∞)]
+            @test A \ b == [ones(3); zeros(∞)]
+        end
     end
 
     @testset "almost-banded" begin