Thin/Wide QR, Actually call QR

Author: dlfivefifty

src/BlockBandedMatrices.jl

@@ -26,7 +26,7 @@ import LazyArrays: AbstractStridedLayout, ColumnMajor, @lazymul, MatMulMatAdd, M
                     triangularlayout, UpperTriangularLayout, TriangularLayout, MatMulVec, MatLdivVec,
                     triangulardata, subarraylayout, _copyto!, @lazyldiv, @lazylmul,
                     ArrayMulArrayStyle, AbstractColumnMajor, DenseColumnMajor, ColumnMajor,
-                    DiagonalLayout, MatMulMat, apply!
+                    DiagonalLayout, MatMulMat, apply!, materialize!
 
 import BlockArrays: BlockSizes, nblocks, blocksize, blockcheckbounds, global2blockindex,
                         Block, BlockSlice, getblock, unblock, setblock!, globalrange,

src/blockskylineqr.jl

@@ -13,31 +13,38 @@ _apply_ql!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayo
 apply_ql!(A, τ, B) = _apply_ql!(MemoryLayout(A), MemoryLayout(τ), MemoryLayout(B), A, τ, B)
 
 function qr!(A::BlockBandedMatrix)
-    bs = BlockSizes((cumulsizes(blocksizes(A),1),))
-    τ = PseudoBlockVector{Float64}(undef, bs)
     l,u = blockbandwidths(A)
-    N,M = nblocks(A)
-    for K = 1:N
-        KR = Block.(K:min(K+l,N))
+    M,N = nblocks(A)
+    bs = M < N ? BlockSizes((cumulsizes(blocksizes(A),1),)) : BlockSizes((cumulsizes(blocksizes(A),2),))
+    τ = PseudoBlockVector{Float64}(undef, bs)
+    for K = 1:min(N,M)
+        KR = Block.(K:min(K+l,M))
         V = view(A,KR,Block(K))
         t = view(τ,Block(K))
         qrf!(V,t)
-        for J = K+1:min(K+u,M)
+        for J = K+1:min(K+u,N)
             apply_qr!(V, t, view(A,KR,Block(J)))
         end
     end
     QR(A,τ.blocks)
 end
 
 function ql!(A::BlockBandedMatrix)
-    bs = BlockSizes((cumulsizes(blocksizes(A),1),))
-    τ = PseudoBlockVector{Float64}(undef, bs)
     l,u = blockbandwidths(A)
-    N,M = nblocks(A)
-    for K = N:-1:1
-        KR = Block.(max(K-u,1):K)
+    M,N = nblocks(A)
+
+    bs = if M < N
+        throw(ArgumentError("Wide block-QL not implented"))
+    else
+        BlockSizes((cumulsizes(blocksizes(A),2),))
+    end
+    τ = PseudoBlockVector{Float64}(undef, bs)
+    
+    for K = N:-1:max(N - M + 1,1)
+        μ = M+K-N
+        KR = Block.(max(K-u,1):μ)
         V = view(A,KR,Block(K))
-        t = view(τ,Block(K))
+        t = view(τ,Block(K-N+min(M,N)))
         qlf!(V,t)
         for J = K-1:-1:max(K-l,1)
             apply_ql!(V, t, view(A,KR,Block(J)))
