Skip to content

Commit bbaf308

Browse files
dlfivefiftydlfivefifty
authored
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
1 parent cf5824f commit bbaf308

File tree

5 files changed

+344
-3
lines changed

5 files changed

+344
-3
lines changed

‎.github/workflows/ci.yml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,4 +49,5 @@ jobs:
4949
- uses: julia-actions/julia-processcoverage@v1
5050
- uses: codecov/codecov-action@v4
5151
with:
52+
token: ${{ secrets.CODECOV_TOKEN }}
5253
file: lcov.info

‎Project.toml‎

Lines changed: 11 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,19 +1,28 @@
11
name = "MatrixFactorizations"
22
uuid = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
3-
version = "2.2"
3+
version = "3.0.0"
44

55
[deps]
66
ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
7+
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
78
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
89
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
910
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
1011

12+
[weakdeps]
13+
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
14+
15+
[extensions]
16+
MatrixFactorizationsBandedMatricesExt = "BandedMatrices"
17+
1118
[compat]
1219
ArrayLayouts = "1.9.2"
20+
BandedMatrices = "1.6"
1321
julia = "1.9"
1422

1523
[extras]
24+
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
1625
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
1726

1827
[targets]
19-
test = ["Test"]
28+
test = ["BandedMatrices", "Test"]

‎ext/MatrixFactorizationsBandedMatricesExt.jl‎

