From bbaf3084fe60ad43f610851bee778dd504e95bf0 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 18 May 2024 08:06:05 +0100 Subject: [PATCH] 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 ba8705608dc28cc96145aa83f09978065de97cd6. * v3.0 * Increase coverage --- .github/workflows/ci.yml | 1 + Project.toml | 13 +- ext/MatrixFactorizationsBandedMatricesExt.jl | 187 +++++++++++++++++++ test/runtests.jl | 4 +- test/test_banded.jl | 142 ++++++++++++++ 5 files changed, 344 insertions(+), 3 deletions(-) create mode 100644 ext/MatrixFactorizationsBandedMatricesExt.jl create mode 100644 test/test_banded.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ffd1adf..d1d413b 100644 --- a/.github/workflows/ci.yml +++ b/.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 diff --git a/Project.toml b/Project.toml index 34fcd23..419d07b 100644 --- a/Project.toml +++ b/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"] diff --git a/ext/MatrixFactorizationsBandedMatricesExt.jl b/ext/MatrixFactorizationsBandedMatricesExt.jl new file mode 100644 index 0000000..c8c8a0d --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 07b0493..93fc0dd 100644 --- a/test/runtests.jl +++ b/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") \ No newline at end of file +include("test_reversecholesky.jl") + +include("test_banded.jl") \ No newline at end of file diff --git a/test/test_banded.jl b/test/test_banded.jl new file mode 100644 index 0000000..ea17704 --- /dev/null +++ b/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 \ No newline at end of file