v3.0 (#65)

* BandedMatrices Extension (#64)

* BandedMatrices Extension

* Update Project.toml

* Add BandedMatrix extension and tests

* Update MatrixFactorizationsBandedMatricesExt.jl

* fix tests

* v2.3

* Update ci.yml

* Revert "v2.3"

This reverts commit ba87056.

* v3.0

* Increase coverage

.github/workflows/ci.yml

@@ -49,4 +49,5 @@ jobs:
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v4
         with:
+          token: ${{ secrets.CODECOV_TOKEN }}
           file: lcov.info

Project.toml

@@ -1,19 +1,28 @@
 name = "MatrixFactorizations"
 uuid = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
-version = "2.2"
+version = "3.0.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
+[weakdeps]
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+
+[extensions]
+MatrixFactorizationsBandedMatricesExt = "BandedMatrices"
+
 [compat]
 ArrayLayouts = "1.9.2"
+BandedMatrices = "1.6"
 julia = "1.9"
 
 [extras]
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["BandedMatrices", "Test"]

ext/MatrixFactorizationsBandedMatricesExt.jl

@@ -0,0 +1,187 @@
+module MatrixFactorizationsBandedMatricesExt
+using BandedMatrices, MatrixFactorizations, LinearAlgebra
+using MatrixFactorizations.ArrayLayouts
+import MatrixFactorizations: ql, ql!, QLPackedQLayout, AdjQLPackedQLayout, QL
+import ArrayLayouts: materialize!, reflector!, reflectorApply!
+import LinearAlgebra: ldiv!
+using BandedMatrices: bandeddata
+using Base: require_one_based_indexing
+###
+# QL
+###
+
+
+ql(A::BandedMatrix{T}) where T = ql!(BandedMatrix{float(T)}(A, (max(bandwidth(A,1),bandwidth(A,1)+bandwidth(A,2)+size(A,1)-size(A,2)),bandwidth(A,2))))
+
+ql!(A::BandedMatrix) = banded_ql!(A)
+
+function banded_ql!(L::BandedMatrix{T}) where T
+    D = bandeddata(L)
+    l,u = bandwidths(L)
+    ν = l+u+1
+    m,n=size(L)
+    τ = zeros(T, min(m,n))
+
+    for k = n:-1:max((n - m + 1 + (T<:Real)),1)
+        μ = m+k-n
+        x = view(D,u+1+μ-k:-1:max(1,u-k+2), k)
+        τk = reflector!(x)
+        τ[k-n+min(m,n)] = τk
+        N = length(x)
+        for j = max(1,μ-l):k-1
+            reflectorApply!(x, τk, view(D, u+1+μ-j:-1:u+2+μ-j-N,j))
+        end
+    end
+    QL(L, τ)
+end
+
+function materialize!(M::Lmul{<:QLPackedQLayout{<:AbstractBandedLayout}})
+    A,B = M.A,M.B
+    require_one_based_indexing(B)
+    mA, nA = size(A.factors)
+    mB, nB = size(B,1), size(B,2)
+    if mA != mB
+        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
+    end
+    Afactors = A.factors
+    l,u = bandwidths(Afactors)
+    D = bandeddata(Afactors)
+    for k = max(nA - mA + 1,1):nA
+        μ = mA+k-nA
+        for j = 1:nB
+            vBj = B[μ,j]
+            for i = max(1,k-u):μ-1
+                vBj += conj(D[i-k+u+1,k])*B[i,j]
+            end
+            vBj = A.τ[k-nA+min(mA,nA)]*vBj
+            B[μ,j] -= vBj
+            for i = max(1,k-u):μ-1
+                B[i,j] -= D[i-k+u+1,k]*vBj
+            end
+        end
+    end
+    B
+end
+
+function materialize!(M::Lmul{<:AdjQLPackedQLayout{<:AbstractBandedLayout}})
+    adjA,B = M.A,M.B
+    require_one_based_indexing(B)
+    A = parent(adjA)
+    mA, nA = size(A.factors)
+    mB, nB = size(B,1), size(B,2)
+    if mA != mB
+        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
+    end
+    Afactors = A.factors
+    l,u = bandwidths(Afactors)
+    D = bandeddata(Afactors)
+    @inbounds begin
+        for k = nA:-1:max(nA - mA + 1,1)
+            μ = mA+k-nA
+            for j = 1:nB
+                vBj = B[μ,j]
+                for i = max(1,k-u):μ-1
+                    vBj += conj(D[i-k+u+1,k])*B[i,j]
+                end
+                vBj = conj(A.τ[k-nA+min(mA,nA)])*vBj
+                B[μ,j] -= vBj
+                for i = max(1,k-u):μ-1
+                    B[i,j] -= D[i-k+u+1,k]*vBj
+                end
+            end
+        end
+    end
+    B
+end
+
+### QBc/QcBc
+function materialize!(M::Rmul{<:Any,<:QLPackedQLayout{<:AbstractBandedLayout}})
+    A,Q = M.A,M.B
+    mQ, nQ = size(Q.factors)
+    mA, nA = size(A,1), size(A,2)
+    if nA != mQ
+        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
+    end
+    Qfactors = Q.factors
+    l,u = bandwidths(Qfactors)
+    D = Qfactors.data
+    @inbounds begin
+        for k = nQ:-1:max(nQ - mQ + 1,1)
+            μ = mQ+k-nQ
+            for i = 1:mA
+                vAi = A[i,μ]
+                for j = max(1,k-u):μ-1
+                    vAi += A[i,j]*D[j-k+u+1,k]
+                end
+                vAi = vAi*Q.τ[k-nQ+min(mQ,nQ)]
+                A[i,μ] -= vAi
+                for j = max(1,k-u):μ-1
+                    A[i,j] -= vAi*conj(D[j-k+u+1,k])
+                end
+            end
+        end
+    end
+    A
+end
+
+### AQc
+function materialize!(M::Rmul{<:Any,<:AdjQLPackedQLayout{<:AbstractBandedLayout}})
+    A,adjQ = M.A,M.B
+    Q = parent(adjQ)
+    mQ, nQ = size(Q.factors)
+    mA, nA = size(A,1), size(A,2)
+    if nA != mQ
+        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
+    end
+    Qfactors = Q.factors
+    l,u = bandwidths(Qfactors)
+    D = Qfactors.data
+    @inbounds begin
+        for k = max(nQ - mQ + 1,1):nQ
+            μ = mQ+k-nQ
+            for i = 1:mA
+                vAi = A[i,μ]
+                for j = max(1,k-u):μ-1
+                    vAi += A[i,j]*D[j-k+u+1,k]
+                end
+                vAi = vAi*conj(Q.τ[k-nQ+min(mQ,nQ)])
+                A[i,μ] -= vAi
+                for j = max(1,k-u):μ-1
+                    A[i,j] -= vAi*conj(D[j-k+u+1,k])
+                end
+            end
+        end
+    end
+    A
+end
+
+
+
+
+function _banded_widerect_ldiv!(A::QL, B)
+    error("Not implemented")
+end
+function _banded_longrect_ldiv!(A::QL, B)
+    error("Not implemented")
+end
+function _banded_square_ldiv!(A::QL, B)
+    L = A.factors
+    lmul!(adjoint(A.Q), B)
+    B .= Ldiv(LowerTriangular(L), B)
+    B
+end
+
+for Typ in (:StridedVector, :StridedMatrix, :AbstractVecOrMat)
+    @eval function ldiv!(A::QL{T,<:BandedMatrix}, B::$Typ{T}) where T
+        m, n = size(A)
+        if m == n
+            _banded_square_ldiv!(A, B)
+        elseif n > m
+            _banded_widerect_ldiv!(A, B)
+        else
+            _banded_longrect_ldiv!(A, B)
+        end
+    end
+end
+
+end # module

test/runtests.jl

@@ -29,4 +29,6 @@ include("test_qr.jl")
 include("test_ql.jl")
 include("test_rq.jl")
 include("test_polar.jl")
-include("test_reversecholesky.jl")
+include("test_reversecholesky.jl")
+
+include("test_banded.jl")

test/test_banded.jl

@@ -0,0 +1,142 @@
+module BandedMatrixFactorizationTests
+using MatrixFactorizations, LinearAlgebra, BandedMatrices, Test
+
+@testset "QL tests" begin
+    for T in (Float64,ComplexF64,Float32,ComplexF32)
+        A=brand(T,10,10,3,2)
+        Q,L=ql(A)
+        @test ql(A).factors ≈ ql!(Matrix(A)).factors
+        @test ql(A).τ ≈ ql!(Matrix(A)).τ
+        @test Matrix(Q)*Matrix(L) ≈ A
+        b=rand(T,10)
+        @test mul!(similar(b),Q,mul!(similar(b),Q',b)) ≈ b
+        for j=1:size(A,2)
+            @test Q' * A[:,j] ≈ L[:,j]
+        end
+
+        A=brand(T,14,10,3,2)
+        Q,L=ql(A)
+        @test ql(A).factors ≈ ql!(Matrix(A)).factors
+        @test ql(A).τ ≈ ql!(Matrix(A)).τ
+        @test_broken Matrix(Q)*Matrix(L) ≈ A
+
+        for k=1:size(A,1),j=1:size(A,2)
+            @test Q[k,j] ≈ Matrix(Q)[k,j]
+        end
+
+        A=brand(T,10,14,3,2)
+        Q,L=ql(A)
+        @test ql(A).factors ≈ ql!(Matrix(A)).factors
+        @test ql(A).τ ≈ ql!(Matrix(A)).τ
+        @test Matrix(Q)*Matrix(L) ≈ A
+
+        for k=1:size(Q,1),j=1:size(Q,2)
+            @test Q[k,j] ≈ Matrix(Q)[k,j]
+        end
+
+        A=brand(T,10,14,3,6)
+        Q,L=ql(A)
+        @test ql(A).factors ≈ ql!(Matrix(A)).factors
+        @test ql(A).τ ≈ ql!(Matrix(A)).τ
+        @test Matrix(Q)*Matrix(L) ≈ A
+
+        for k=1:size(Q,1),j=1:size(Q,2)
+            @test Q[k,j] ≈ Matrix(Q)[k,j]
+        end
+
+        A=brand(T,100,100,3,4)
+        @test ql(A).factors ≈ ql!(Matrix(A)).factors
+        @test ql(A).τ ≈ ql!(Matrix(A)).τ
+        b=rand(T,100)
+        @test ql(A)\b ≈ Matrix(A)\b
+        b=rand(T,100,2)
+        @test ql(A)\b ≈ Matrix(A)\b
+        @test_throws DimensionMismatch ql(A) \ randn(3)
+        @test_throws DimensionMismatch ql(A).Q'randn(3)
+
+        A=brand(T,102,100,3,4)
+        @test ql(A).factors ≈ ql!(Matrix(A)).factors
+        @test ql(A).τ ≈ ql!(Matrix(A)).τ
+        b=rand(T,102)
+        @test_broken ql(A)\b ≈ Matrix(A)\b
+        b=rand(T,102,2)
+        @test_broken ql(A)\b ≈ Matrix(A)\b
+        @test_throws DimensionMismatch ql(A) \ randn(3)
+        @test_throws DimensionMismatch ql(A).Q'randn(3)
+
+        A=brand(T,100,102,3,4)
+        @test ql(A).factors ≈ ql!(Matrix(A)).factors
+        @test ql(A).τ ≈ ql!(Matrix(A)).τ
+        b=rand(T,100)
+        @test_broken ql(A)\b ≈ Matrix(A)\b
+
+        A = LinearAlgebra.Tridiagonal(randn(T,99), randn(T,100), randn(T,99))
+        @test ql(A).factors ≈ ql!(Matrix(A)).factors
+        @test ql(A).τ ≈ ql!(Matrix(A)).τ
+        b=rand(T,100)
+        @test ql(A)\b ≈ Matrix(A)\b
+        b=rand(T,100,2)
+        @test ql(A)\b ≈ Matrix(A)\b
+        @test_throws DimensionMismatch ql(A) \ randn(3)
+        @test_throws DimensionMismatch ql(A).Q'randn(3)
+    end
+
+    @testset "lmul!/rmul!" begin
+        for T in (Float32, Float64, ComplexF32, ComplexF64)
+            A = brand(T,100,100,3,4)
+            Q,R = qr(A)
+            x = randn(T,100)
+            b = randn(T,100,2)
+            @test lmul!(Q, copy(x)) ≈ Matrix(Q)*x
+            @test lmul!(Q, copy(b)) ≈ Matrix(Q)*b
+            @test lmul!(Q', copy(x)) ≈ Matrix(Q)'*x
+            @test lmul!(Q', copy(b)) ≈ Matrix(Q)'*b
+            c = randn(T,2,100)
+            @test rmul!(copy(c), Q) ≈ c*Matrix(Q)
+            @test rmul!(copy(c), Q') ≈ c*Matrix(Q')
+
+            A = brand(T,100,100,3,4)
+            Q,L = ql(A)
+            x = randn(T,100)
+            b = randn(T,100,2)
+            @test lmul!(Q, copy(x)) ≈ Matrix(Q)*x
+            @test lmul!(Q, copy(b)) ≈ Matrix(Q)*b
+            @test lmul!(Q', copy(x)) ≈ Matrix(Q)'*x
+            @test lmul!(Q', copy(b)) ≈ Matrix(Q)'*b
+            c = randn(T,2,100)
+            @test rmul!(copy(c), Q) ≈ c*Matrix(Q)
+            @test rmul!(copy(c), Q') ≈ c*Matrix(Q')
+        end
+    end
+
+    @testset "Mixed types" begin
+        A=brand(10,10,3,2)
+        b=rand(ComplexF64,10)
+        Q,L=ql(A)
+        @test L\(Q'*b) ≈ ql(A)\b ≈ Matrix(A)\b
+        @test Q*L ≈ A
+
+        A=brand(ComplexF64,10,10,3,2)
+        b=rand(10)
+        Q,L=ql(A)
+        @test Q*L ≈ A
+        @test L\(Q'*b) ≈ ql(A)\b ≈ Matrix(A)\b
+
+        A = BandedMatrix{Int}(undef, (2,1), (4,4))
+        A.data .= 1:length(A.data)
+        Q, L = ql(A)
+        @test_broken Q*L ≈ A
+    end
+
+    @testset "bounds" begin
+        A = brand(100,100,3,4)
+        @test_throws DimensionMismatch ql(A) \ randn(50)
+        @test_throws DimensionMismatch ldiv!(ql(A), randn(50))
+        Q = ql(A).Q
+        @test_throws DimensionMismatch lmul!(Q, randn(50))
+        @test_throws DimensionMismatch lmul!(Q', randn(50))
+        @test_throws DimensionMismatch rmul!(randn(50)', Q)
+        @test_throws DimensionMismatch rmul!(randn(50)', Q')
+    end
+end
+end