Lines changed: 187 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,187 @@
1+
module MatrixFactorizationsBandedMatricesExt
2+
using BandedMatrices, MatrixFactorizations, LinearAlgebra
3+
using MatrixFactorizations.ArrayLayouts
4+
import MatrixFactorizations: ql, ql!, QLPackedQLayout, AdjQLPackedQLayout, QL
5+
import ArrayLayouts: materialize!, reflector!, reflectorApply!
6+
import LinearAlgebra: ldiv!
7+
using BandedMatrices: bandeddata
8+
using Base: require_one_based_indexing
9+
###
10+
# QL
11+
###
12+
13+
14+
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))))
15+
16+
ql!(A::BandedMatrix) = banded_ql!(A)
17+
18+
function banded_ql!(L::BandedMatrix{T}) where T
19+
D = bandeddata(L)
20+
l,u = bandwidths(L)
21+
ν = l+u+1
22+
m,n=size(L)
23+
τ = zeros(T, min(m,n))
24+
25+
for k = n:-1:max((n - m + 1 + (T<:Real)),1)
26+
μ = m+k-n
27+
x = view(D,u+1+μ-k:-1:max(1,u-k+2), k)
28+
τk = reflector!(x)
29+
τ[k-n+min(m,n)] = τk
30+
N = length(x)
31+
for j = max(1,μ-l):k-1
32+
reflectorApply!(x, τk, view(D, u+1+μ-j:-1:u+2+μ-j-N,j))
33+
end
34+
end
35+
QL(L, τ)
36+
end
37+
38+
function materialize!(M::Lmul{<:QLPackedQLayout{<:AbstractBandedLayout}})
39+
A,B = M.A,M.B
40+
require_one_based_indexing(B)
41+
mA, nA = size(A.factors)
42+
mB, nB = size(B,1), size(B,2)
43+
if mA != mB
44+
throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
45+
end
46+
Afactors = A.factors
47+
l,u = bandwidths(Afactors)
48+
D = bandeddata(Afactors)
49+
for k = max(nA - mA + 1,1):nA
50+
μ = mA+k-nA
51+
for j = 1:nB
52+
vBj = B[μ,j]
53+
for i = max(1,k-u):μ-1
54+
vBj += conj(D[i-k+u+1,k])*B[i,j]
55+
end
56+
vBj = A.τ[k-nA+min(mA,nA)]*vBj
57+
B[μ,j] -= vBj
58+
for i = max(1,k-u):μ-1
59+
B[i,j] -= D[i-k+u+1,k]*vBj
60+
end
61+
end
62+
end
63+
B
64+
end
65+
66+
function materialize!(M::Lmul{<:AdjQLPackedQLayout{<:AbstractBandedLayout}})
67+
adjA,B = M.A,M.B
68+
require_one_based_indexing(B)
69+
A = parent(adjA)
70+
mA, nA = size(A.factors)
71+
mB, nB = size(B,1), size(B,2)
72+
if mA != mB
73+
throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
74+
end
75+
Afactors = A.factors
76+
l,u = bandwidths(Afactors)
77+
D = bandeddata(Afactors)
78+
@inbounds begin
79+
for k = nA:-1:max(nA - mA + 1,1)
80+
μ = mA+k-nA
81+
for j = 1:nB
82+
vBj = B[μ,j]
83+
for i = max(1,k-u):μ-1
84+
vBj += conj(D[i-k+u+1,k])*B[i,j]
85+
end
86+
vBj = conj(A.τ[k-nA+min(mA,nA)])*vBj
87+
B[μ,j] -= vBj
88+
for i = max(1,k-u):μ-1
89+
B[i,j] -= D[i-k+u+1,k]*vBj
90+
end
91+
end
92+
end
93+
end
94+
B
95+
end
96+
97+
### QBc/QcBc
98+
function materialize!(M::Rmul{<:Any,<:QLPackedQLayout{<:AbstractBandedLayout}})
99+
A,Q = M.A,M.B
100+
mQ, nQ = size(Q.factors)
101+
mA, nA = size(A,1), size(A,2)
102+
if nA != mQ
103+
throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
104+
end
105+
Qfactors = Q.factors
106+
l,u = bandwidths(Qfactors)
107+
D = Qfactors.data
108+
@inbounds begin
109+
for k = nQ:-1:max(nQ - mQ + 1,1)
110+
μ = mQ+k-nQ
111+
for i = 1:mA
112+
vAi = A[i,μ]
113+
for j = max(1,k-u):μ-1
114+
vAi += A[i,j]*D[j-k+u+1,k]
115+
end
116+
vAi = vAi*Q.τ[k-nQ+min(mQ,nQ)]
117+
A[i,μ] -= vAi
118+
for j = max(1,k-u):μ-1
119+
A[i,j] -= vAi*conj(D[j-k+u+1,k])
120+
end
121+
end
122+
end
123+
end
124+
A
125+
end
126+
127+
### AQc
128+
function materialize!(M::Rmul{<:Any,<:AdjQLPackedQLayout{<:AbstractBandedLayout}})
129+
A,adjQ = M.A,M.B
130+
Q = parent(adjQ)
131+
mQ, nQ = size(Q.factors)
132+
mA, nA = size(A,1), size(A,2)
133+
if nA != mQ
134+
throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
135+
end
136+
Qfactors = Q.factors
137+
l,u = bandwidths(Qfactors)
138+
D = Qfactors.data
139+
@inbounds begin
140+
for k = max(nQ - mQ + 1,1):nQ
141+
μ = mQ+k-nQ
142+
for i = 1:mA
143+
vAi = A[i,μ]
144+
for j = max(1,k-u):μ-1
145+
vAi += A[i,j]*D[j-k+u+1,k]
146+
end
147+
vAi = vAi*conj(Q.τ[k-nQ+min(mQ,nQ)])
148+
A[i,μ] -= vAi
149+
for j = max(1,k-u):μ-1
150+
A[i,j] -= vAi*conj(D[j-k+u+1,k])
151+
end
152+
end
153+
end
154+
end
155+
A
156+
end
157+
158+
159+
160+
161+
function _banded_widerect_ldiv!(A::QL, B)
162+
error("Not implemented")
163+
end
164+
function _banded_longrect_ldiv!(A::QL, B)
165+
error("Not implemented")
166+
end
167+
function _banded_square_ldiv!(A::QL, B)
168+
L = A.factors
169+
lmul!(adjoint(A.Q), B)
170+
B .= Ldiv(LowerTriangular(L), B)
171+
B
172+
end
173+
174+
for Typ in (:StridedVector, :StridedMatrix, :AbstractVecOrMat)
175+
@eval function ldiv!(A::QL{T,<:BandedMatrix}, B::$Typ{T}) where T
176+
m, n = size(A)
177+
if m == n
178+
_banded_square_ldiv!(A, B)
179+
elseif n > m
180+
_banded_widerect_ldiv!(A, B)
181+
else
182+
_banded_longrect_ldiv!(A, B)
183+
end
184+
end
185+
end
186+
187+
end # module

