Faster QR (#75)

* Faster QR

* Update Project.toml

* Update Project.toml

Author: dlfivefifty

src/AbstractBlockBandedMatrix.jl

@@ -79,9 +79,10 @@ blockbanded_blockrowstart(A, i::BlockRange) = blockbanded_blockrowstart(A, minim
 blockbanded_blockcolstop(A, i::BlockRange) = blockbanded_blockcolstop(A, maximum(i))
 blockbanded_blockrowstop(A, i::BlockRange) = blockbanded_blockrowstop(A, maximum(i))
 
+const AllBlockBandedLayout{UPLO,UNIT} = Union{AbstractBlockBandedLayout,TriangularLayout{UPLO,UNIT,<:AbstractBlockBandedLayout}}
 
-@inline blockcolsupport(::AbstractBlockBandedLayout, A, i) = blockbanded_blockcolstart(A,i):blockbanded_blockcolstop(A,i)
-@inline blockrowsupport(::AbstractBlockBandedLayout, A, i) = blockbanded_blockrowstart(A,i):blockbanded_blockrowstop(A,i)
+@inline blockcolsupport(::AllBlockBandedLayout, A, i) = blockbanded_blockcolstart(A,i):blockbanded_blockcolstop(A,i)
+@inline blockrowsupport(::AllBlockBandedLayout, A, i) = blockbanded_blockrowstart(A,i):blockbanded_blockrowstop(A,i)
 
 
 # length of i-the column/row

src/BlockBandedMatrices.jl

@@ -32,7 +32,7 @@ import BlockArrays: blocksize, blockcheckbounds, BlockedUnitRange, blockisequal,
                         Block, BlockSlice, getblock, unblock, setblock!, block, blockindex,
                         _blocklengths2blocklasts, BlockIndexRange, sizes_from_blocks, BlockSlice1,
                         blockcolsupport, blockrowsupport, blockcolstart, blockcolstop, blockrowstart, blockrowstop,
-                        AbstractBlockLayout, BlockLayout, blocks
+                        AbstractBlockLayout, BlockLayout, blocks, hasmatchingblocks
 
 import BandedMatrices: isbanded, bandwidths, bandwidth, banded_getindex, colrange,
                         inbands_setindex!, inbands_getindex, banded_setindex!,

src/blockskylineqr.jl

@@ -109,15 +109,38 @@ function materialize!(Mul::MatLmulVec{<:AdjQLPackedQLayout{<:AbstractBlockBanded
     Mul.B
 end
 
+function materialize!(Mul::MatLmulMat{<:AdjQRPackedQLayout{<:AbstractBlockBandedLayout}})
+    adjQ,Bin = Mul.A,Mul.B
+    Q = parent(adjQ)
+    A = Q.factors
+    l,u = blockbandwidths(A)
+    N,M = blocksize(A)
+    # impose block structure
+    ax1 = (axes(A,1),)
+    τ  = PseudoBlockArray(Q.τ, ax1)
+    B = PseudoBlockArray(Bin, (axes(A,1),axes(Bin,2)))
+    for K = 1:min(N,M), J = 1:blocksize(Bin,2)
+        KR = Block.(K:min(K+l,N))
+        V = view(A,KR,Block(K))
+        t = view(τ,Block(K))
+        apply_qr!(V, t, view(B,KR,Block(J)))
+    end
+    Bin
+end
+
 # avoid LinearALgebra Strided obsession 
 
 for Typ in (:StridedVector, :StridedMatrix, :AbstractVector, :AbstractMatrix, :LayoutMatrix)
     @eval function ldiv!(A::QR{<:Any,<:BlockSkylineMatrix}, B::$Typ)
         lmul!(adjoint(A.Q), B)
         M,N = blocksize(A.factors)
-        MN = min(M,N)
-        V = view(A.factors,Block.(1:MN), Block.(1:MN))
-        materialize!(Ldiv(UpperTriangular(V), view(B,1:size(V,1),:)))
+        if M == N
+            materialize!(Ldiv(UpperTriangular(A.factors), view(B,1:size(A.factors,1),:)))
+        else
+            MN = min(M,N)
+            V = view(A.factors,Block.(1:MN), Block.(1:MN))
+            materialize!(Ldiv(UpperTriangular(V), view(B,1:size(V,1),:)))
+        end
         B
     end
 end

src/triblockbanded.jl

@@ -10,7 +10,7 @@ isbandedblockbanded(A::AbstractTriangular) =
 function blockbandwidths(A::Union{UpperTriangular,UnitUpperTriangular}) 
     P = parent(A)
     if hasmatchingblocks(P)
-        (min(0,blockbandwidths(P,1)), blockbandwidth(P,2))
+        (min(0,blockbandwidth(P,1)), blockbandwidth(P,2))
     else
         blockbandwidths(P)
     end

test/test_blockskylineqr.jl

@@ -1,4 +1,5 @@
-using BlockBandedMatrices, LinearAlgebra, MatrixFactorizations
+using BlockBandedMatrices, BlockArrays, LinearAlgebra, MatrixFactorizations, Test
+import BlockBandedMatrices: blockcolsupport
 
 @testset "BlockBandedMatrix QR/QL" begin
     @testset "Square QR" begin
@@ -118,7 +119,7 @@ using BlockBandedMatrices, LinearAlgebra, MatrixFactorizations
         N = 5
         A = BlockBandedMatrix{Float64}(undef, 1:N,1:N+1, (2,1))
         A.data .= randn.()
-        
+
         @test_throws ArgumentError ql(A)
     end
 
@@ -148,6 +149,34 @@ using BlockBandedMatrices, LinearAlgebra, MatrixFactorizations
             end
         end
     end
+
+    @testset "Fast QR \\" begin
+        A = [1. 2; 3 4]; A = A + A'
+        B = [5. 6; 7 8]
+        N = 10_000;
+        T = BlockBandedMatrix(mortar(Tridiagonal(fill(Matrix(B'),N-1), fill(zeros(2,2),N), fill(B,N-1))))
+        z = 1+2im
+        @time F = qr(T - z*I);
+        @test blockbandwidths(UpperTriangular(F.factors)) == (0,2)
+        @test blockcolsupport(UpperTriangular(F.factors),Block(4)) == Block.(2:4)
+        V = view(F.factors,Block.(1:N), Block.(1:N))
+        @test blocksize(V) == (N,N)
+        @test @allocated(blocksize(V)) ≤ 40
+        @test blockbandwidths(V) == (1,2)
+        @test @allocated(blockbandwidths(V)) ≤ 40
+        R = UpperTriangular(V)
+        @test axes(R) == axes(V) == axes(F.factors)
+        @test blocksize(R) == (N,N)
+        @test @allocated(blocksize(R)) ≤ 40
+        @test blockbandwidths(R) == (0,2)
+        @test @allocated(blockbandwidths(R)) ≤ 40
+        @test blockcolsupport(R,Block(4)) == Block.(2:4)
+        @test @allocated(blockcolsupport(R,Block(4))) ≤ 40
+
+        b = [1; zeros(size(T,1)-1)]
+        B = [Matrix(I,2,2); zeros(size(T,1)-2,2)]
+        @test ((T - z*I)\b)[1] ≈ (F\b)[1] ≈ (F \ B)[1,1] ≈ ((T - z*I)\B)[1,1] ≈ -0.1309123477325813 + 0.28471699370329884im 
+    end
 end