Skip to content

Commit c38b382

Browse files
mtfishmanmtfishman
authored
Implement left_orth/right_orth (#122)
1 parent 13d0332 commit c38b382

File tree

5 files changed

+165
-4
lines changed

5 files changed

+165
-4
lines changed

‎Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "BlockSparseArrays"
22
uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
33
authors = ["ITensor developers <[email protected]> and contributors"]
4-
version = "0.6.4"
4+
version = "0.6.5"
55

66
[deps]
77
Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

‎src/BlockSparseArrays.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,5 +49,6 @@ include("factorizations/truncation.jl")
4949
include("factorizations/qr.jl")
5050
include("factorizations/lq.jl")
5151
include("factorizations/polar.jl")
52+
include("factorizations/orthnull.jl")
5253

5354
end

‎src/factorizations/orthnull.jl

Lines changed: 124 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,124 @@
1+
using MatrixAlgebraKit:
2+
MatrixAlgebraKit,
3+
left_orth!,
4+
left_polar!,
5+
lq_compact!,
6+
qr_compact!,
7+
right_orth!,
8+
right_polar!,
9+
select_algorithm,
10+
svd_compact!
11+
12+
function MatrixAlgebraKit.initialize_output(
13+
::typeof(left_orth!), A::AbstractBlockSparseMatrix
14+
)
15+
return nothing
16+
end
17+
function MatrixAlgebraKit.check_input(::typeof(left_orth!), A::AbstractBlockSparseMatrix, F)
18+
!isnothing(F) && throw(
19+
ArgumentError(
20+
"`left_orth!` on block sparse matrices does not support specifying the output"
21+
),
22+
)
23+
return nothing
24+
end
25+
26+
function MatrixAlgebraKit.left_orth_qr!(A::AbstractBlockSparseMatrix, F, alg)
27+
!isnothing(F) && throw(
28+
ArgumentError(
29+
"`left_orth!` on block sparse matrices does not support specifying the output"
30+
),
31+
)
32+
alg′ = select_algorithm(qr_compact!, A, alg)
33+
return qr_compact!(A, alg′)
34+
end
35+
function MatrixAlgebraKit.left_orth_polar!(A::AbstractBlockSparseMatrix, F, alg)
36+
!isnothing(F) && throw(
37+
ArgumentError(
38+
"`left_orth!` on block sparse matrices does not support specifying the output"
39+
),
40+
)
41+
alg′ = select_algorithm(left_polar!, A, alg)
42+
return left_polar!(A, alg′)
43+
end
44+
function MatrixAlgebraKit.left_orth_svd!(
45+
A::AbstractBlockSparseMatrix, F, alg, trunc::Nothing=nothing
46+
)
47+
!isnothing(F) && throw(
48+
ArgumentError(
49+
"`left_orth!` on block sparse matrices does not support specifying the output"
50+
),
51+
)
52+
alg′ = select_algorithm(svd_compact!, A, alg)
53+
U, S, Vᴴ = svd_compact!(A, alg′)
54+
return U, S * Vᴴ
55+
end
56+
function MatrixAlgebraKit.left_orth_svd!(A::AbstractBlockSparseMatrix, F, alg, trunc)
57+
!isnothing(F) && throw(
58+
ArgumentError(
59+
"`left_orth!` on block sparse matrices does not support specifying the output"
60+
),
61+
)
62+
alg′ = select_algorithm(svd_compact!, A, alg)
63+
alg_trunc = select_algorithm(svd_trunc!, A, alg′; trunc)
64+
U, S, Vᴴ = svd_trunc!(A, alg_trunc)
65+
return U, S * Vᴴ
66+
end
67+
68+
function MatrixAlgebraKit.initialize_output(
69+
::typeof(right_orth!), A::AbstractBlockSparseMatrix
70+
)
71+
return nothing
72+
end
73+
function MatrixAlgebraKit.check_input(
74+
::typeof(right_orth!), A::AbstractBlockSparseMatrix, F::Nothing
75+
)
76+
!isnothing(F) && throw(
77+
ArgumentError(
78+
"`right_orth!` on block sparse matrices does not support specifying the output"
79+
),
80+
)
81+
return nothing
82+
end
83+
84+
function MatrixAlgebraKit.right_orth_lq!(A::AbstractBlockSparseMatrix, F, alg)
85+
!isnothing(F) && throw(
86+
ArgumentError(
87+
"`right_orth!` on block sparse matrices does not support specifying the output"
88+
),
89+
)
90+
alg′ = select_algorithm(lq_compact!, A, alg)
91+
return lq_compact!(A, alg′)
92+
end
93+
function MatrixAlgebraKit.right_orth_polar!(A::AbstractBlockSparseMatrix, F, alg)
94+
!isnothing(F) && throw(
95+
ArgumentError(
96+
"`right_orth!` on block sparse matrices does not support specifying the output"
97+
),
98+
)
99+
alg′ = select_algorithm(right_polar!, A, alg)
100+
return right_polar!(A, alg′)
101+
end
102+
function MatrixAlgebraKit.right_orth_svd!(
103+
A::AbstractBlockSparseMatrix, F, alg, trunc::Nothing=nothing
104+
)
105+
!isnothing(F) && throw(
106+
ArgumentError(
107+
"`right_orth!` on block sparse matrices does not support specifying the output"
108+
),
109+
)
110+
alg′ = select_algorithm(svd_compact!, A, alg)
111+
U, S, Vᴴ = svd_compact!(A, alg′)
112+
return U * S, Vᴴ
113+
end
114+
function MatrixAlgebraKit.right_orth_svd!(A::AbstractBlockSparseMatrix, F, alg, trunc)
115+
!isnothing(F) && throw(
116+
ArgumentError(
117+
"`right_orth!` on block sparse matrices does not support specifying the output"
118+
),
119+
)
120+
alg′ = select_algorithm(svd_compact!, A, alg)
121+
alg_trunc = select_algorithm(svd_trunc!, A, alg′; trunc)
122+
U, S, Vᴴ = svd_trunc!(A, alg_trunc)
123+
return U * S, Vᴴ
124+
end

‎src/factorizations/qr.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
using MatrixAlgebraKit: MatrixAlgebraKit, qr_compact!, qr_full!
1+
using MatrixAlgebraKit: MatrixAlgebraKit, lq_compact!, lq_full!, qr_compact!, qr_full!
22

33
# TODO: this is a hardcoded for now to get around this function not being defined in the
44
# type domain

‎test/test_factorizations.jl

Lines changed: 38 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,13 @@
11
using BlockArrays: Block, BlockedMatrix, BlockedVector, blocks, mortar
22
using BlockSparseArrays: BlockSparseArray, BlockDiagonal, eachblockstoredindex
33
using MatrixAlgebraKit:
4+
left_orth,
45
left_polar,
56
lq_compact,
67
lq_full,
78
qr_compact,
89
qr_full,
10+
right_orth,
911
right_polar,
1012
svd_compact,
1113
svd_full,
@@ -166,7 +168,7 @@ end
166168
end
167169
end
168170

169-
@testset "qr_compact" for T in (Float32, Float64, ComplexF32, ComplexF64)
171+
@testset "qr_compact (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
170172
for i in [2, 3], j in [2, 3], k in [2, 3], l in [2, 3]
171173
A = BlockSparseArray{T}(undef, ([i, j], [k, l]))
172174
A[Block(1, 1)] = randn(T, i, k)
@@ -177,7 +179,7 @@ end
177179
end
178180
end
179181

180-
@testset "qr_full" for T in (Float32, Float64, ComplexF32, ComplexF64)
182+
@testset "qr_full (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
181183
for i in [2, 3], j in [2, 3], k in [2, 3], l in [2, 3]
182184
A = BlockSparseArray{T}(undef, ([i, j], [k, l]))
183185
A[Block(1, 1)] = randn(T, i, k)
@@ -237,3 +239,37 @@ end
237239
@test C * U ≈ A
238240
@test Matrix(U * U') ≈ LinearAlgebra.I
239241
end
242+
243+
@testset "left_orth (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
244+
A = BlockSparseArray{T}(undef, ([3, 4], [2, 3]))
245+
A[Block(1, 1)] = randn(T, 3, 2)
246+
A[Block(2, 2)] = randn(T, 4, 3)
247+
248+
for kind in (:polar, :qr, :svd)
249+
U, C = left_orth(A; kind)
250+
@test U * C ≈ A
251+
@test Matrix(U'U) ≈ LinearAlgebra.I
252+
end
253+
254+
U, C = left_orth(A; trunc=(; maxrank=2))
255+
@test size(U, 2) ≤ 2
256+
@test size(C, 1) ≤ 2
257+
@test Matrix(U'U) ≈ LinearAlgebra.I
258+
end
259+
260+
@testset "right_orth (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
261+
A = BlockSparseArray{T}(undef, ([2, 3], [3, 4]))
262+
A[Block(1, 1)] = randn(T, 2, 3)
263+
A[Block(2, 2)] = randn(T, 3, 4)
264+
265+
for kind in (:lq, :polar, :svd)
266+
C, U = right_orth(A; kind)
267+
@test C * U ≈ A
268+
@test Matrix(U * U') ≈ LinearAlgebra.I
269+
end
270+
271+
C, U = right_orth(A; trunc=(; maxrank=2))
272+
@test size(C, 2) ≤ 2
273+
@test size(U, 1) ≤ 2
274+
@test Matrix(U * U') ≈ LinearAlgebra.I
275+
end

0 commit comments

Comments
 (0)