‎test/runtests.jl‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,4 +29,6 @@ include("test_qr.jl")
2929
include("test_ql.jl")
3030
include("test_rq.jl")
3131
include("test_polar.jl")
32-
include("test_reversecholesky.jl")
32+
include("test_reversecholesky.jl")
33+
34+
include("test_banded.jl")

‎test/test_banded.jl‎

Lines changed: 142 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,142 @@
1+
module BandedMatrixFactorizationTests
2+
using MatrixFactorizations, LinearAlgebra, BandedMatrices, Test
3+
4+
@testset "QL tests" begin
5+
for T in (Float64,ComplexF64,Float32,ComplexF32)
6+
A=brand(T,10,10,3,2)
7+
Q,L=ql(A)
8+
@test ql(A).factors ≈ ql!(Matrix(A)).factors
9+
@test ql(A).τ ≈ ql!(Matrix(A)).τ
10+
@test Matrix(Q)*Matrix(L) ≈ A
11+
b=rand(T,10)
12+
@test mul!(similar(b),Q,mul!(similar(b),Q',b)) ≈ b
13+
for j=1:size(A,2)
14+
@test Q' * A[:,j] ≈ L[:,j]
15+
end
16+
17+
A=brand(T,14,10,3,2)
18+
Q,L=ql(A)
19+
@test ql(A).factors ≈ ql!(Matrix(A)).factors
20+
@test ql(A).τ ≈ ql!(Matrix(A)).τ
21+
@test_broken Matrix(Q)*Matrix(L) ≈ A
22+
23+
for k=1:size(A,1),j=1:size(A,2)
24+
@test Q[k,j] ≈ Matrix(Q)[k,j]
25+
end
26+
27+
A=brand(T,10,14,3,2)
28+
Q,L=ql(A)
29+
@test ql(A).factors ≈ ql!(Matrix(A)).factors
30+
@test ql(A).τ ≈ ql!(Matrix(A)).τ
31+
@test Matrix(Q)*Matrix(L) ≈ A
32+
33+
for k=1:size(Q,1),j=1:size(Q,2)
34+
@test Q[k,j] ≈ Matrix(Q)[k,j]
35+
end
36+
37+
A=brand(T,10,14,3,6)
38+
Q,L=ql(A)
39+
@test ql(A).factors ≈ ql!(Matrix(A)).factors
40+
@test ql(A).τ ≈ ql!(Matrix(A)).τ
41+
@test Matrix(Q)*Matrix(L) ≈ A
42+
43+
for k=1:size(Q,1),j=1:size(Q,2)
44+
@test Q[k,j] ≈ Matrix(Q)[k,j]
45+
end
46+
47+
A=brand(T,100,100,3,4)
48+
@test ql(A).factors ≈ ql!(Matrix(A)).factors
49+
@test ql(A).τ ≈ ql!(Matrix(A)).τ
50+
b=rand(T,100)
51+
@test ql(A)\b ≈ Matrix(A)\b
52+
b=rand(T,100,2)
53+
@test ql(A)\b ≈ Matrix(A)\b
54+
@test_throws DimensionMismatch ql(A) \ randn(3)
55+
@test_throws DimensionMismatch ql(A).Q'randn(3)
56+
57+
A=brand(T,102,100,3,4)
58+
@test ql(A).factors ≈ ql!(Matrix(A)).factors
59+
@test ql(A).τ ≈ ql!(Matrix(A)).τ
60+
b=rand(T,102)
61+
@test_broken ql(A)\b ≈ Matrix(A)\b
62+
b=rand(T,102,2)
63+
@test_broken ql(A)\b ≈ Matrix(A)\b
64+
@test_throws DimensionMismatch ql(A) \ randn(3)
65+
@test_throws DimensionMismatch ql(A).Q'randn(3)
66+
67+
A=brand(T,100,102,3,4)
68+
@test ql(A).factors ≈ ql!(Matrix(A)).factors
69+
@test ql(A).τ ≈ ql!(Matrix(A)).τ
70+
b=rand(T,100)
71+
@test_broken ql(A)\b ≈ Matrix(A)\b
72+
73+
A = LinearAlgebra.Tridiagonal(randn(T,99), randn(T,100), randn(T,99))
74+
@test ql(A).factors ≈ ql!(Matrix(A)).factors
75+
@test ql(A).τ ≈ ql!(Matrix(A)).τ
76+
b=rand(T,100)
77+
@test ql(A)\b ≈ Matrix(A)\b
78+
b=rand(T,100,2)
79+
@test ql(A)\b ≈ Matrix(A)\b
80+
@test_throws DimensionMismatch ql(A) \ randn(3)
81+
@test_throws DimensionMismatch ql(A).Q'randn(3)
82+
end
83+
84+
@testset "lmul!/rmul!" begin
85+
for T in (Float32, Float64, ComplexF32, ComplexF64)
86+
A = brand(T,100,100,3,4)
87+
Q,R = qr(A)
88+
x = randn(T,100)
89+
b = randn(T,100,2)
90+
@test lmul!(Q, copy(x)) ≈ Matrix(Q)*x
91+
@test lmul!(Q, copy(b)) ≈ Matrix(Q)*b
92+
@test lmul!(Q', copy(x)) ≈ Matrix(Q)'*x
93+
@test lmul!(Q', copy(b)) ≈ Matrix(Q)'*b
94+
c = randn(T,2,100)
95+
@test rmul!(copy(c), Q) ≈ c*Matrix(Q)
96+
@test rmul!(copy(c), Q') ≈ c*Matrix(Q')
97+
98+
A = brand(T,100,100,3,4)
99+
Q,L = ql(A)
100+
x = randn(T,100)
101+
b = randn(T,100,2)
102+
@test lmul!(Q, copy(x)) ≈ Matrix(Q)*x
103+
@test lmul!(Q, copy(b)) ≈ Matrix(Q)*b
104+
@test lmul!(Q', copy(x)) ≈ Matrix(Q)'*x
105+
@test lmul!(Q', copy(b)) ≈ Matrix(Q)'*b
106+
c = randn(T,2,100)
107+
@test rmul!(copy(c), Q) ≈ c*Matrix(Q)
108+
@test rmul!(copy(c), Q') ≈ c*Matrix(Q')
109+
end
110+
end
111+
112+
@testset "Mixed types" begin
113+
A=brand(10,10,3,2)
114+
b=rand(ComplexF64,10)
115+
Q,L=ql(A)
116+
@test L\(Q'*b) ≈ ql(A)\b ≈ Matrix(A)\b
117+
@test Q*L ≈ A
118+
119+
A=brand(ComplexF64,10,10,3,2)
120+
b=rand(10)
121+
Q,L=ql(A)
122+
@test Q*L ≈ A
123+
@test L\(Q'*b) ≈ ql(A)\b ≈ Matrix(A)\b
124+
125+
A = BandedMatrix{Int}(undef, (2,1), (4,4))
126+
A.data .= 1:length(A.data)
127+
Q, L = ql(A)
128+
@test_broken Q*L ≈ A
129+
end
130+
131+
@testset "bounds" begin
132+
A = brand(100,100,3,4)
133+
@test_throws DimensionMismatch ql(A) \ randn(50)
134+
@test_throws DimensionMismatch ldiv!(ql(A), randn(50))
135+
Q = ql(A).Q
136+
@test_throws DimensionMismatch lmul!(Q, randn(50))
137+
@test_throws DimensionMismatch lmul!(Q', randn(50))
138+
@test_throws DimensionMismatch rmul!(randn(50)', Q)
139+
@test_throws DimensionMismatch rmul!(randn(50)', Q')
140+
end
141+
end
142+
end

0 commit comments

Comments
 (0)