Skip to content

Commit 34222ac

Browse files
dlfivefiftydlfivefifty
authored
BandedMatrices Extension (#64)
* BandedMatrices Extension * Update Project.toml * Add BandedMatrix extension and tests
1 parent cf5824f commit 34222ac

File tree

4 files changed

+322
-3
lines changed

4 files changed

+322
-3
lines changed

‎Project.toml

+11-2
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-dev"
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

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

‎test/runtests.jl

+3-1
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

+131
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,131 @@
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 Q*L ≈ A
129+
end
130+
end
131+
end

0 commit comments

Comments
 (0)