@@ -58,7 +65,7 @@ function lmul!(adjQ::Adjoint{<:Any,<:QRPackedQ{<:Any,<:BlockSkylineMatrix}}, Bin
     bs = BlockSizes((cumulsizes(blocksizes(A),1),))
     τ  = PseudoBlockArray(Q.τ, bs)
     B = PseudoBlockArray(Bin, bs)
-    for K = 1:N
+    for K = 1:min(N,M)
         KR = Block.(K:min(K+l,N))
         V = view(A,KR,Block(K))
         t = view(τ,Block(K))
@@ -87,14 +94,23 @@ end
 
 # avoid LinearALgebra Strided obsession 
 
+for Typ in (:StridedVector, :StridedMatrix, :AbstractVector, :AbstractMatrix)
+    @eval function ldiv!(A::QR{<:Any,<:BlockSkylineMatrix}, B::$Typ)
+        lmul!(adjoint(A.Q), B)
+        apply!(\, UpperTriangular(A.factors), B)
+    end
+end
+
+
+
 function ldiv!(A::QL{<:Any,<:BlockSkylineMatrix}, B::AbstractVector)
     lmul!(adjoint(A.Q), B)
-    B .= Ldiv(LowerTriangular(A.factors), B)
+    apply!(\, LowerTriangular(A.factors), B)
 end
 
 function ldiv!(A::QL{<:Any,<:BlockSkylineMatrix}, B::AbstractMatrix)
     lmul!(adjoint(A.Q), B)
-    B .= Ldiv(LowerTriangular(A.factors), B)
+    apply!(\, LowerTriangular(A.factors), B)
 end
 
 \(A::AbstractBlockBandedMatrix, b::AbstractVecOrMat) = qr(A)\b  # use QR because LU would be a _mess_ to implement

src/triblockbanded.jl

@@ -82,32 +82,25 @@ function _matchingblocks_triangular_mul!(::Val{'L'}, UNIT, A, dest)
     dest
 end
 
-@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
-         M::MatMulVec{<:TriangularLayout{UPLO,UNIT,<:AbstractBlockBandedLayout},
+@inline function materialize!(M::MatMulVec{<:TriangularLayout{UPLO,UNIT,<:AbstractBlockBandedLayout},
                                    <:AbstractStridedLayout,T,T}) where {UPLO,UNIT,T<:BlasFloat}
     U,x = M.args
-    @boundscheck size(U,1) == size(dest,1) || throw(BoundsError())
+    @boundscheck size(U,1) == size(x,1) || throw(BoundsError())
     if hasmatchingblocks(U)
-        x ≡ dest || copyto!(dest, x)
-        _matchingblocks_triangular_mul!(Val(UPLO), Val(UNIT), triangulardata(U), dest)
+        _matchingblocks_triangular_mul!(Val(UPLO), Val(UNIT), triangulardata(U), x)
     else # use default
-        if x ≡ dest
-            materialize!(MulAdd(MemoryLayout(dest), BandedBlockBandedColumnMajor(), MemoryLayout(x),
-                                one(T), U, copy(x), zero(T), dest))
-        else
-            materialize!(MulAdd(MemoryLayout(dest), BandedBlockBandedColumnMajor(), MemoryLayout(x),
-                                one(T), U, x, zero(T), dest))
-        end
+        materialize!(MulAdd(MemoryLayout(dest), BandedBlockBandedColumnMajor(), MemoryLayout(x),
+                            one(T), U, copy(x), zero(T), dest))
     end
 end
 
 
 
-@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
-         M::MatLdivVec{<:TriangularLayout{'U',UNIT,<:AbstractBlockBandedLayout},
-                                   <:AbstractStridedLayout}) where {T,UNIT}
-    U,x = M.args
-    x ≡ dest || copyto!(dest, x)
+@inline function materialize!(M::MatLdivVec{<:TriangularLayout{'U',UNIT,<:AbstractBlockBandedLayout},
+                                   <:AbstractStridedLayout}) where UNIT
+    U,dest = M.args
+    T = eltype(dest)
+
     A = triangulardata(U)
     @assert hasmatchingblocks(A)
 
@@ -122,7 +115,7 @@ end
     for K = N:-1:1
         b_2 = view(b, Block(K))
         Ũ = _triangular_matrix(Val('U'), Val(UNIT), view(A, Block(K,K)))
-        b_2 .= Ldiv(Ũ, b_2)
+        apply!(\, Ũ, b_2)
 
         if K ≥ 2
             KR = blockcolstart(A, K):Block(K-1)
@@ -135,11 +128,10 @@ end
     dest
 end
 
-@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
-         M::MatLdivVec{<:TriangularLayout{'L',UNIT,<:AbstractBlockBandedLayout},
-                                   <:AbstractStridedLayout}) where {UNIT,T}
-    L,x = M.args
-    x ≡ dest || copyto!(dest, x)
+@inline function materialize!(M::MatLdivVec{<:TriangularLayout{'L',UNIT,<:AbstractBlockBandedLayout},
+                                   <:AbstractStridedLayout}) where UNIT
+    L,dest = M.args
+    T = eltype(dest)
     A = triangulardata(L)
     @assert hasmatchingblocks(A)
 

test/test_blockskylineqr.jl

@@ -1,49 +1,126 @@
 using BlockBandedMatrices, LinearAlgebra, MatrixFactorizations
 
-@testset "BlockBandedMatrix QR" begin
-    N = 5
-    A = BlockBandedMatrix{Float64}(undef, (1:N,1:N), (2,1))
-    A.data .= randn.()
-
-    F = qr(A)
-    @test F.factors ≈ qr(Matrix(A)).factors
-    @test F.τ ≈ LinearAlgebra.qrfactUnblocked!(Matrix(A)).τ
-
-    Q,R = F
-    Q̃,R̃ = qr(Matrix(A))
-    @test R ≈ R̃
-    @test Q ≈ Q̃
-    
-    b = randn(size(A,1))
-    @test Q'b ≈ Q̃'b
-    @test Q*b ≈ Q̃*b
-    
-    @test F\b ≈ ldiv!(F, copy(b)) ≈ Matrix(A)\b ≈ A\b
-end
+@testset "BlockBandedMatrix QR/QL" begin
+    @testset "Square QR" begin
+        N = 5
+        A = BlockBandedMatrix{Float64}(undef, (1:N,1:N), (2,1))
+        A.data .= randn.()
+
+        F = qr(A)
+        @test F.factors ≈ qr(Matrix(A)).factors
+        @test F.τ ≈ LinearAlgebra.qrfactUnblocked!(Matrix(A)).τ
+
+        Q,R = F
+        Q̃,R̃ = qr(Matrix(A))
+        @test R ≈ R̃
+        @test Q ≈ Q̃
+        
+        b = randn(size(A,1))
+        @test Q'b ≈ Q̃'b
+        @test Q*b ≈ Q̃*b
+        
+        @test F\b ≈ ldiv!(F, copy(b)) ≈ Matrix(A)\b ≈ A\b
+    end
+
+    @testset "Thin QR" begin
+        N = 5
+        A = BlockBandedMatrix{Float64}(undef, (1:N+1,1:N), (2,1))
+        A.data .= randn.()
+
+        F = qr(A)
+        @test F.factors ≈ qr(Matrix(A)).factors
+        @test F.τ ≈ LinearAlgebra.qrfactUnblocked!(Matrix(A)).τ
+
+        Q,R = F
+        Q̃,R̃ = qr(Matrix(A))
+        @test R ≈ R̃
+        @test Q ≈ Q̃
+        
+        b = randn(size(A,1))
+        @test Q'b ≈ Q̃'b
+        @test Q*b ≈ Q̃*b
+        
+        @test_broken F\b ≈ ldiv!(F, copy(b))[1:15] ≈ Matrix(A)\b ≈ A\b
+    end
+
+    @testset "Wide QR" begin
+        N = 5
+        A = BlockBandedMatrix{Float64}(undef, (1:N,1:N+1), (2,1))
+        A.data .= randn.()
+
+        F = qr(A)
+        @test F.factors ≈ qr(Matrix(A)).factors
+        @test F.τ ≈ LinearAlgebra.qrfactUnblocked!(Matrix(A)).τ
+
+        Q,R = F
+        Q̃,R̃ = qr(Matrix(A))
+        @test R ≈ R̃
+        @test Q ≈ Q̃
+        
+        b = randn(size(A,1))
+        @test Q'b ≈ Q̃'b
+        @test Q*b ≈ Q̃*b
+
+        @test_throws DimensionMismatch ldiv!(F, copy(b))
+        
+        @test_broken F\b ≈ ldiv!(F, copy(b)) ≈ Matrix(A)\b ≈ A\b
+    end
 
-@testset "BlockBandedMatrix QL" begin
-    N = 5
-    A = BlockBandedMatrix{Float64}(undef, (1:N,1:N), (2,1))
-    A.data .= randn.()
-
-    F = ql(A)
-    @test F.factors ≈ ql(Matrix(A)).factors
-    @test F.τ ≈ MatrixFactorizations.qlfactUnblocked!(Matrix(A)).τ
-
-    Q,L = F
-    Q̃,L̃ = ql(Matrix(A))
-    @test L ≈ L̃
-    @test Q ≈ Q̃
-    
-    b = randn(size(A,1))
-    @test Q'b ≈ Q̃'b
-    @test Q*b ≈ Q̃*b
-
-    b = randn(size(A,1))
-    @test F\b ≈ ldiv!(F, copy(b)) ≈ Matrix(A)\b ≈ A\b
+    @testset "Square QL" begin
+        N = 5
+        A = BlockBandedMatrix{Float64}(undef, (1:N,1:N), (2,1))
+        A.data .= randn.()
+
+        F = ql(A)
+        @test F.factors ≈ ql(Matrix(A)).factors
+        @test F.τ ≈ MatrixFactorizations.qlfactUnblocked!(Matrix(A)).τ
+
+        Q,L = F
+        Q̃,L̃ = ql(Matrix(A))
+        @test L ≈ L̃
+        @test Q ≈ Q̃
+        
+        b = randn(size(A,1))
+        @test Q'b ≈ Q̃'b
+        @test Q*b ≈ Q̃*b
+
+        b = randn(size(A,1))
+        @test F\b ≈ ldiv!(F, copy(b)) ≈ Matrix(A)\b ≈ A\b
+    end
+
+    @testset "Thin QL" begin
+        N = 5
+        A = BlockBandedMatrix{Float64}(undef, (1:N+1,1:N), (2,1))
+        A.data .= randn.()
+
+        # F = ql(A)
+        # @test F.factors ≈ ql(Matrix(A)).factors
+        # @test F.τ ≈ MatrixFactorizations.qlfactUnblocked!(Matrix(A)).τ
+
+        # Q,L = F
+        # Q̃,L̃ = ql(Matrix(A))
+        # @test L ≈ L̃
+        # @test Q ≈ Q̃
+        
+        # b = randn(size(A,1))
+        # @test Q'b ≈ Q̃'b
+        # @test Q*b ≈ Q̃*b
+
+        # b = randn(size(A,1))
+        # @test F\b ≈ ldiv!(F, copy(b)) ≈ Matrix(A)\b ≈ A\b
+    end
+
+    @testset "Wide QL" begin
+        N = 5
+        A = BlockBandedMatrix{Float64}(undef, (1:N,1:N+1), (2,1))
+        A.data .= randn.()
+
+        @test_throws ArgumentError ql(A)
+    end
 end
 
 
+
 # N = 500;
 # A = BlockBandedMatrix{Float64}(undef, (1:N,1:N), (2,1));
 # A.data .= randn.();