From 48c1ef6a7cd31badb06cdbd850fb16343ce68a32 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 6 Nov 2025 09:28:36 -0500 Subject: [PATCH 01/14] add testsuite --- test/testsuite/TestSuite.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 test/testsuite/TestSuite.jl diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl new file mode 100644 index 00000000..9adc4560 --- /dev/null +++ b/test/testsuite/TestSuite.jl @@ -0,0 +1,11 @@ +# Based on the design of GPUArrays.jl + +""" + TestSuite + +Suite of tests that may be used for all packages inheriting from MatrixAlgebraKit. + +""" +module TestSuite + +end From a118e5dcef8d0133825c235aaecd80cd4df22318 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 6 Nov 2025 10:46:13 -0500 Subject: [PATCH 02/14] add testsuite macro --- test/testsuite/TestSuite.jl | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 9adc4560..74a9cc34 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -8,4 +8,18 @@ Suite of tests that may be used for all packages inheriting from MatrixAlgebraKi """ module TestSuite +const tests = Dict() + +macro testsuite(name, ex) + safe_name = lowercase(replace(replace(name, " " => "_"), "/" => "_")) + fn = Symbol("test_", safe_name) + return quote + $(esc(fn))(AT; eltypes = supported_eltypes(AT, $(esc(fn)))) = $(esc(ex))(AT, eltypes) + @assert !haskey(tests, $name) "testsuite already exists" + tests[$name] = $fn + end +end + +end + end From 6dc2081e02d09ceeb660b79653a2bc4280517572 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 6 Nov 2025 10:47:14 -0500 Subject: [PATCH 03/14] add various test utility --- test/testsuite/TestSuite.jl | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 74a9cc34..34e5c927 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -8,6 +8,11 @@ Suite of tests that may be used for all packages inheriting from MatrixAlgebraKi """ module TestSuite +using Test, TestExtras +using MatrixAlgebraKit +using MatrixAlgebraKit: diagview +using LinearAlgebra: norm, istriu + const tests = Dict() macro testsuite(name, ex) @@ -20,6 +25,25 @@ macro testsuite(name, ex) end end +testargs_summary(args...) = string(args) + +instantiate_matrix(::Type{T}, size) where {T <: Number} = randn(T, size) +instantiate_matrix(::Type{AT}, size) where {AT <: Array} = randn(eltype(AT), size) + +function has_positive_diagonal(A) + T = eltype(A) + return if T <: Real + all(≥(zero(T)), diagview(A)) + else + all(≥(zero(real(T))), real(diagview(A))) && + all(≈(zero(real(T))), imag(diagview(A))) + end +end +isleftnull(N, A; kwargs...) = isapprox(norm(N' * A), 0; kwargs...) + +# TODO: actually make this a test +macro testinferred(ex) + return esc(:(@inferred $ex)) end end From 9b47a363029b6991a2a85ede85ade5babb36f35f Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 6 Nov 2025 10:47:21 -0500 Subject: [PATCH 04/14] add QR test suite --- test/testsuite/TestSuite.jl | 2 + test/testsuite/qr.jl | 154 ++++++++++++++++++++++++++++++++++++ 2 files changed, 156 insertions(+) create mode 100644 test/testsuite/qr.jl diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 34e5c927..b7e86afe 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -46,4 +46,6 @@ macro testinferred(ex) return esc(:(@inferred $ex)) end +include("qr.jl") + end diff --git a/test/testsuite/qr.jl b/test/testsuite/qr.jl new file mode 100644 index 00000000..b7f2e82e --- /dev/null +++ b/test/testsuite/qr.jl @@ -0,0 +1,154 @@ +function test_qr(::Type{T}, sz; kwargs...) where {T} + summary_str = testargs_summary(T, sz) + return @testset "qr $summary_str" begin + test_qr_compact(T, sz; kwargs...) + test_qr_full(T, sz; kwargs...) + test_qr_null(T, sz; kwargs...) + end +end + +function test_qr_compact( + ::Type{T}, sz; + test_positive = true, test_pivoted = true, test_blocksize = true, kwargs... + ) where {T <: Number} + summary_str = testargs_summary(T, sz) + return @testset "qr_compact! $summary_str" begin + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + + # does the elementary functionality work + Q, R = @testinferred qr_compact(A) + @test Q * R ≈ A + @test isisometric(Q) + @test istriu(R) + @test A == Ac + + # can I pass in outputs? + Q2, R2 = @testinferred qr_compact!(deepcopy(A), (Q, R)) + @test Q2 * R2 ≈ A + @test isisometric(Q2) + @test istriu(R2) + + # do we support `positive = true`? + if test_positive + Qpos, Rpos = @testinferred qr_compact(A; positive = true) + @test Qpos * Rpos ≈ A + @test isisometric(Qpos) + @test istriu(Rpos) + @test has_positive_diagonal(Rpos) + else + @test_throws ArgumentError qr_compact(A; positive = true) + end + + # do we support `pivoted = true`? + if test_pivoted + Qpiv, Rpiv = @testinferred qr_compact(A; pivoted = true) + @test Qpiv * Rpiv ≈ A + @test isisometric(Qpos) + else + @test_throws ArgumentError qr_compact(A; pivoted = true) + end + + # do we support `blocksize = Int`? + if test_blocksize + Qblocked, Rblocked = @testinferred qr_compact(A; blocksize = 2) + @test Qblocked * Rblocked ≈ A + @test isisometric(Qblocked) + else + @test_throws ArgumentError qr_compact(A; blocksize = 2) + end + end +end + +function test_qr_full( + ::Type{T}, sz; + test_positive = true, test_pivoted = true, test_blocksize = true, kwargs... + ) where {T <: Number} + summary_str = testargs_summary(T, sz) + return @testset "qr_full! $summary_str" begin + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + + # does the elementary functionality work + Q, R = @testinferred qr_full(A) + @test Q * R ≈ A + @test isunitary(Q) + @test istriu(R) + @test A == Ac + + # can I pass in outputs? + Q2, R2 = @testinferred qr_full!(deepcopy(A), (Q, R)) + @test Q2 * R2 ≈ A + @test isunitary(Q2) + @test istriu(R2) + + # do we support `positive = true`? + if test_positive + Qpos, Rpos = @testinferred qr_full(A; positive = true) + @test Qpos * Rpos ≈ A + @test isunitary(Qpos) + @test istriu(Rpos) + @test has_positive_diagonal(Rpos) + else + @test_throws ArgumentError qr_full(A; positive = true) + end + + # do we support `pivoted = true`? + if test_pivoted + Qpiv, Rpiv = @testinferred qr_full(A; pivoted = true) + @test Qpiv * Rpiv ≈ A + @test isunitary(Qpos) + else + @test_throws ArgumentError qr_full(A; pivoted = true) + end + + # do we support `blocksize = Int`? + if test_blocksize + Qblocked, Rblocked = @testinferred qr_full(A; blocksize = 2) + @test Qblocked * Rblocked ≈ A + @test isunitary(Qblocked) + else + @test_throws ArgumentError qr_full(A; blocksize = 2) + end + end +end + +function test_qr_null( + ::Type{T}, sz; + test_pivoted = true, test_blocksize = true, kwargs... + ) where {T <: Number} + summary_str = testargs_summary(T, sz) + return @testset "qr_null! $summary_str" begin + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + + # does the elementary functionality work + N = @testinferred qr_null(A) + @test isleftnull(N, A) + @test isisometric(N) + @test A == Ac + + # can I pass in outputs? + N2 = @testinferred qr_null!(deepcopy(A), N) + @test isleftnull(N2, A) + @test isisometric(N2) + + # do we support `pivoted = true`? + if test_pivoted + Npiv = @testinferred qr_null(A; pivoted = true) + @test isleftnull(Npiv, A) + @test isisometric(Npiv) + else + @test_throws ArgumentError qr_null(A; pivoted = true) + end + + # do we support `blocksize = Int`? + if test_blocksize + Nblocked = @testinferred qr_null(A; blocksize = 2) + @test isleftnull(Nblocked, A) + @test isisometric(Nblocked) + else + @test_throws ArgumentError qr_null(A; blocksize = 2) + end + end +end From 280f38ddd304f73c0d3219caeef3a27b981f4d41 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 6 Nov 2025 11:52:21 -0500 Subject: [PATCH 05/14] add tolerances --- test/testsuite/TestSuite.jl | 6 +++- test/testsuite/qr.jl | 62 ++++++++++++++++++++----------------- 2 files changed, 39 insertions(+), 29 deletions(-) diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index b7e86afe..c5233fb8 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -30,6 +30,9 @@ testargs_summary(args...) = string(args) instantiate_matrix(::Type{T}, size) where {T <: Number} = randn(T, size) instantiate_matrix(::Type{AT}, size) where {AT <: Array} = randn(eltype(AT), size) +precision(::Type{T}) where {T <: Number} = sqrt(eps(real(T))) +precision(::Type{T}) where {T} = precision(eltype(T)) + function has_positive_diagonal(A) T = eltype(A) return if T <: Real @@ -39,7 +42,8 @@ function has_positive_diagonal(A) all(≈(zero(real(T))), imag(diagview(A))) end end -isleftnull(N, A; kwargs...) = isapprox(norm(N' * A), 0; kwargs...) +isleftnull(N, A; atol::Real = 0, rtol::Real = precision(eltype(A))) = + isapprox(norm(A' * N), 0; atol = max(atol, norm(A) * rtol)) # TODO: actually make this a test macro testinferred(ex) diff --git a/test/testsuite/qr.jl b/test/testsuite/qr.jl index b7f2e82e..5f2c0a72 100644 --- a/test/testsuite/qr.jl +++ b/test/testsuite/qr.jl @@ -1,4 +1,4 @@ -function test_qr(::Type{T}, sz; kwargs...) where {T} +function test_qr(T::Type, sz; kwargs...) summary_str = testargs_summary(T, sz) return @testset "qr $summary_str" begin test_qr_compact(T, sz; kwargs...) @@ -8,9 +8,11 @@ function test_qr(::Type{T}, sz; kwargs...) where {T} end function test_qr_compact( - ::Type{T}, sz; - test_positive = true, test_pivoted = true, test_blocksize = true, kwargs... - ) where {T <: Number} + T::Type, sz; + test_positive = true, test_pivoted = true, test_blocksize = true, + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) summary_str = testargs_summary(T, sz) return @testset "qr_compact! $summary_str" begin A = instantiate_matrix(T, sz) @@ -19,21 +21,21 @@ function test_qr_compact( # does the elementary functionality work Q, R = @testinferred qr_compact(A) @test Q * R ≈ A - @test isisometric(Q) + @test isisometric(Q; atol, rtol) @test istriu(R) @test A == Ac # can I pass in outputs? Q2, R2 = @testinferred qr_compact!(deepcopy(A), (Q, R)) @test Q2 * R2 ≈ A - @test isisometric(Q2) + @test isisometric(Q2; atol, rtol) @test istriu(R2) # do we support `positive = true`? if test_positive Qpos, Rpos = @testinferred qr_compact(A; positive = true) @test Qpos * Rpos ≈ A - @test isisometric(Qpos) + @test isisometric(Qpos; atol, rtol) @test istriu(Rpos) @test has_positive_diagonal(Rpos) else @@ -44,7 +46,7 @@ function test_qr_compact( if test_pivoted Qpiv, Rpiv = @testinferred qr_compact(A; pivoted = true) @test Qpiv * Rpiv ≈ A - @test isisometric(Qpos) + @test isisometric(Qpos; atol, rtol) else @test_throws ArgumentError qr_compact(A; pivoted = true) end @@ -53,7 +55,7 @@ function test_qr_compact( if test_blocksize Qblocked, Rblocked = @testinferred qr_compact(A; blocksize = 2) @test Qblocked * Rblocked ≈ A - @test isisometric(Qblocked) + @test isisometric(Qblocked; atol, rtol) else @test_throws ArgumentError qr_compact(A; blocksize = 2) end @@ -61,9 +63,11 @@ function test_qr_compact( end function test_qr_full( - ::Type{T}, sz; - test_positive = true, test_pivoted = true, test_blocksize = true, kwargs... - ) where {T <: Number} + T::Type, sz; + test_positive = true, test_pivoted = true, test_blocksize = true, + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) summary_str = testargs_summary(T, sz) return @testset "qr_full! $summary_str" begin A = instantiate_matrix(T, sz) @@ -72,21 +76,21 @@ function test_qr_full( # does the elementary functionality work Q, R = @testinferred qr_full(A) @test Q * R ≈ A - @test isunitary(Q) + @test isunitary(Q; atol, rtol) @test istriu(R) @test A == Ac # can I pass in outputs? Q2, R2 = @testinferred qr_full!(deepcopy(A), (Q, R)) @test Q2 * R2 ≈ A - @test isunitary(Q2) + @test isunitary(Q2; atol, rtol) @test istriu(R2) # do we support `positive = true`? if test_positive Qpos, Rpos = @testinferred qr_full(A; positive = true) @test Qpos * Rpos ≈ A - @test isunitary(Qpos) + @test isunitary(Qpos; atol, rtol) @test istriu(Rpos) @test has_positive_diagonal(Rpos) else @@ -97,7 +101,7 @@ function test_qr_full( if test_pivoted Qpiv, Rpiv = @testinferred qr_full(A; pivoted = true) @test Qpiv * Rpiv ≈ A - @test isunitary(Qpos) + @test isunitary(Qpos; atol, rtol) else @test_throws ArgumentError qr_full(A; pivoted = true) end @@ -106,7 +110,7 @@ function test_qr_full( if test_blocksize Qblocked, Rblocked = @testinferred qr_full(A; blocksize = 2) @test Qblocked * Rblocked ≈ A - @test isunitary(Qblocked) + @test isunitary(Qblocked; atol, rtol) else @test_throws ArgumentError qr_full(A; blocksize = 2) end @@ -114,9 +118,11 @@ function test_qr_full( end function test_qr_null( - ::Type{T}, sz; - test_pivoted = true, test_blocksize = true, kwargs... - ) where {T <: Number} + T::Type, sz; + test_pivoted = true, test_blocksize = true, + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) summary_str = testargs_summary(T, sz) return @testset "qr_null! $summary_str" begin A = instantiate_matrix(T, sz) @@ -124,20 +130,20 @@ function test_qr_null( # does the elementary functionality work N = @testinferred qr_null(A) - @test isleftnull(N, A) - @test isisometric(N) + @test isleftnull(N, A; atol, rtol) + @test isisometric(N; atol, rtol) @test A == Ac # can I pass in outputs? N2 = @testinferred qr_null!(deepcopy(A), N) - @test isleftnull(N2, A) - @test isisometric(N2) + @test isleftnull(N2, A; atol, rtol) + @test isisometric(N2; atol, rtol) # do we support `pivoted = true`? if test_pivoted Npiv = @testinferred qr_null(A; pivoted = true) - @test isleftnull(Npiv, A) - @test isisometric(Npiv) + @test isleftnull(Npiv, A; atol, rtol) + @test isisometric(Npiv; atol, rtol) else @test_throws ArgumentError qr_null(A; pivoted = true) end @@ -145,8 +151,8 @@ function test_qr_null( # do we support `blocksize = Int`? if test_blocksize Nblocked = @testinferred qr_null(A; blocksize = 2) - @test isleftnull(Nblocked, A) - @test isisometric(Nblocked) + @test isleftnull(Nblocked, A; atol, rtol) + @test isisometric(Nblocked; atol, rtol) else @test_throws ArgumentError qr_null(A; blocksize = 2) end From caaeca95400926e0f26bdff0fdb7061e5e99caeb Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 6 Nov 2025 11:52:48 -0500 Subject: [PATCH 06/14] add StableRNG --- test/testsuite/TestSuite.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index c5233fb8..bd9496a6 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -12,6 +12,7 @@ using Test, TestExtras using MatrixAlgebraKit using MatrixAlgebraKit: diagview using LinearAlgebra: norm, istriu +using Random, StableRNGs const tests = Dict() @@ -27,8 +28,11 @@ end testargs_summary(args...) = string(args) -instantiate_matrix(::Type{T}, size) where {T <: Number} = randn(T, size) -instantiate_matrix(::Type{AT}, size) where {AT <: Array} = randn(eltype(AT), size) +const rng = StableRNG(123) +seed_rng!(seed) = Random.seed!(rng, seed) + +instantiate_matrix(::Type{T}, size) where {T <: Number} = randn(rng, T, size) +instantiate_matrix(::Type{AT}, size) where {AT <: Array} = randn(rng, eltype(AT), size) precision(::Type{T}) where {T <: Number} = sqrt(eps(real(T))) precision(::Type{T}) where {T} = precision(eltype(T)) From 3d36a3256d7ae10ffc360f67ca2a981e99aaedd9 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 6 Nov 2025 11:53:04 -0500 Subject: [PATCH 07/14] fix error strictness --- test/testsuite/qr.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/test/testsuite/qr.jl b/test/testsuite/qr.jl index 5f2c0a72..c354c2e4 100644 --- a/test/testsuite/qr.jl +++ b/test/testsuite/qr.jl @@ -39,7 +39,7 @@ function test_qr_compact( @test istriu(Rpos) @test has_positive_diagonal(Rpos) else - @test_throws ArgumentError qr_compact(A; positive = true) + @test_throws Exception qr_compact(A; positive = true) end # do we support `pivoted = true`? @@ -48,7 +48,7 @@ function test_qr_compact( @test Qpiv * Rpiv ≈ A @test isisometric(Qpos; atol, rtol) else - @test_throws ArgumentError qr_compact(A; pivoted = true) + @test_throws Exception qr_compact(A; pivoted = true) end # do we support `blocksize = Int`? @@ -57,7 +57,7 @@ function test_qr_compact( @test Qblocked * Rblocked ≈ A @test isisometric(Qblocked; atol, rtol) else - @test_throws ArgumentError qr_compact(A; blocksize = 2) + @test_throws Exception qr_compact(A; blocksize = 2) end end end @@ -94,7 +94,7 @@ function test_qr_full( @test istriu(Rpos) @test has_positive_diagonal(Rpos) else - @test_throws ArgumentError qr_full(A; positive = true) + @test_throws Exception qr_full(A; positive = true) end # do we support `pivoted = true`? @@ -103,7 +103,7 @@ function test_qr_full( @test Qpiv * Rpiv ≈ A @test isunitary(Qpos; atol, rtol) else - @test_throws ArgumentError qr_full(A; pivoted = true) + @test_throws Exception qr_full(A; pivoted = true) end # do we support `blocksize = Int`? @@ -112,7 +112,7 @@ function test_qr_full( @test Qblocked * Rblocked ≈ A @test isunitary(Qblocked; atol, rtol) else - @test_throws ArgumentError qr_full(A; blocksize = 2) + @test_throws Exception qr_full(A; blocksize = 2) end end end @@ -145,7 +145,7 @@ function test_qr_null( @test isleftnull(Npiv, A; atol, rtol) @test isisometric(Npiv; atol, rtol) else - @test_throws ArgumentError qr_null(A; pivoted = true) + @test_throws Exception qr_null(A; pivoted = true) end # do we support `blocksize = Int`? @@ -154,7 +154,7 @@ function test_qr_null( @test isleftnull(Nblocked, A; atol, rtol) @test isisometric(Nblocked; atol, rtol) else - @test_throws ArgumentError qr_null(A; blocksize = 2) + @test_throws Exception qr_null(A; blocksize = 2) end end end From b166b1c4b9c782f5ce6854a238371e1146c8ac53 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 6 Nov 2025 11:53:11 -0500 Subject: [PATCH 08/14] add Diagonal support --- test/testsuite/TestSuite.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index bd9496a6..c75abe1f 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -11,7 +11,7 @@ module TestSuite using Test, TestExtras using MatrixAlgebraKit using MatrixAlgebraKit: diagview -using LinearAlgebra: norm, istriu +using LinearAlgebra: Diagonal, norm, istriu using Random, StableRNGs const tests = Dict() @@ -33,6 +33,7 @@ seed_rng!(seed) = Random.seed!(rng, seed) instantiate_matrix(::Type{T}, size) where {T <: Number} = randn(rng, T, size) instantiate_matrix(::Type{AT}, size) where {AT <: Array} = randn(rng, eltype(AT), size) +instantiate_matrix(::Type{AT}, size) where {AT <: Diagonal} = Diagonal(randn(rng, eltype(AT), size)) precision(::Type{T}) where {T <: Number} = sqrt(eps(real(T))) precision(::Type{T}) where {T} = precision(eltype(T)) From 56db38b19961da3ca081392a2298b7f40496faf6 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 6 Nov 2025 11:53:21 -0500 Subject: [PATCH 09/14] replace testsuite --- test/qr.jl | 432 +++++++++++++++++++++++++++-------------------------- 1 file changed, 223 insertions(+), 209 deletions(-) diff --git a/test/qr.jl b/test/qr.jl index c4f0c9d6..ed8c255f 100644 --- a/test/qr.jl +++ b/test/qr.jl @@ -7,217 +7,231 @@ using LinearAlgebra: diag, I, Diagonal BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) GenericFloats = (Float16, BigFloat, Complex{BigFloat}) -@testset "qr_compact! and qr_null! for T = $T" for T in BLASFloats - rng = StableRNG(123) - m = 54 - for n in (37, m, 63) - minmn = min(m, n) - A = randn(rng, T, m, n) - Q, R = @constinferred qr_compact(A) - @test Q isa Matrix{T} && size(Q) == (m, minmn) - @test R isa Matrix{T} && size(R) == (minmn, n) - @test Q * R ≈ A - N = @constinferred qr_null(A) - @test N isa Matrix{T} && size(N) == (m, m - minmn) - @test isisometric(Q) - @test maximum(abs, A' * N) < eps(real(T))^(2 / 3) - @test isisometric(N) +@isdefined(TestSuite) || include("testsuite/TestSuite.jl") +using .TestSuite - Ac = similar(A) - Q2, R2 = @constinferred qr_compact!(copy!(Ac, A), (Q, R)) - @test Q2 === Q - @test R2 === R - N2 = @constinferred qr_null!(copy!(Ac, A), N) - @test N2 === N - - Q2 = similar(Q) - noR = similar(A, minmn, 0) - qr_compact!(copy!(Ac, A), (Q2, noR)) - @test Q == Q2 - - # unblocked algorithm - qr_compact!(copy!(Ac, A), (Q, R); blocksize = 1) - @test Q * R ≈ A - @test isisometric(Q) - qr_compact!(copy!(Ac, A), (Q2, noR); blocksize = 1) - @test Q == Q2 - qr_compact!(copy!(Ac, A), (Q2, noR); blocksize = 1) - qr_null!(copy!(Ac, A), N; blocksize = 1) - @test maximum(abs, A' * N) < eps(real(T))^(2 / 3) - @test isisometric(N) - if n <= m - qr_compact!(copy!(Q2, A), (Q2, noR); blocksize = 1) # in-place Q - @test Q ≈ Q2 - @test_throws ArgumentError qr_compact!(copy!(Q2, A), (Q2, R); blocksize = 1) - @test_throws ArgumentError qr_compact!(copy!(Q2, A), (Q2, noR); positive = true) - @test_throws ArgumentError qr_compact!(copy!(Q2, A), (Q2, noR); blocksize = 8) - end - # other blocking - qr_compact!(copy!(Ac, A), (Q, R); blocksize = 8) - @test Q * R ≈ A - @test isisometric(Q) - qr_compact!(copy!(Ac, A), (Q2, noR); blocksize = 8) - @test Q == Q2 - qr_null!(copy!(Ac, A), N; blocksize = 8) - @test maximum(abs, A' * N) < eps(real(T))^(2 / 3) - @test isisometric(N) - - # pivoted - qr_compact!(copy!(Ac, A), (Q, R); pivoted = true) - @test Q * R ≈ A - @test Q' * Q ≈ I - qr_compact!(copy!(Ac, A), (Q2, noR); pivoted = true) - @test Q == Q2 - # positive - qr_compact!(copy!(Ac, A), (Q, R); positive = true) - @test Q * R ≈ A - @test isisometric(Q) - @test all(>=(zero(real(T))), real(diag(R))) - qr_compact!(copy!(Ac, A), (Q2, noR); positive = true) - @test Q == Q2 - # positive and blocksize 1 - qr_compact!(copy!(Ac, A), (Q, R); positive = true, blocksize = 1) - @test Q * R ≈ A - @test isisometric(Q) - @test all(>=(zero(real(T))), real(diag(R))) - qr_compact!(copy!(Ac, A), (Q2, noR); positive = true, blocksize = 1) - @test Q == Q2 - # positive and pivoted - qr_compact!(copy!(Ac, A), (Q, R); positive = true, pivoted = true) - @test Q * R ≈ A - @test isisometric(Q) - if n <= m - # the following test tries to find the diagonal element (in order to test positivity) - # before the column permutation. This only works if all columns have a diagonal - # element - for j in 1:n - i = findlast(!iszero, view(R, :, j)) - @test real(R[i, j]) >= zero(real(T)) - end - end - qr_compact!(copy!(Ac, A), (Q2, noR); positive = true, pivoted = true) - @test Q == Q2 - end +m = 54 +for T in BLASFloats, n in (37, m, 63) + TestSuite.seed_rng!(123) + TestSuite.test_qr(T, (m, n)) end - -@testset "qr_full! for T = $T" for T in BLASFloats - rng = StableRNG(123) - m = 54 - for n in (37, m, 63) - minmn = min(m, n) - A = randn(rng, T, m, n) - Q, R = qr_full(A) - @test Q isa Matrix{T} && size(Q) == (m, m) - @test R isa Matrix{T} && size(R) == (m, n) - @test Q * R ≈ A - @test isunitary(Q) - - Ac = similar(A) - Q2 = similar(Q) - noR = similar(A, m, 0) - Q2, R2 = @constinferred qr_full!(copy!(Ac, A), (Q, R)) - @test Q2 === Q - @test R2 === R - @test Q * R ≈ A - @test isunitary(Q) - qr_full!(copy!(Ac, A), (Q2, noR)) - @test Q == Q2 - - # unblocked algorithm - qr_full!(copy!(Ac, A), (Q, R); blocksize = 1) - @test Q * R ≈ A - @test isunitary(Q) - qr_full!(copy!(Ac, A), (Q2, noR); blocksize = 1) - @test Q == Q2 - if n == m - qr_full!(copy!(Q2, A), (Q2, noR); blocksize = 1) # in-place Q - @test Q ≈ Q2 - end - # other blocking - qr_full!(copy!(Ac, A), (Q, R); blocksize = 8) - @test Q * R ≈ A - @test isunitary(Q) - qr_full!(copy!(Ac, A), (Q2, noR); blocksize = 8) - @test Q == Q2 - # pivoted - qr_full!(copy!(Ac, A), (Q, R); pivoted = true) - @test Q * R ≈ A - @test isunitary(Q) - qr_full!(copy!(Ac, A), (Q2, noR); pivoted = true) - @test Q == Q2 - # positive - qr_full!(copy!(Ac, A), (Q, R); positive = true) - @test Q * R ≈ A - @test isunitary(Q) - @test all(>=(zero(real(T))), real(diag(R))) - qr_full!(copy!(Ac, A), (Q2, noR); positive = true) - @test Q == Q2 - # positive and blocksize 1 - qr_full!(copy!(Ac, A), (Q, R); positive = true, blocksize = 1) - @test Q * R ≈ A - @test isunitary(Q) - @test all(>=(zero(real(T))), real(diag(R))) - qr_full!(copy!(Ac, A), (Q2, noR); positive = true, blocksize = 1) - @test Q == Q2 - # positive and pivoted - qr_full!(copy!(Ac, A), (Q, R); positive = true, pivoted = true) - @test Q * R ≈ A - @test isunitary(Q) - if n <= m - # the following test tries to find the diagonal element (in order to test positivity) - # before the column permutation. This only works if all columns have a diagonal - # element - for j in 1:n - i = findlast(!iszero, view(R, :, j)) - @test real(R[i, j]) >= zero(real(T)) - end - end - qr_full!(copy!(Ac, A), (Q2, noR); positive = true, pivoted = true) - @test Q == Q2 - end +for T in (BLASFloats..., GenericFloats...) + AT = Diagonal{T, Vector{T}} + TestSuite.test_qr(AT, m; test_pivoted = false, test_blocksize = false) end -@testset "qr_compact, qr_full and qr_null for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) - rng = StableRNG(123) - atol = eps(real(T))^(3 / 4) - for m in (54, 0) - Ad = randn(rng, T, m) - A = Diagonal(Ad) - # compact - Q, R = @constinferred qr_compact(A) - @test Q isa Diagonal{T} && size(Q) == (m, m) - @test R isa Diagonal{T} && size(R) == (m, m) - @test Q * R ≈ A - @test isunitary(Q) - - # compact and positive - Qp, Rp = @constinferred qr_compact(A; positive = true) - @test Qp isa Diagonal{T} && size(Qp) == (m, m) - @test Rp isa Diagonal{T} && size(Rp) == (m, m) - @test Qp * Rp ≈ A - @test isunitary(Qp) - @test all(≥(zero(real(T))), real(diag(Rp))) && - all(≈(zero(real(T)); atol), imag(diag(Rp))) - - # full - Q, R = @constinferred qr_full(A) - @test Q isa Diagonal{T} && size(Q) == (m, m) - @test R isa Diagonal{T} && size(R) == (m, m) - @test Q * R ≈ A - @test isunitary(Q) - - # full and positive - Qp, Rp = @constinferred qr_full(A; positive = true) - @test Qp isa Diagonal{T} && size(Qp) == (m, m) - @test Rp isa Diagonal{T} && size(Rp) == (m, m) - @test Qp * Rp ≈ A - @test isunitary(Qp) - @test all(≥(zero(real(T))), real(diag(Rp))) && - all(≈(zero(real(T)); atol), imag(diag(Rp))) - - # null - N = @constinferred qr_null(A) - @test N isa AbstractMatrix{T} && size(N) == (m, 0) - end -end +# @testset "qr_compact! and qr_null! for T = $T" for T in BLASFloats +# rng = StableRNG(123) +# m = 54 +# for n in (37, m, 63) +# minmn = min(m, n) +# A = randn(rng, T, m, n) +# Q, R = @constinferred qr_compact(A) +# @test Q isa Matrix{T} && size(Q) == (m, minmn) +# @test R isa Matrix{T} && size(R) == (minmn, n) +# @test Q * R ≈ A +# N = @constinferred qr_null(A) +# @test N isa Matrix{T} && size(N) == (m, m - minmn) +# @test isisometric(Q) +# @test maximum(abs, A' * N) < eps(real(T))^(2 / 3) +# @test isisometric(N) +# +# Ac = similar(A) +# Q2, R2 = @constinferred qr_compact!(copy!(Ac, A), (Q, R)) +# @test Q2 === Q +# @test R2 === R +# N2 = @constinferred qr_null!(copy!(Ac, A), N) +# @test N2 === N +# +# Q2 = similar(Q) +# noR = similar(A, minmn, 0) +# qr_compact!(copy!(Ac, A), (Q2, noR)) +# @test Q == Q2 +# +# # unblocked algorithm +# qr_compact!(copy!(Ac, A), (Q, R); blocksize = 1) +# @test Q * R ≈ A +# @test isisometric(Q) +# qr_compact!(copy!(Ac, A), (Q2, noR); blocksize = 1) +# @test Q == Q2 +# qr_compact!(copy!(Ac, A), (Q2, noR); blocksize = 1) +# qr_null!(copy!(Ac, A), N; blocksize = 1) +# @test maximum(abs, A' * N) < eps(real(T))^(2 / 3) +# @test isisometric(N) +# if n <= m +# qr_compact!(copy!(Q2, A), (Q2, noR); blocksize = 1) # in-place Q +# @test Q ≈ Q2 +# @test_throws ArgumentError qr_compact!(copy!(Q2, A), (Q2, R); blocksize = 1) +# @test_throws ArgumentError qr_compact!(copy!(Q2, A), (Q2, noR); positive = true) +# @test_throws ArgumentError qr_compact!(copy!(Q2, A), (Q2, noR); blocksize = 8) +# end +# # other blocking +# qr_compact!(copy!(Ac, A), (Q, R); blocksize = 8) +# @test Q * R ≈ A +# @test isisometric(Q) +# qr_compact!(copy!(Ac, A), (Q2, noR); blocksize = 8) +# @test Q == Q2 +# qr_null!(copy!(Ac, A), N; blocksize = 8) +# @test maximum(abs, A' * N) < eps(real(T))^(2 / 3) +# @test isisometric(N) +# +# # pivoted +# qr_compact!(copy!(Ac, A), (Q, R); pivoted = true) +# @test Q * R ≈ A +# @test Q' * Q ≈ I +# qr_compact!(copy!(Ac, A), (Q2, noR); pivoted = true) +# @test Q == Q2 +# # positive +# qr_compact!(copy!(Ac, A), (Q, R); positive = true) +# @test Q * R ≈ A +# @test isisometric(Q) +# @test all(>=(zero(real(T))), real(diag(R))) +# qr_compact!(copy!(Ac, A), (Q2, noR); positive = true) +# @test Q == Q2 +# # positive and blocksize 1 +# qr_compact!(copy!(Ac, A), (Q, R); positive = true, blocksize = 1) +# @test Q * R ≈ A +# @test isisometric(Q) +# @test all(>=(zero(real(T))), real(diag(R))) +# qr_compact!(copy!(Ac, A), (Q2, noR); positive = true, blocksize = 1) +# @test Q == Q2 +# # positive and pivoted +# qr_compact!(copy!(Ac, A), (Q, R); positive = true, pivoted = true) +# @test Q * R ≈ A +# @test isisometric(Q) +# if n <= m +# # the following test tries to find the diagonal element (in order to test positivity) +# # before the column permutation. This only works if all columns have a diagonal +# # element +# for j in 1:n +# i = findlast(!iszero, view(R, :, j)) +# @test real(R[i, j]) >= zero(real(T)) +# end +# end +# qr_compact!(copy!(Ac, A), (Q2, noR); positive = true, pivoted = true) +# @test Q == Q2 +# end +# end +# +# @testset "qr_full! for T = $T" for T in BLASFloats +# rng = StableRNG(123) +# m = 54 +# for n in (37, m, 63) +# minmn = min(m, n) +# A = randn(rng, T, m, n) +# Q, R = qr_full(A) +# @test Q isa Matrix{T} && size(Q) == (m, m) +# @test R isa Matrix{T} && size(R) == (m, n) +# @test Q * R ≈ A +# @test isunitary(Q) +# +# Ac = similar(A) +# Q2 = similar(Q) +# noR = similar(A, m, 0) +# Q2, R2 = @constinferred qr_full!(copy!(Ac, A), (Q, R)) +# @test Q2 === Q +# @test R2 === R +# @test Q * R ≈ A +# @test isunitary(Q) +# qr_full!(copy!(Ac, A), (Q2, noR)) +# @test Q == Q2 +# +# # unblocked algorithm +# qr_full!(copy!(Ac, A), (Q, R); blocksize = 1) +# @test Q * R ≈ A +# @test isunitary(Q) +# qr_full!(copy!(Ac, A), (Q2, noR); blocksize = 1) +# @test Q == Q2 +# if n == m +# qr_full!(copy!(Q2, A), (Q2, noR); blocksize = 1) # in-place Q +# @test Q ≈ Q2 +# end +# # other blocking +# qr_full!(copy!(Ac, A), (Q, R); blocksize = 8) +# @test Q * R ≈ A +# @test isunitary(Q) +# qr_full!(copy!(Ac, A), (Q2, noR); blocksize = 8) +# @test Q == Q2 +# # pivoted +# qr_full!(copy!(Ac, A), (Q, R); pivoted = true) +# @test Q * R ≈ A +# @test isunitary(Q) +# qr_full!(copy!(Ac, A), (Q2, noR); pivoted = true) +# @test Q == Q2 +# # positive +# qr_full!(copy!(Ac, A), (Q, R); positive = true) +# @test Q * R ≈ A +# @test isunitary(Q) +# @test all(>=(zero(real(T))), real(diag(R))) +# qr_full!(copy!(Ac, A), (Q2, noR); positive = true) +# @test Q == Q2 +# # positive and blocksize 1 +# qr_full!(copy!(Ac, A), (Q, R); positive = true, blocksize = 1) +# @test Q * R ≈ A +# @test isunitary(Q) +# @test all(>=(zero(real(T))), real(diag(R))) +# qr_full!(copy!(Ac, A), (Q2, noR); positive = true, blocksize = 1) +# @test Q == Q2 +# # positive and pivoted +# qr_full!(copy!(Ac, A), (Q, R); positive = true, pivoted = true) +# @test Q * R ≈ A +# @test isunitary(Q) +# if n <= m +# # the following test tries to find the diagonal element (in order to test positivity) +# # before the column permutation. This only works if all columns have a diagonal +# # element +# for j in 1:n +# i = findlast(!iszero, view(R, :, j)) +# @test real(R[i, j]) >= zero(real(T)) +# end +# end +# qr_full!(copy!(Ac, A), (Q2, noR); positive = true, pivoted = true) +# @test Q == Q2 +# end +# end +# +# @testset "qr_compact, qr_full and qr_null for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) +# rng = StableRNG(123) +# atol = eps(real(T))^(3 / 4) +# for m in (54, 0) +# Ad = randn(rng, T, m) +# A = Diagonal(Ad) +# +# # compact +# Q, R = @constinferred qr_compact(A) +# @test Q isa Diagonal{T} && size(Q) == (m, m) +# @test R isa Diagonal{T} && size(R) == (m, m) +# @test Q * R ≈ A +# @test isunitary(Q) +# +# # compact and positive +# Qp, Rp = @constinferred qr_compact(A; positive = true) +# @test Qp isa Diagonal{T} && size(Qp) == (m, m) +# @test Rp isa Diagonal{T} && size(Rp) == (m, m) +# @test Qp * Rp ≈ A +# @test isunitary(Qp) +# @test all(≥(zero(real(T))), real(diag(Rp))) && +# all(≈(zero(real(T)); atol), imag(diag(Rp))) +# +# # full +# Q, R = @constinferred qr_full(A) +# @test Q isa Diagonal{T} && size(Q) == (m, m) +# @test R isa Diagonal{T} && size(R) == (m, m) +# @test Q * R ≈ A +# @test isunitary(Q) +# +# # full and positive +# Qp, Rp = @constinferred qr_full(A; positive = true) +# @test Qp isa Diagonal{T} && size(Qp) == (m, m) +# @test Rp isa Diagonal{T} && size(Rp) == (m, m) +# @test Qp * Rp ≈ A +# @test isunitary(Qp) +# @test all(≥(zero(real(T))), real(diag(Rp))) && +# all(≈(zero(real(T)); atol), imag(diag(Rp))) +# +# # null +# N = @constinferred qr_null(A) +# @test N isa AbstractMatrix{T} && size(N) == (m, 0) +# end +# end From 61b2efa2abecca825fcf8f413f7f6c011ce7aea3 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Wed, 12 Nov 2025 08:55:43 -0500 Subject: [PATCH 10/14] Test GPU QR in testsuite --- Project.toml | 4 +- src/implementations/qr.jl | 8 ++- test/amd/qr.jl | 123 ------------------------------------ test/cuda/qr.jl | 123 ------------------------------------ test/qr.jl | 9 +++ test/runtests.jl | 6 -- test/testsuite/TestSuite.jl | 5 ++ test/testsuite/qr.jl | 12 +++- 8 files changed, 32 insertions(+), 258 deletions(-) delete mode 100644 test/amd/qr.jl delete mode 100644 test/cuda/qr.jl diff --git a/Project.toml b/Project.toml index d7044a33..01252e07 100644 --- a/Project.toml +++ b/Project.toml @@ -30,6 +30,7 @@ GenericLinearAlgebra = "0.3.19" GenericSchur = "0.5.6" JET = "0.9, 0.10" LinearAlgebra = "1" +Random = "1" SafeTestsets = "0.1" StableRNGs = "1" Test = "1" @@ -43,6 +44,7 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -50,4 +52,4 @@ TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["Aqua", "JET", "SafeTestsets", "Test", "TestExtras", "ChainRulesCore", "ChainRulesTestUtils", "StableRNGs", "Zygote", "CUDA", "AMDGPU", "GenericLinearAlgebra", "GenericSchur"] +test = ["Aqua", "JET", "SafeTestsets", "Test", "TestExtras", "ChainRulesCore", "ChainRulesTestUtils", "StableRNGs", "Zygote", "CUDA", "AMDGPU", "GenericLinearAlgebra", "GenericSchur", "Random"] diff --git a/src/implementations/qr.jl b/src/implementations/qr.jl index a3c9d1f5..c3937b27 100644 --- a/src/implementations/qr.jl +++ b/src/implementations/qr.jl @@ -270,10 +270,12 @@ function _gpu_unmqr!( end function _gpu_qr!( - A::AbstractMatrix, Q::AbstractMatrix, R::AbstractMatrix; positive = false, blocksize = 1 + A::AbstractMatrix, Q::AbstractMatrix, R::AbstractMatrix; pivoted = false, positive = false, blocksize = 1 ) blocksize > 1 && throw(ArgumentError("CUSOLVER/ROCSOLVER does not provide a blocked implementation for a QR decomposition")) + pivoted && + throw(ArgumentError("CUSOLVER/ROCSOLVER does not provide a pivoted implementation for a QR decomposition")) m, n = size(A) minmn = min(m, n) computeR = length(R) > 0 @@ -309,10 +311,12 @@ function _gpu_qr!( end function _gpu_qr_null!( - A::AbstractMatrix, N::AbstractMatrix; positive = false, blocksize = 1 + A::AbstractMatrix, N::AbstractMatrix; positive = false, blocksize = 1, pivoted = false ) blocksize > 1 && throw(ArgumentError("CUSOLVER/ROCSOLVER does not provide a blocked implementation for a QR decomposition")) + pivoted && + throw(ArgumentError("CUSOLVER/ROCSOLVER does not provide a pivoted implementation for a QR decomposition")) m, n = size(A) minmn = min(m, n) fill!(N, zero(eltype(N))) diff --git a/test/amd/qr.jl b/test/amd/qr.jl deleted file mode 100644 index 542f5858..00000000 --- a/test/amd/qr.jl +++ /dev/null @@ -1,123 +0,0 @@ -using MatrixAlgebraKit -using MatrixAlgebraKit: diagview -using Test -using TestExtras -using StableRNGs -using AMDGPU - -include(joinpath("..", "utilities.jl")) - -eltypes = (Float32, Float64, ComplexF32, ComplexF64) - -@testset "qr_compact! and qr_null! for T = $T" for T in eltypes - rng = StableRNG(123) - m = 54 - for n in (37, m, 63) - minmn = min(m, n) - A = ROCArray(randn(rng, T, m, n)) - Q, R = @constinferred qr_compact(A) - @test Q isa ROCMatrix{T} && size(Q) == (m, minmn) - @test R isa ROCMatrix{T} && size(R) == (minmn, n) - @test Q * R ≈ A - N = @constinferred qr_null(A) - @test N isa ROCMatrix{T} && size(N) == (m, m - minmn) - @test isapproxone(Q' * Q) - @test maximum(abs, A' * N) < eps(real(T))^(2 / 3) - @test isapproxone(N' * N) - - Ac = similar(A) - Q2, R2 = @constinferred qr_compact!(copy!(Ac, A), (Q, R)) - @test Q2 === Q - @test R2 === R - N2 = @constinferred qr_null!(copy!(Ac, A), N) - @test N2 === N - - # noR - Q2 = similar(Q) - noR = similar(A, minmn, 0) - qr_compact!(copy!(Ac, A), (Q2, noR)) - @test Q == Q2 - - # positive - qr_compact!(copy!(Ac, A), (Q, R); positive = true) - @test Q * R ≈ A - @test isapproxone(Q' * Q) - @test all(>=(zero(real(T))), real(diagview(R))) - qr_compact!(copy!(Ac, A), (Q2, noR); positive = true) - @test Q == Q2 - - # explicit blocksize - qr_compact!(copy!(Ac, A), (Q, R); blocksize = 1) - @test Q * R ≈ A - @test isapproxone(Q' * Q) - qr_compact!(copy!(Ac, A), (Q2, noR); blocksize = 1) - @test Q == Q2 - qr_compact!(copy!(Ac, A), (Q2, noR); blocksize = 1) - qr_null!(copy!(Ac, A), N; blocksize = 1) - @test maximum(abs, A' * N) < eps(real(T))^(2 / 3) - @test isapproxone(N' * N) - if n <= m - qr_compact!(copy!(Q2, A), (Q2, noR); blocksize = 1) # in-place Q - @test Q ≈ Q2 - # these do not work because of the in-place Q - @test_throws ArgumentError qr_compact!(copy!(Q2, A), (Q2, R2)) - @test_throws ArgumentError qr_compact!(copy!(Q2, A), (Q2, noR); positive = true) - end - # no blocked CUDA - @test_throws ArgumentError qr_compact!(copy!(Ac, A), (Q2, R); blocksize = 8) - end -end - -@testset "qr_full! for T = $T" for T in eltypes - rng = StableRNG(123) - m = 63 - for n in (37, m, 63) - minmn = min(m, n) - A = ROCArray(randn(rng, T, m, n)) - Q, R = qr_full(A) - @test Q isa ROCMatrix{T} && size(Q) == (m, m) - @test R isa ROCMatrix{T} && size(R) == (m, n) - @test Q * R ≈ A - @test isapproxone(Q' * Q) - - Ac = similar(A) - Q2 = similar(Q) - noR = similar(A, m, 0) - Q2, R2 = @constinferred qr_full!(copy!(Ac, A), (Q, R)) - @test Q2 === Q - @test R2 === R - @test Q * R ≈ A - @test isapproxone(Q' * Q) - qr_full!(copy!(Ac, A), (Q2, noR)) - @test Q == Q2 - - # noR - noR = similar(A, m, 0) - Q2 = similar(Q) - qr_full!(copy!(Ac, A), (Q2, noR)) - @test Q == Q2 - - # positive - qr_full!(copy!(Ac, A), (Q, R); positive = true) - @test Q * R ≈ A - @test isapproxone(Q' * Q) - @test all(>=(zero(real(T))), real(diagview(R))) - qr_full!(copy!(Ac, A), (Q2, noR); positive = true) - @test Q == Q2 - - # explicit blocksize - qr_full!(copy!(Ac, A), (Q, R); blocksize = 1) - @test Q * R ≈ A - @test isapproxone(Q' * Q) - qr_full!(copy!(Ac, A), (Q2, noR); blocksize = 1) - @test Q == Q2 - if n == m - qr_full!(copy!(Q2, A), (Q2, noR); blocksize = 1) # in-place Q - @test Q ≈ Q2 - @test_throws ArgumentError qr_full!(copy!(Q2, A), (Q2, R2)) - @test_throws ArgumentError qr_full!(copy!(Q2, A), (Q2, noR); positive = true) - end - # no blocked CUDA - @test_throws ArgumentError qr_full!(copy!(Ac, A), (Q, R); blocksize = 8) - end -end diff --git a/test/cuda/qr.jl b/test/cuda/qr.jl deleted file mode 100644 index 93e4b748..00000000 --- a/test/cuda/qr.jl +++ /dev/null @@ -1,123 +0,0 @@ -using MatrixAlgebraKit -using MatrixAlgebraKit: diagview -using Test -using TestExtras -using StableRNGs -using CUDA - -include(joinpath("..", "utilities.jl")) - -eltypes = (Float32, Float64, ComplexF32, ComplexF64) - -@testset "qr_compact! and qr_null! for T = $T" for T in eltypes - rng = StableRNG(123) - m = 54 - for n in (37, m, 63) - minmn = min(m, n) - A = CuArray(randn(rng, T, m, n)) - Q, R = @constinferred qr_compact(A) - @test Q isa CuMatrix{T} && size(Q) == (m, minmn) - @test R isa CuMatrix{T} && size(R) == (minmn, n) - @test Q * R ≈ A - N = @constinferred qr_null(A) - @test N isa CuMatrix{T} && size(N) == (m, m - minmn) - @test isapproxone(Q' * Q) - @test maximum(abs, A' * N) < eps(real(T))^(2 / 3) - @test isapproxone(N' * N) - - Ac = similar(A) - Q2, R2 = @constinferred qr_compact!(copy!(Ac, A), (Q, R)) - @test Q2 === Q - @test R2 === R - N2 = @constinferred qr_null!(copy!(Ac, A), N) - @test N2 === N - - # noR - Q2 = similar(Q) - noR = similar(A, minmn, 0) - qr_compact!(copy!(Ac, A), (Q2, noR)) - @test Q == Q2 - - # positive - qr_compact!(copy!(Ac, A), (Q, R); positive = true) - @test Q * R ≈ A - @test isapproxone(Q' * Q) - @test all(>=(zero(real(T))), real(diagview(R))) - qr_compact!(copy!(Ac, A), (Q2, noR); positive = true) - @test Q == Q2 - - # explicit blocksize - qr_compact!(copy!(Ac, A), (Q, R); blocksize = 1) - @test Q * R ≈ A - @test isapproxone(Q' * Q) - qr_compact!(copy!(Ac, A), (Q2, noR); blocksize = 1) - @test Q == Q2 - qr_compact!(copy!(Ac, A), (Q2, noR); blocksize = 1) - qr_null!(copy!(Ac, A), N; blocksize = 1) - @test maximum(abs, A' * N) < eps(real(T))^(2 / 3) - @test isapproxone(N' * N) - if n <= m - qr_compact!(copy!(Q2, A), (Q2, noR); blocksize = 1) # in-place Q - @test Q ≈ Q2 - # these do not work because of the in-place Q - @test_throws ArgumentError qr_compact!(copy!(Q2, A), (Q2, R2)) - @test_throws ArgumentError qr_compact!(copy!(Q2, A), (Q2, noR); positive = true) - end - # no blocked CUDA - @test_throws ArgumentError qr_compact!(copy!(Ac, A), (Q2, R); blocksize = 8) - end -end - -@testset "qr_full! for T = $T" for T in eltypes - rng = StableRNG(123) - m = 63 - for n in (37, m, 63) - minmn = min(m, n) - A = CuArray(randn(rng, T, m, n)) - Q, R = qr_full(A) - @test Q isa CuMatrix{T} && size(Q) == (m, m) - @test R isa CuMatrix{T} && size(R) == (m, n) - @test Q * R ≈ A - @test isapproxone(Q' * Q) - - Ac = similar(A) - Q2 = similar(Q) - noR = similar(A, m, 0) - Q2, R2 = @constinferred qr_full!(copy!(Ac, A), (Q, R)) - @test Q2 === Q - @test R2 === R - @test Q * R ≈ A - @test isapproxone(Q' * Q) - qr_full!(copy!(Ac, A), (Q2, noR)) - @test Q == Q2 - - # noR - noR = similar(A, m, 0) - Q2 = similar(Q) - qr_full!(copy!(Ac, A), (Q2, noR)) - @test Q == Q2 - - # positive - qr_full!(copy!(Ac, A), (Q, R); positive = true) - @test Q * R ≈ A - @test isapproxone(Q' * Q) - @test all(>=(zero(real(T))), real(diagview(R))) - qr_full!(copy!(Ac, A), (Q2, noR); positive = true) - @test Q == Q2 - - # explicit blocksize - qr_full!(copy!(Ac, A), (Q, R); blocksize = 1) - @test Q * R ≈ A - @test isapproxone(Q' * Q) - qr_full!(copy!(Ac, A), (Q2, noR); blocksize = 1) - @test Q == Q2 - if n == m - qr_full!(copy!(Q2, A), (Q2, noR); blocksize = 1) # in-place Q - @test Q ≈ Q2 - @test_throws ArgumentError qr_full!(copy!(Q2, A), (Q2, R2)) - @test_throws ArgumentError qr_full!(copy!(Q2, A), (Q2, noR); positive = true) - end - # no blocked CUDA - @test_throws ArgumentError qr_full!(copy!(Ac, A), (Q, R); blocksize = 8) - end -end diff --git a/test/qr.jl b/test/qr.jl index ed8c255f..d7f974be 100644 --- a/test/qr.jl +++ b/test/qr.jl @@ -3,6 +3,7 @@ using Test using TestExtras using StableRNGs using LinearAlgebra: diag, I, Diagonal +using CUDA, AMDGPU BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) GenericFloats = (Float16, BigFloat, Complex{BigFloat}) @@ -14,6 +15,14 @@ m = 54 for T in BLASFloats, n in (37, m, 63) TestSuite.seed_rng!(123) TestSuite.test_qr(T, (m, n)) + if CUDA.functional() + TestSuite.test_qr(CuMatrix{T}, (m, n); test_pivoted = false, test_blocksize = false) + TestSuite.test_qr(Diagonal{T, CuVector{T}}, m; test_pivoted = false, test_blocksize = false) + end + if AMDGPU.functional() + TestSuite.test_qr(ROCMatrix{T}, (m, n); test_pivoted = false, test_blocksize = false) + TestSuite.test_qr(Diagonal{T, ROCVector{T}}, m; test_pivoted = false, test_blocksize = false) + end end for T in (BLASFloats..., GenericFloats...) AT = Diagonal{T, Vector{T}} diff --git a/test/runtests.jl b/test/runtests.jl index ec255538..2345fa80 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -57,9 +57,6 @@ end using CUDA if CUDA.functional() - @safetestset "CUDA QR" begin - include("cuda/qr.jl") - end @safetestset "CUDA LQ" begin include("cuda/lq.jl") end @@ -85,9 +82,6 @@ end using AMDGPU if AMDGPU.functional() - @safetestset "AMDGPU QR" begin - include("amd/qr.jl") - end @safetestset "AMDGPU LQ" begin include("amd/lq.jl") end diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index c75abe1f..b7adb959 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -13,6 +13,7 @@ using MatrixAlgebraKit using MatrixAlgebraKit: diagview using LinearAlgebra: Diagonal, norm, istriu using Random, StableRNGs +using AMDGPU, CUDA const tests = Dict() @@ -33,7 +34,11 @@ seed_rng!(seed) = Random.seed!(rng, seed) instantiate_matrix(::Type{T}, size) where {T <: Number} = randn(rng, T, size) instantiate_matrix(::Type{AT}, size) where {AT <: Array} = randn(rng, eltype(AT), size) +instantiate_matrix(::Type{AT}, size) where {AT <: CuArray} = CuArray(randn(rng, eltype(AT), size)) +instantiate_matrix(::Type{AT}, size) where {AT <: ROCArray} = ROCArray(randn(rng, eltype(AT), size)) instantiate_matrix(::Type{AT}, size) where {AT <: Diagonal} = Diagonal(randn(rng, eltype(AT), size)) +instantiate_matrix(::Type{AT}, size) where {T, AT <: Diagonal{T, <:CuVector}} = Diagonal(CuArray(randn(rng, eltype(AT), size))) +instantiate_matrix(::Type{AT}, size) where {T, AT <: Diagonal{T, <:ROCVector}} = Diagonal(ROCArray(randn(rng, eltype(AT), size))) precision(::Type{T}) where {T <: Number} = sqrt(eps(real(T))) precision(::Type{T}) where {T} = precision(eltype(T)) diff --git a/test/testsuite/qr.jl b/test/testsuite/qr.jl index c354c2e4..d381c046 100644 --- a/test/testsuite/qr.jl +++ b/test/testsuite/qr.jl @@ -1,9 +1,15 @@ function test_qr(T::Type, sz; kwargs...) summary_str = testargs_summary(T, sz) return @testset "qr $summary_str" begin - test_qr_compact(T, sz; kwargs...) - test_qr_full(T, sz; kwargs...) - test_qr_null(T, sz; kwargs...) + if T <: Union{CuArray, ROCArray} + test_qr_compact(T, sz; pivoted = false, kwargs...) + test_qr_full(T, sz; pivoted = false, kwargs...) + test_qr_null(T, sz; pivoted = false, kwargs...) + else + test_qr_compact(T, sz; kwargs...) + test_qr_full(T, sz; kwargs...) + test_qr_null(T, sz; kwargs...) + end end end From cf7839a85c78e9555ea69c47bcdf46c2760555ad Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Thu, 13 Nov 2025 09:58:10 +0100 Subject: [PATCH 11/14] LQ for TestSuite --- test/lq.jl | 21 +++++ test/testsuite/TestSuite.jl | 6 +- test/testsuite/lq.jl | 160 ++++++++++++++++++++++++++++++++++++ test/testsuite/qr.jl | 12 +-- 4 files changed, 189 insertions(+), 10 deletions(-) create mode 100644 test/testsuite/lq.jl diff --git a/test/lq.jl b/test/lq.jl index 8de5a582..1efd239e 100644 --- a/test/lq.jl +++ b/test/lq.jl @@ -4,10 +4,30 @@ using TestExtras using StableRNGs using LinearAlgebra: diag, I, Diagonal using MatrixAlgebraKit: LQViaTransposedQR, LAPACK_HouseholderQR +using CUDA, AMDGPU BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) GenericFloats = (Float16, BigFloat, Complex{BigFloat}) +m = 54 +for T in BLASFloats, n in (37, m, 63) + TestSuite.seed_rng!(123) + TestSuite.test_lq(T, (m, n); test_pivoted = false) + if CUDA.functional() + TestSuite.test_lq(CuMatrix{T}, (m, n); test_pivoted = false, test_blocksize = false) + TestSuite.test_lq(Diagonal{T, CuVector{T}}, m; test_pivoted = false, test_blocksize = false) + end + if AMDGPU.functional() + TestSuite.test_lq(ROCMatrix{T}, (m, n); test_pivoted = false, test_blocksize = false) + TestSuite.test_lq(Diagonal{T, ROCVector{T}}, m; test_pivoted = false, test_blocksize = false) + end +end +for T in (BLASFloats..., GenericFloats...) + AT = Diagonal{T, Vector{T}} + TestSuite.test_lq(AT, m; test_pivoted = false, test_blocksize = false) +end + +#= @testset "lq_compact! for T = $T" for T in BLASFloats rng = StableRNG(123) m = 54 @@ -253,3 +273,4 @@ end @test N isa AbstractMatrix{T} && size(N) == (0, m) end end +=# diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index b7adb959..b62f0c2d 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -11,7 +11,7 @@ module TestSuite using Test, TestExtras using MatrixAlgebraKit using MatrixAlgebraKit: diagview -using LinearAlgebra: Diagonal, norm, istriu +using LinearAlgebra: Diagonal, norm, istriu, istril using Random, StableRNGs using AMDGPU, CUDA @@ -55,11 +55,15 @@ end isleftnull(N, A; atol::Real = 0, rtol::Real = precision(eltype(A))) = isapprox(norm(A' * N), 0; atol = max(atol, norm(A) * rtol)) +isrightnull(Nᴴ, A; atol::Real = 0, rtol::Real = precision(eltype(A))) = + isapprox(norm(A * Nᴴ'), 0; atol = max(atol, norm(A) * rtol)) + # TODO: actually make this a test macro testinferred(ex) return esc(:(@inferred $ex)) end include("qr.jl") +include("lq.jl") end diff --git a/test/testsuite/lq.jl b/test/testsuite/lq.jl new file mode 100644 index 00000000..5ad41c3b --- /dev/null +++ b/test/testsuite/lq.jl @@ -0,0 +1,160 @@ +function test_lq(T::Type, sz; kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "lq $summary_str" begin + test_lq_compact(T, sz; kwargs...) + test_lq_full(T, sz; kwargs...) + test_lq_null(T, sz; kwargs...) + end +end + +function test_lq_compact( + T::Type, sz; + test_positive = true, test_pivoted = true, test_blocksize = true, + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "lq_compact! $summary_str" begin + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + + # does the elementary functionality work + L, Q = @testinferred lq_compact(A) + @test L * Q ≈ A + @test isisometric(Q; side = :right, atol, rtol) + @test istril(L) + @test A == Ac + + # can I pass in outputs? + L2, Q2 = @testinferred lq_compact!(deepcopy(A), (L, Q)) + @test L2 * Q2 ≈ A + @test isisometric(Q2; side = :right, atol, rtol) + @test istril(L2) + + # do we support `positive = true`? + if test_positive + Lpos, Qpos = @testinferred lq_compact(A; positive = true) + @test Lpos * Qpos ≈ A + @test isisometric(Qpos; side = :right, atol, rtol) + @test istril(Lpos) + @test has_positive_diagonal(Lpos) + else + @test_throws Exception lq_compact(A; positive = true) + end + + # do we support `pivoted = true`? + if test_pivoted + Lpiv, Qpiv = @testinferred lq_compact(A; pivoted = true) + @test Lpiv * Qpiv ≈ A + @test isisometric(Qpos; side = :right, atol, rtol) + else + @test_throws Exception lq_compact(A; pivoted = true) + end + + # do we support `blocksize = Int`? + if test_blocksize + Lblocked, Qblocked = @testinferred lq_compact(A; blocksize = 2) + @test Lblocked * Qblocked ≈ A + @test isisometric(Qblocked; side = :right, atol, rtol) + else + @test_throws Exception lq_compact(A; blocksize = 2) + end + end +end + +function test_lq_full( + T::Type, sz; + test_positive = true, test_pivoted = true, test_blocksize = true, + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "lq_full! $summary_str" begin + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + + # does the elementary functionality work + L, Q = @testinferred lq_full(A) + @test L * Q ≈ A + @test isunitary(Q; atol, rtol) + @test istril(L) + @test A == Ac + + # can I pass in outputs? + L2, Q2 = @testinferred lq_full!(deepcopy(A), (L, Q)) + @test L2 * Q2 ≈ A + @test isunitary(Q2; atol, rtol) + @test istril(L2) + + # do we support `positive = true`? + if test_positive + Lpos, Qpos = @testinferred lq_full(A; positive = true) + @test Lpos * Qpos ≈ A + @test isunitary(Qpos; atol, rtol) + @test istril(Lpos) + @test has_positive_diagonal(Lpos) + else + @test_throws Exception lq_full(A; positive = true) + end + + # do we support `pivoted = true`? + if test_pivoted + Lpiv, Qpiv = @testinferred lq_full(A; pivoted = true) + @test Lpiv * Qpiv ≈ A + @test isunitary(Qpos; atol, rtol) + else + @test_throws Exception lq_full(A; pivoted = true) + end + + # do we support `blocksize = Int`? + if test_blocksize + Lblocked, Qblocked = @testinferred lq_full(A; blocksize = 2) + @test Lblocked * Qblocked ≈ A + @test isunitary(Qblocked; atol, rtol) + else + @test_throws Exception lq_full(A; blocksize = 2) + end + end +end + +function test_lq_null( + T::Type, sz; + test_pivoted = true, test_blocksize = true, + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "lq_null! $summary_str" begin + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + + # does the elementary functionality work + Nᴴ = @testinferred lq_null(A) + @test isrightnull(Nᴴ, A; atol, rtol) + @test isisometric(Nᴴ; side = :right, atol, rtol) + @test A == Ac + + # can I pass in outputs? + Nᴴ2 = @testinferred lq_null!(deepcopy(A), Nᴴ) + @test isrightnull(Nᴴ2, A; atol, rtol) + @test isisometric(Nᴴ2; side = :right, atol, rtol) + + # do we support `pivoted = true`? + if test_pivoted + Nᴴpiv = @testinferred lq_null(A; pivoted = true) + @test isrightnull(Nᴴpiv, A; atol, rtol) + @test isisometric(Nᴴpiv; side = :right, atol, rtol) + #else # DISABLE for now as lq_null does support pivoting... + # @test_throws Exception lq_null(A; pivoted = true) + end + + # do we support `blocksize = Int`? + if test_blocksize + Nᴴblocked = @testinferred lq_null(A; blocksize = 2) + @test isrightnull(Nᴴblocked, A; atol, rtol) + @test isisometric(Nᴴblocked; side = :right, atol, rtol) + else + @test_throws Exception lq_null(A; blocksize = 2) + end + end +end diff --git a/test/testsuite/qr.jl b/test/testsuite/qr.jl index d381c046..c354c2e4 100644 --- a/test/testsuite/qr.jl +++ b/test/testsuite/qr.jl @@ -1,15 +1,9 @@ function test_qr(T::Type, sz; kwargs...) summary_str = testargs_summary(T, sz) return @testset "qr $summary_str" begin - if T <: Union{CuArray, ROCArray} - test_qr_compact(T, sz; pivoted = false, kwargs...) - test_qr_full(T, sz; pivoted = false, kwargs...) - test_qr_null(T, sz; pivoted = false, kwargs...) - else - test_qr_compact(T, sz; kwargs...) - test_qr_full(T, sz; kwargs...) - test_qr_null(T, sz; kwargs...) - end + test_qr_compact(T, sz; kwargs...) + test_qr_full(T, sz; kwargs...) + test_qr_null(T, sz; kwargs...) end end From ff4fc12439543d6f09b1b63686fd84830308ac6f Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Wed, 19 Nov 2025 05:49:58 -0500 Subject: [PATCH 12/14] Move polar to testsuite --- test/polar.jl | 95 ++++++++----------------------------- test/testsuite/TestSuite.jl | 1 + test/testsuite/polar.jl | 81 +++++++++++++++++++++++++++++++ 3 files changed, 103 insertions(+), 74 deletions(-) create mode 100644 test/testsuite/polar.jl diff --git a/test/polar.jl b/test/polar.jl index 8087149e..6ac7b321 100644 --- a/test/polar.jl +++ b/test/polar.jl @@ -4,80 +4,27 @@ using TestExtras using StableRNGs using LinearAlgebra: LinearAlgebra, I, isposdef -@testset "left_polar! for T = $T" for T in (Float32, Float64, ComplexF32, ComplexF64) - rng = StableRNG(123) - m = 54 - @testset "size ($m, $n)" for n in (37, m) - k = min(m, n) - if LinearAlgebra.LAPACK.version() < v"3.12.0" - svdalgs = (LAPACK_DivideAndConquer(), LAPACK_QRIteration(), LAPACK_Bisection()) - else - svdalgs = (LAPACK_DivideAndConquer(), LAPACK_QRIteration(), LAPACK_Bisection(), LAPACK_Jacobi()) - end - algs = (PolarViaSVD.(svdalgs)..., PolarNewton()) - @testset "algorithm $alg" for alg in algs - A = randn(rng, T, m, n) - - W, P = left_polar(A; alg) - @test W isa Matrix{T} && size(W) == (m, n) - @test P isa Matrix{T} && size(P) == (n, n) - @test W * P ≈ A - @test isisometric(W) - @test isposdef(P) - - Ac = similar(A) - W2, P2 = @constinferred left_polar!(copy!(Ac, A), (W, P), alg) - @test W2 === W - @test P2 === P - @test W * P ≈ A - @test isisometric(W) - @test isposdef(P) - - noP = similar(P, (0, 0)) - W2, P2 = @constinferred left_polar!(copy!(Ac, A), (W, noP), alg) - @test P2 === noP - @test W2 === W - @test isisometric(W) - P = W' * A # compute P explicitly to verify W correctness - @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) - @test isposdef(project_hermitian!(P)) - end +BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) +GenericFloats = (Float16, BigFloat, Complex{BigFloat}) + +@isdefined(TestSuite) || include("testsuite/TestSuite.jl") +using .TestSuite + +m = 54 +for T in BLASFloats, n in (37, m, 63) + TestSuite.seed_rng!(123) + TestSuite.test_polar(T, (m, n)) + if CUDA.functional() + TestSuite.test_polar(CuMatrix{T}, (m, n); test_pivoted = false, test_blocksize = false) + TestSuite.test_polar(Diagonal{T, CuVector{T}}, m; test_pivoted = false, test_blocksize = false) end -end - -@testset "right_polar! for T = $T" for T in (Float32, Float64, ComplexF32, ComplexF64) - rng = StableRNG(123) - n = 54 - @testset "size ($m, $n)" for m in (37, n) - k = min(m, n) - svdalgs = (LAPACK_DivideAndConquer(), LAPACK_QRIteration(), LAPACK_Bisection()) - algs = (PolarViaSVD.(svdalgs)..., PolarNewton()) - @testset "algorithm $alg" for alg in algs - A = randn(rng, T, m, n) - - P, Wᴴ = right_polar(A; alg) - @test Wᴴ isa Matrix{T} && size(Wᴴ) == (m, n) - @test P isa Matrix{T} && size(P) == (m, m) - @test P * Wᴴ ≈ A - @test isisometric(Wᴴ; side = :right) - @test isposdef(P) - - Ac = similar(A) - P2, Wᴴ2 = @constinferred right_polar!(copy!(Ac, A), (P, Wᴴ), alg) - @test P2 === P - @test Wᴴ2 === Wᴴ - @test P * Wᴴ ≈ A - @test isisometric(Wᴴ; side = :right) - @test isposdef(P) - - noP = similar(P, (0, 0)) - P2, Wᴴ2 = @constinferred right_polar!(copy!(Ac, A), (noP, Wᴴ), alg) - @test P2 === noP - @test Wᴴ2 === Wᴴ - @test isisometric(Wᴴ; side = :right) - P = A * Wᴴ' # compute P explicitly to verify W correctness - @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) - @test isposdef(project_hermitian!(P)) - end + if AMDGPU.functional() + TestSuite.test_polar(ROCMatrix{T}, (m, n); test_pivoted = false, test_blocksize = false) + TestSuite.test_polar(Diagonal{T, ROCVector{T}}, m; test_pivoted = false, test_blocksize = false) end end +for T in (BLASFloats..., GenericFloats...) + AT = Diagonal{T, Vector{T}} + TestSuite.test_polar(AT, m; test_pivoted = false, test_blocksize = false) +end + diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index b62f0c2d..765d7d25 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -65,5 +65,6 @@ end include("qr.jl") include("lq.jl") +include("polar.jl") end diff --git a/test/testsuite/polar.jl b/test/testsuite/polar.jl new file mode 100644 index 00000000..67eece8b --- /dev/null +++ b/test/testsuite/polar.jl @@ -0,0 +1,81 @@ +function test_polar(T::Type, sz; kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "polar $summary_str" begin + test_left_polar(T, sz; kwargs...) + test_right_polar(T, sz; kwargs...) + end +end + +function test_left_polar( + T::Type, sz; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "left_polar! $summary_str" begin + algs = (PolarViaSVD(), PolarNewton()) + @testset "algorithm $alg" for alg in algs + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + W, P = left_polar(A; alg) + @test eltype(W) == eltype(A) && size(W) == (size(A, 1), size(A, 2)) + @test eltype(P) == eltype(A) && size(P) == (size(A, 2), size(A, 2)) + @test W * P ≈ A + @test isisometric(W) + @test isposdef(P) + + W2, P2 = @constinferred left_polar!(Ac, (W, P), alg) + @test W2 === W + @test P2 === P + @test W * P ≈ A + @test isisometric(W) + @test isposdef(P) + + noP = similar(P, (0, 0)) + W2, P2 = @constinferred left_polar!(copy!(Ac, A), (W, noP), alg) + @test P2 === noP + @test W2 === W + @test isisometric(W) + P = W' * A # compute P explicitly to verify W correctness + @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) + @test isposdef(project_hermitian!(P)) + end + end +end + +function test_right_polar( + T::Type, sz; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "right_polar! $summary_str" begin + algs = (PolarViaSVD(), PolarNewton()) + @testset "algorithm $alg" for alg in algs + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + P, Wᴴ = right_polar(A; alg) + @test eltype(Wᴴ) == eltype(A) && size(Wᴴ) == (size(A, 1), size(A, 2)) + @test eltype(P) == eltype(A) && size(P) == (size(A, 1), size(A, 1)) + @test P * Wᴴ ≈ A + @test isisometric(Wᴴ; side = :right) + @test isposdef(P) + + P2, Wᴴ2 = @constinferred right_polar!(Ac, (P, Wᴴ), alg) + @test P2 === P + @test Wᴴ2 === Wᴴ + @test P * Wᴴ ≈ A + @test isisometric(Wᴴ; side = :right) + @test isposdef(P) + + noP = similar(P, (0, 0)) + P2, Wᴴ2 = @constinferred right_polar!(copy!(Ac, A), (noP, Wᴴ), alg) + @test P2 === noP + @test Wᴴ2 === Wᴴ + @test isisometric(Wᴴ; side = :right) + P = A * Wᴴ' # compute P explicitly to verify W correctness + @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) + @test isposdef(project_hermitian!(P)) + end + end +end From 2d8cb84b8defd28fdd04de3e1b145b60455afa6b Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Wed, 19 Nov 2025 05:51:20 -0500 Subject: [PATCH 13/14] Purge QR/LQ --- test/amd/lq.jl | 117 ---------------------- test/cuda/lq.jl | 117 ---------------------- test/lq.jl | 248 ----------------------------------------------- test/qr.jl | 216 ----------------------------------------- test/runtests.jl | 6 -- 5 files changed, 704 deletions(-) delete mode 100644 test/amd/lq.jl delete mode 100644 test/cuda/lq.jl diff --git a/test/amd/lq.jl b/test/amd/lq.jl deleted file mode 100644 index a3dc8e53..00000000 --- a/test/amd/lq.jl +++ /dev/null @@ -1,117 +0,0 @@ -using MatrixAlgebraKit -using MatrixAlgebraKit: diagview -using Test -using TestExtras -using StableRNGs -using AMDGPU - -include(joinpath("..", "utilities.jl")) - -@testset "lq_compact! for T = $T" for T in (Float32, Float64, ComplexF32, ComplexF64) - rng = StableRNG(123) - m = 54 - for n in (37, m, 63) - minmn = min(m, n) - A = ROCArray(randn(rng, T, m, n)) - L, Q = @constinferred lq_compact(A) - @test L isa ROCMatrix{T} && size(L) == (m, minmn) - @test Q isa ROCMatrix{T} && size(Q) == (minmn, n) - @test L * Q ≈ A - @test isapproxone(Q * Q') - Nᴴ = @constinferred lq_null(A) - @test Nᴴ isa ROCMatrix{T} && size(Nᴴ) == (n - minmn, n) - @test maximum(abs, A * Nᴴ') < eps(real(T))^(2 / 3) - @test isapproxone(Nᴴ * Nᴴ') - - Ac = similar(A) - L2, Q2 = @constinferred lq_compact!(copy!(Ac, A), (L, Q)) - @test L2 === L - @test Q2 === Q - Nᴴ2 = @constinferred lq_null!(copy!(Ac, A), Nᴴ) - @test Nᴴ2 === Nᴴ - - # noL - noL = similar(A, 0, minmn) - Q2 = similar(Q) - lq_compact!(copy!(Ac, A), (noL, Q2)) - @test Q == Q2 - - # positive - lq_compact!(copy!(Ac, A), (L, Q); positive = true) - @test L * Q ≈ A - @test isapproxone(Q * Q') - @test all(>=(zero(real(T))), real(diagview(L))) - lq_compact!(copy!(Ac, A), (noL, Q2); positive = true) - @test Q == Q2 - - # explicit blocksize - lq_compact!(copy!(Ac, A), (L, Q); blocksize = 1) - @test L * Q ≈ A - @test isapproxone(Q * Q') - lq_compact!(copy!(Ac, A), (noL, Q2); blocksize = 1) - @test Q == Q2 - lq_null!(copy!(Ac, A), Nᴴ; blocksize = 1) - @test maximum(abs, A * Nᴴ') < eps(real(T))^(2 / 3) - @test isapproxone(Nᴴ * Nᴴ') - if m <= n - lq_compact!(copy!(Q2, A), (noL, Q2); blocksize = 1) # in-place Q - @test Q ≈ Q2 - # these do not work because of the in-place Q - @test_throws ArgumentError lq_compact!(copy!(Q2, A), (L, Q2); blocksize = 1) - @test_throws ArgumentError lq_compact!(copy!(Q2, A), (noL, Q2); positive = true) - end - # no blocked CUDA - @test_throws ArgumentError lq_compact!(copy!(Q2, A), (L, Q2); blocksize = 8) - end -end - -@testset "lq_full! for T = $T" for T in (Float32, Float64, ComplexF32, ComplexF64) - rng = StableRNG(123) - m = 54 - for n in (37, m, 63) - minmn = min(m, n) - A = ROCArray(randn(rng, T, m, n)) - L, Q = lq_full(A) - @test L isa ROCMatrix{T} && size(L) == (m, n) - @test Q isa ROCMatrix{T} && size(Q) == (n, n) - @test L * Q ≈ A - @test isapproxone(Q * Q') - - Ac = similar(A) - L2, Q2 = @constinferred lq_full!(copy!(Ac, A), (L, Q)) - @test L2 === L - @test Q2 === Q - @test L * Q ≈ A - @test isapproxone(Q * Q') - - # noL - noL = similar(A, 0, n) - Q2 = similar(Q) - lq_full!(copy!(Ac, A), (noL, Q2)) - @test Q == Q2 - - # positive - lq_full!(copy!(Ac, A), (L, Q); positive = true) - @test L * Q ≈ A - @test isapproxone(Q * Q') - @test all(>=(zero(real(T))), real(diagview(L))) - lq_full!(copy!(Ac, A), (noL, Q2); positive = true) - @test Q == Q2 - - # explicit blocksize - lq_full!(copy!(Ac, A), (L, Q); blocksize = 1) - @test L * Q ≈ A - @test isapproxone(Q * Q') - lq_full!(copy!(Ac, A), (noL, Q2); blocksize = 1) - @test Q == Q2 - if n == m - lq_full!(copy!(Q2, A), (noL, Q2); blocksize = 1) # in-place Q - @test Q ≈ Q2 - # these do not work because of the in-place Q - @test_throws ArgumentError lq_full!(copy!(Q2, A), (L, Q2); blocksize = 1) - @test_throws ArgumentError lq_full!(copy!(Q2, A), (noL, Q2); positive = true) - end - # no blocked CUDA - @test_throws ArgumentError lq_full!(copy!(Ac, A), (L, Q); blocksize = 8) - end -end diff --git a/test/cuda/lq.jl b/test/cuda/lq.jl deleted file mode 100644 index 9a162c6b..00000000 --- a/test/cuda/lq.jl +++ /dev/null @@ -1,117 +0,0 @@ -using MatrixAlgebraKit -using MatrixAlgebraKit: diagview -using Test -using TestExtras -using StableRNGs -using CUDA - -include(joinpath("..", "utilities.jl")) - -@testset "lq_compact! for T = $T" for T in (Float32, Float64, ComplexF32, ComplexF64) - rng = StableRNG(123) - m = 54 - for n in (37, m, 63) - minmn = min(m, n) - A = CuArray(randn(rng, T, m, n)) - L, Q = @constinferred lq_compact(A) - @test L isa CuMatrix{T} && size(L) == (m, minmn) - @test Q isa CuMatrix{T} && size(Q) == (minmn, n) - @test L * Q ≈ A - @test isapproxone(Q * Q') - Nᴴ = @constinferred lq_null(A) - @test Nᴴ isa CuMatrix{T} && size(Nᴴ) == (n - minmn, n) - @test maximum(abs, A * Nᴴ') < eps(real(T))^(2 / 3) - @test isapproxone(Nᴴ * Nᴴ') - - Ac = similar(A) - L2, Q2 = @constinferred lq_compact!(copy!(Ac, A), (L, Q)) - @test L2 === L - @test Q2 === Q - Nᴴ2 = @constinferred lq_null!(copy!(Ac, A), Nᴴ) - @test Nᴴ2 === Nᴴ - - # noL - noL = similar(A, 0, minmn) - Q2 = similar(Q) - lq_compact!(copy!(Ac, A), (noL, Q2)) - @test Q == Q2 - - # positive - lq_compact!(copy!(Ac, A), (L, Q); positive = true) - @test L * Q ≈ A - @test isapproxone(Q * Q') - @test all(>=(zero(real(T))), real(diagview(L))) - lq_compact!(copy!(Ac, A), (noL, Q2); positive = true) - @test Q == Q2 - - # explicit blocksize - lq_compact!(copy!(Ac, A), (L, Q); blocksize = 1) - @test L * Q ≈ A - @test isapproxone(Q * Q') - lq_compact!(copy!(Ac, A), (noL, Q2); blocksize = 1) - @test Q == Q2 - lq_null!(copy!(Ac, A), Nᴴ; blocksize = 1) - @test maximum(abs, A * Nᴴ') < eps(real(T))^(2 / 3) - @test isapproxone(Nᴴ * Nᴴ') - if m <= n - lq_compact!(copy!(Q2, A), (noL, Q2); blocksize = 1) # in-place Q - @test Q ≈ Q2 - # these do not work because of the in-place Q - @test_throws ArgumentError lq_compact!(copy!(Q2, A), (L, Q2); blocksize = 1) - @test_throws ArgumentError lq_compact!(copy!(Q2, A), (noL, Q2); positive = true) - end - # no blocked CUDA - @test_throws ArgumentError lq_compact!(copy!(Q2, A), (L, Q2); blocksize = 8) - end -end - -@testset "lq_full! for T = $T" for T in (Float32, Float64, ComplexF32, ComplexF64) - rng = StableRNG(123) - m = 54 - for n in (37, m, 63) - minmn = min(m, n) - A = CuArray(randn(rng, T, m, n)) - L, Q = lq_full(A) - @test L isa CuMatrix{T} && size(L) == (m, n) - @test Q isa CuMatrix{T} && size(Q) == (n, n) - @test L * Q ≈ A - @test isapproxone(Q * Q') - - Ac = similar(A) - L2, Q2 = @constinferred lq_full!(copy!(Ac, A), (L, Q)) - @test L2 === L - @test Q2 === Q - @test L * Q ≈ A - @test isapproxone(Q * Q') - - # noL - noL = similar(A, 0, n) - Q2 = similar(Q) - lq_full!(copy!(Ac, A), (noL, Q2)) - @test Q == Q2 - - # positive - lq_full!(copy!(Ac, A), (L, Q); positive = true) - @test L * Q ≈ A - @test isapproxone(Q * Q') - @test all(>=(zero(real(T))), real(diagview(L))) - lq_full!(copy!(Ac, A), (noL, Q2); positive = true) - @test Q == Q2 - - # explicit blocksize - lq_full!(copy!(Ac, A), (L, Q); blocksize = 1) - @test L * Q ≈ A - @test isapproxone(Q * Q') - lq_full!(copy!(Ac, A), (noL, Q2); blocksize = 1) - @test Q == Q2 - if n == m - lq_full!(copy!(Q2, A), (noL, Q2); blocksize = 1) # in-place Q - @test Q ≈ Q2 - # these do not work because of the in-place Q - @test_throws ArgumentError lq_full!(copy!(Q2, A), (L, Q2); blocksize = 1) - @test_throws ArgumentError lq_full!(copy!(Q2, A), (noL, Q2); positive = true) - end - # no blocked CUDA - @test_throws ArgumentError lq_full!(copy!(Ac, A), (L, Q); blocksize = 8) - end -end diff --git a/test/lq.jl b/test/lq.jl index 1efd239e..eb9fd773 100644 --- a/test/lq.jl +++ b/test/lq.jl @@ -26,251 +26,3 @@ for T in (BLASFloats..., GenericFloats...) AT = Diagonal{T, Vector{T}} TestSuite.test_lq(AT, m; test_pivoted = false, test_blocksize = false) end - -#= -@testset "lq_compact! for T = $T" for T in BLASFloats - rng = StableRNG(123) - m = 54 - for n in (37, m, 63) - minmn = min(m, n) - A = randn(rng, T, m, n) - L, Q = @constinferred lq_compact(A) - @test L isa Matrix{T} && size(L) == (m, minmn) - @test Q isa Matrix{T} && size(Q) == (minmn, n) - @test L * Q ≈ A - @test isisometric(Q; side = :right) - Nᴴ = @constinferred lq_null(A) - @test Nᴴ isa Matrix{T} && size(Nᴴ) == (n - minmn, n) - @test maximum(abs, A * Nᴴ') < eps(real(T))^(2 / 3) - @test isisometric(Nᴴ; side = :right) - - Ac = similar(A) - L2, Q2 = @constinferred lq_compact!(copy!(Ac, A), (L, Q)) - @test L2 === L - @test Q2 === Q - Nᴴ2 = @constinferred lq_null!(copy!(Ac, A), Nᴴ) - @test Nᴴ2 === Nᴴ - - noL = similar(A, 0, minmn) - Q2 = similar(Q) - lq_compact!(copy!(Ac, A), (noL, Q2)) - @test Q == Q2 - - # Transposed QR algorithm - qr_alg = LAPACK_HouseholderQR() - lq_alg = LQViaTransposedQR(qr_alg) - L2, Q2 = @constinferred lq_compact!(copy!(Ac, A), (L, Q), lq_alg) - @test L2 === L - @test Q2 === Q - Nᴴ2 = @constinferred lq_null!(copy!(Ac, A), Nᴴ, lq_alg) - @test Nᴴ2 === Nᴴ - noL = similar(A, 0, minmn) - Q2 = similar(Q) - lq_compact!(copy!(Ac, A), (noL, Q2), lq_alg) - @test Q == Q2 - - # unblocked algorithm - lq_compact!(copy!(Ac, A), (L, Q); blocksize = 1) - @test L * Q ≈ A - @test isisometric(Q; side = :right) - lq_compact!(copy!(Ac, A), (noL, Q2); blocksize = 1) - @test Q == Q2 - lq_null!(copy!(Ac, A), Nᴴ; blocksize = 1) - @test maximum(abs, A * Nᴴ') < eps(real(T))^(2 / 3) - @test isisometric(Nᴴ; side = :right) - if m <= n - lq_compact!(copy!(Q2, A), (noL, Q2); blocksize = 1) # in-place Q - @test Q ≈ Q2 - @test_throws ArgumentError lq_compact!(copy!(Q2, A), (L, Q2); blocksize = 1) - @test_throws ArgumentError lq_compact!(copy!(Q2, A), (noL, Q2); positive = true) - @test_throws ArgumentError lq_compact!(copy!(Q2, A), (noL, Q2); blocksize = 8) - end - lq_compact!(copy!(Ac, A), (L, Q); blocksize = 8) - @test L * Q ≈ A - @test isisometric(Q; side = :right) - lq_compact!(copy!(Ac, A), (noL, Q2); blocksize = 8) - @test Q == Q2 - lq_null!(copy!(Ac, A), Nᴴ; blocksize = 8) - @test maximum(abs, A * Nᴴ') < eps(real(T))^(2 / 3) - @test isisometric(Nᴴ; side = :right) - @test Nᴴ * Nᴴ' ≈ I - - qr_alg = LAPACK_HouseholderQR(; blocksize = 1) - lq_alg = LQViaTransposedQR(qr_alg) - lq_compact!(copy!(Ac, A), (L, Q), lq_alg) - @test L * Q ≈ A - @test Q * Q' ≈ I - lq_compact!(copy!(Ac, A), (noL, Q2), lq_alg) - @test Q == Q2 - lq_null!(copy!(Ac, A), Nᴴ, lq_alg) - @test maximum(abs, A * Nᴴ') < eps(real(T))^(2 / 3) - @test Nᴴ * Nᴴ' ≈ I - - # pivoted - @test_throws ArgumentError lq_compact!(copy!(Ac, A), (L, Q); pivoted = true) - - # positive - lq_compact!(copy!(Ac, A), (L, Q); positive = true) - @test L * Q ≈ A - @test isisometric(Q; side = :right) - @test all(>=(zero(real(T))), real(diag(L))) - lq_compact!(copy!(Ac, A), (noL, Q2); positive = true) - @test Q == Q2 - - # positive and blocksize 1 - lq_compact!(copy!(Ac, A), (L, Q); positive = true, blocksize = 1) - @test L * Q ≈ A - @test isisometric(Q; side = :right) - @test all(>=(zero(real(T))), real(diag(L))) - lq_compact!(copy!(Ac, A), (noL, Q2); positive = true, blocksize = 1) - @test Q == Q2 - qr_alg = LAPACK_HouseholderQR(; positive = true, blocksize = 1) - lq_alg = LQViaTransposedQR(qr_alg) - lq_compact!(copy!(Ac, A), (L, Q), lq_alg) - @test L * Q ≈ A - @test Q * Q' ≈ I - @test all(>=(zero(real(T))), real(diag(L))) - lq_compact!(copy!(Ac, A), (noL, Q2), lq_alg) - @test Q == Q2 - end -end - -@testset "lq_full! for T = $T" for T in BLASFloats - rng = StableRNG(123) - m = 54 - for n in (37, m, 63) - minmn = min(m, n) - A = randn(rng, T, m, n) - L, Q = lq_full(A) - @test L isa Matrix{T} && size(L) == (m, n) - @test Q isa Matrix{T} && size(Q) == (n, n) - @test L * Q ≈ A - @test isunitary(Q) - - Ac = similar(A) - L2, Q2 = @constinferred lq_full!(copy!(Ac, A), (L, Q)) - @test L2 === L - @test Q2 === Q - @test L * Q ≈ A - @test isunitary(Q) - - noL = similar(A, 0, n) - Q2 = similar(Q) - lq_full!(copy!(Ac, A), (noL, Q2)) - @test Q == Q2 - - # Transposed QR algorithm - qr_alg = LAPACK_HouseholderQR() - lq_alg = LQViaTransposedQR(qr_alg) - L2, Q2 = @constinferred lq_full!(copy!(Ac, A), (L, Q), lq_alg) - @test L2 === L - @test Q2 === Q - @test L * Q ≈ A - @test Q * Q' ≈ I - noL = similar(A, 0, n) - Q2 = similar(Q) - lq_full!(copy!(Ac, A), (noL, Q2), lq_alg) - @test Q == Q2 - - # unblocked algorithm - lq_full!(copy!(Ac, A), (L, Q); blocksize = 1) - @test L * Q ≈ A - @test isunitary(Q) - lq_full!(copy!(Ac, A), (noL, Q2); blocksize = 1) - @test Q == Q2 - if n == m - lq_full!(copy!(Q2, A), (noL, Q2); blocksize = 1) # in-place Q - @test Q ≈ Q2 - end - qr_alg = LAPACK_HouseholderQR(; blocksize = 1) - lq_alg = LQViaTransposedQR(qr_alg) - lq_full!(copy!(Ac, A), (L, Q), lq_alg) - @test L * Q ≈ A - @test Q * Q' ≈ I - lq_full!(copy!(Ac, A), (noL, Q2), lq_alg) - @test Q == Q2 - if n == m - lq_full!(copy!(Q2, A), (noL, Q2), lq_alg) # in-place Q - @test Q ≈ Q2 - end - - # other blocking - lq_full!(copy!(Ac, A), (L, Q); blocksize = 18) - @test L * Q ≈ A - @test isunitary(Q) - lq_full!(copy!(Ac, A), (noL, Q2); blocksize = 18) - @test Q == Q2 - # pivoted - @test_throws ArgumentError lq_full!(copy!(Ac, A), (L, Q); pivoted = true) - # positive - lq_full!(copy!(Ac, A), (L, Q); positive = true) - @test L * Q ≈ A - @test isunitary(Q) - @test all(>=(zero(real(T))), real(diag(L))) - lq_full!(copy!(Ac, A), (noL, Q2); positive = true) - @test Q == Q2 - - qr_alg = LAPACK_HouseholderQR(; positive = true) - lq_alg = LQViaTransposedQR(qr_alg) - lq_full!(copy!(Ac, A), (L, Q), lq_alg) - @test L * Q ≈ A - @test Q * Q' ≈ I - @test all(>=(zero(real(T))), real(diag(L))) - lq_full!(copy!(Ac, A), (noL, Q2), lq_alg) - @test Q == Q2 - - # positive and blocksize 1 - lq_full!(copy!(Ac, A), (L, Q); positive = true, blocksize = 1) - @test L * Q ≈ A - @test isunitary(Q) - @test all(>=(zero(real(T))), real(diag(L))) - lq_full!(copy!(Ac, A), (noL, Q2); positive = true, blocksize = 1) - @test Q == Q2 - end -end - -@testset "lq_compact, lq_full and lq_null for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) - rng = StableRNG(123) - atol = eps(real(T))^(3 / 4) - for m in (54, 0) - Ad = randn(rng, T, m) - A = Diagonal(Ad) - - # compact - L, Q = @constinferred lq_compact(A) - @test Q isa Diagonal{T} && size(Q) == (m, m) - @test L isa Diagonal{T} && size(L) == (m, m) - @test L * Q ≈ A - @test isunitary(Q) - - # compact and positive - Lp, Qp = @constinferred lq_compact(A; positive = true) - @test Qp isa Diagonal{T} && size(Qp) == (m, m) - @test Lp isa Diagonal{T} && size(Lp) == (m, m) - @test Lp * Qp ≈ A - @test isunitary(Qp) - @test all(≥(zero(real(T))), real(diag(Lp))) && - all(≈(zero(real(T)); atol), imag(diag(Lp))) - - # full - L, Q = @constinferred lq_full(A) - @test Q isa Diagonal{T} && size(Q) == (m, m) - @test L isa Diagonal{T} && size(L) == (m, m) - @test L * Q ≈ A - @test isunitary(Q) - - # full and positive - Lp, Qp = @constinferred lq_full(A; positive = true) - @test Qp isa Diagonal{T} && size(Qp) == (m, m) - @test Lp isa Diagonal{T} && size(Lp) == (m, m) - @test Lp * Qp ≈ A - @test isunitary(Qp) - @test all(≥(zero(real(T))), real(diag(Lp))) && - all(≈(zero(real(T)); atol), imag(diag(Lp))) - - # null - N = @constinferred lq_null(A) - @test N isa AbstractMatrix{T} && size(N) == (0, m) - end -end -=# diff --git a/test/qr.jl b/test/qr.jl index d7f974be..44af61f0 100644 --- a/test/qr.jl +++ b/test/qr.jl @@ -28,219 +28,3 @@ for T in (BLASFloats..., GenericFloats...) AT = Diagonal{T, Vector{T}} TestSuite.test_qr(AT, m; test_pivoted = false, test_blocksize = false) end - - -# @testset "qr_compact! and qr_null! for T = $T" for T in BLASFloats -# rng = StableRNG(123) -# m = 54 -# for n in (37, m, 63) -# minmn = min(m, n) -# A = randn(rng, T, m, n) -# Q, R = @constinferred qr_compact(A) -# @test Q isa Matrix{T} && size(Q) == (m, minmn) -# @test R isa Matrix{T} && size(R) == (minmn, n) -# @test Q * R ≈ A -# N = @constinferred qr_null(A) -# @test N isa Matrix{T} && size(N) == (m, m - minmn) -# @test isisometric(Q) -# @test maximum(abs, A' * N) < eps(real(T))^(2 / 3) -# @test isisometric(N) -# -# Ac = similar(A) -# Q2, R2 = @constinferred qr_compact!(copy!(Ac, A), (Q, R)) -# @test Q2 === Q -# @test R2 === R -# N2 = @constinferred qr_null!(copy!(Ac, A), N) -# @test N2 === N -# -# Q2 = similar(Q) -# noR = similar(A, minmn, 0) -# qr_compact!(copy!(Ac, A), (Q2, noR)) -# @test Q == Q2 -# -# # unblocked algorithm -# qr_compact!(copy!(Ac, A), (Q, R); blocksize = 1) -# @test Q * R ≈ A -# @test isisometric(Q) -# qr_compact!(copy!(Ac, A), (Q2, noR); blocksize = 1) -# @test Q == Q2 -# qr_compact!(copy!(Ac, A), (Q2, noR); blocksize = 1) -# qr_null!(copy!(Ac, A), N; blocksize = 1) -# @test maximum(abs, A' * N) < eps(real(T))^(2 / 3) -# @test isisometric(N) -# if n <= m -# qr_compact!(copy!(Q2, A), (Q2, noR); blocksize = 1) # in-place Q -# @test Q ≈ Q2 -# @test_throws ArgumentError qr_compact!(copy!(Q2, A), (Q2, R); blocksize = 1) -# @test_throws ArgumentError qr_compact!(copy!(Q2, A), (Q2, noR); positive = true) -# @test_throws ArgumentError qr_compact!(copy!(Q2, A), (Q2, noR); blocksize = 8) -# end -# # other blocking -# qr_compact!(copy!(Ac, A), (Q, R); blocksize = 8) -# @test Q * R ≈ A -# @test isisometric(Q) -# qr_compact!(copy!(Ac, A), (Q2, noR); blocksize = 8) -# @test Q == Q2 -# qr_null!(copy!(Ac, A), N; blocksize = 8) -# @test maximum(abs, A' * N) < eps(real(T))^(2 / 3) -# @test isisometric(N) -# -# # pivoted -# qr_compact!(copy!(Ac, A), (Q, R); pivoted = true) -# @test Q * R ≈ A -# @test Q' * Q ≈ I -# qr_compact!(copy!(Ac, A), (Q2, noR); pivoted = true) -# @test Q == Q2 -# # positive -# qr_compact!(copy!(Ac, A), (Q, R); positive = true) -# @test Q * R ≈ A -# @test isisometric(Q) -# @test all(>=(zero(real(T))), real(diag(R))) -# qr_compact!(copy!(Ac, A), (Q2, noR); positive = true) -# @test Q == Q2 -# # positive and blocksize 1 -# qr_compact!(copy!(Ac, A), (Q, R); positive = true, blocksize = 1) -# @test Q * R ≈ A -# @test isisometric(Q) -# @test all(>=(zero(real(T))), real(diag(R))) -# qr_compact!(copy!(Ac, A), (Q2, noR); positive = true, blocksize = 1) -# @test Q == Q2 -# # positive and pivoted -# qr_compact!(copy!(Ac, A), (Q, R); positive = true, pivoted = true) -# @test Q * R ≈ A -# @test isisometric(Q) -# if n <= m -# # the following test tries to find the diagonal element (in order to test positivity) -# # before the column permutation. This only works if all columns have a diagonal -# # element -# for j in 1:n -# i = findlast(!iszero, view(R, :, j)) -# @test real(R[i, j]) >= zero(real(T)) -# end -# end -# qr_compact!(copy!(Ac, A), (Q2, noR); positive = true, pivoted = true) -# @test Q == Q2 -# end -# end -# -# @testset "qr_full! for T = $T" for T in BLASFloats -# rng = StableRNG(123) -# m = 54 -# for n in (37, m, 63) -# minmn = min(m, n) -# A = randn(rng, T, m, n) -# Q, R = qr_full(A) -# @test Q isa Matrix{T} && size(Q) == (m, m) -# @test R isa Matrix{T} && size(R) == (m, n) -# @test Q * R ≈ A -# @test isunitary(Q) -# -# Ac = similar(A) -# Q2 = similar(Q) -# noR = similar(A, m, 0) -# Q2, R2 = @constinferred qr_full!(copy!(Ac, A), (Q, R)) -# @test Q2 === Q -# @test R2 === R -# @test Q * R ≈ A -# @test isunitary(Q) -# qr_full!(copy!(Ac, A), (Q2, noR)) -# @test Q == Q2 -# -# # unblocked algorithm -# qr_full!(copy!(Ac, A), (Q, R); blocksize = 1) -# @test Q * R ≈ A -# @test isunitary(Q) -# qr_full!(copy!(Ac, A), (Q2, noR); blocksize = 1) -# @test Q == Q2 -# if n == m -# qr_full!(copy!(Q2, A), (Q2, noR); blocksize = 1) # in-place Q -# @test Q ≈ Q2 -# end -# # other blocking -# qr_full!(copy!(Ac, A), (Q, R); blocksize = 8) -# @test Q * R ≈ A -# @test isunitary(Q) -# qr_full!(copy!(Ac, A), (Q2, noR); blocksize = 8) -# @test Q == Q2 -# # pivoted -# qr_full!(copy!(Ac, A), (Q, R); pivoted = true) -# @test Q * R ≈ A -# @test isunitary(Q) -# qr_full!(copy!(Ac, A), (Q2, noR); pivoted = true) -# @test Q == Q2 -# # positive -# qr_full!(copy!(Ac, A), (Q, R); positive = true) -# @test Q * R ≈ A -# @test isunitary(Q) -# @test all(>=(zero(real(T))), real(diag(R))) -# qr_full!(copy!(Ac, A), (Q2, noR); positive = true) -# @test Q == Q2 -# # positive and blocksize 1 -# qr_full!(copy!(Ac, A), (Q, R); positive = true, blocksize = 1) -# @test Q * R ≈ A -# @test isunitary(Q) -# @test all(>=(zero(real(T))), real(diag(R))) -# qr_full!(copy!(Ac, A), (Q2, noR); positive = true, blocksize = 1) -# @test Q == Q2 -# # positive and pivoted -# qr_full!(copy!(Ac, A), (Q, R); positive = true, pivoted = true) -# @test Q * R ≈ A -# @test isunitary(Q) -# if n <= m -# # the following test tries to find the diagonal element (in order to test positivity) -# # before the column permutation. This only works if all columns have a diagonal -# # element -# for j in 1:n -# i = findlast(!iszero, view(R, :, j)) -# @test real(R[i, j]) >= zero(real(T)) -# end -# end -# qr_full!(copy!(Ac, A), (Q2, noR); positive = true, pivoted = true) -# @test Q == Q2 -# end -# end -# -# @testset "qr_compact, qr_full and qr_null for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) -# rng = StableRNG(123) -# atol = eps(real(T))^(3 / 4) -# for m in (54, 0) -# Ad = randn(rng, T, m) -# A = Diagonal(Ad) -# -# # compact -# Q, R = @constinferred qr_compact(A) -# @test Q isa Diagonal{T} && size(Q) == (m, m) -# @test R isa Diagonal{T} && size(R) == (m, m) -# @test Q * R ≈ A -# @test isunitary(Q) -# -# # compact and positive -# Qp, Rp = @constinferred qr_compact(A; positive = true) -# @test Qp isa Diagonal{T} && size(Qp) == (m, m) -# @test Rp isa Diagonal{T} && size(Rp) == (m, m) -# @test Qp * Rp ≈ A -# @test isunitary(Qp) -# @test all(≥(zero(real(T))), real(diag(Rp))) && -# all(≈(zero(real(T)); atol), imag(diag(Rp))) -# -# # full -# Q, R = @constinferred qr_full(A) -# @test Q isa Diagonal{T} && size(Q) == (m, m) -# @test R isa Diagonal{T} && size(R) == (m, m) -# @test Q * R ≈ A -# @test isunitary(Q) -# -# # full and positive -# Qp, Rp = @constinferred qr_full(A; positive = true) -# @test Qp isa Diagonal{T} && size(Qp) == (m, m) -# @test Rp isa Diagonal{T} && size(Rp) == (m, m) -# @test Qp * Rp ≈ A -# @test isunitary(Qp) -# @test all(≥(zero(real(T))), real(diag(Rp))) && -# all(≈(zero(real(T)); atol), imag(diag(Rp))) -# -# # null -# N = @constinferred qr_null(A) -# @test N isa AbstractMatrix{T} && size(N) == (m, 0) -# end -# end diff --git a/test/runtests.jl b/test/runtests.jl index 2345fa80..590e04f5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -57,9 +57,6 @@ end using CUDA if CUDA.functional() - @safetestset "CUDA LQ" begin - include("cuda/lq.jl") - end @safetestset "CUDA Projections" begin include("cuda/projections.jl") end @@ -82,9 +79,6 @@ end using AMDGPU if AMDGPU.functional() - @safetestset "AMDGPU LQ" begin - include("amd/lq.jl") - end @safetestset "AMDGPU Projections" begin include("amd/projections.jl") end From c110afc9a2d3be6b32b321e7d00d0847ba6efb7a Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Wed, 19 Nov 2025 06:24:47 -0500 Subject: [PATCH 14/14] Try moving TestExtras around and adding projections --- test/lq.jl | 1 - test/polar.jl | 1 - test/projections.jl | 107 +++++--------------------- test/qr.jl | 1 - test/testsuite/TestSuite.jl | 3 +- test/testsuite/lq.jl | 2 + test/testsuite/polar.jl | 2 + test/testsuite/projections.jl | 136 ++++++++++++++++++++++++++++++++++ test/testsuite/qr.jl | 2 + 9 files changed, 164 insertions(+), 91 deletions(-) create mode 100644 test/testsuite/projections.jl diff --git a/test/lq.jl b/test/lq.jl index eb9fd773..8caa2ca6 100644 --- a/test/lq.jl +++ b/test/lq.jl @@ -1,6 +1,5 @@ using MatrixAlgebraKit using Test -using TestExtras using StableRNGs using LinearAlgebra: diag, I, Diagonal using MatrixAlgebraKit: LQViaTransposedQR, LAPACK_HouseholderQR diff --git a/test/polar.jl b/test/polar.jl index 6ac7b321..1731724a 100644 --- a/test/polar.jl +++ b/test/polar.jl @@ -1,6 +1,5 @@ using MatrixAlgebraKit using Test -using TestExtras using StableRNGs using LinearAlgebra: LinearAlgebra, I, isposdef diff --git a/test/projections.jl b/test/projections.jl index 56bc1a04..d26ae7de 100644 --- a/test/projections.jl +++ b/test/projections.jl @@ -4,93 +4,26 @@ using TestExtras using StableRNGs using LinearAlgebra: LinearAlgebra, Diagonal, norm, normalize! -const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) - -@testset "project_(anti)hermitian! for T = $T" for T in BLASFloats - rng = StableRNG(123) - m = 54 - noisefactor = eps(real(T))^(3 / 4) - for alg in (NativeBlocked(blocksize = 16), NativeBlocked(blocksize = 32), NativeBlocked(blocksize = 64)) - A = randn(rng, T, m, m) - Ah = (A + A') / 2 - Aa = (A - A') / 2 - Ac = copy(A) - - Bh = project_hermitian(A, alg) - @test ishermitian(Bh) - @test Bh ≈ Ah - @test A == Ac - Bh_approx = Bh + noisefactor * Aa - @test !ishermitian(Bh_approx) - @test ishermitian(Bh_approx; rtol = 10 * noisefactor) - - Ba = project_antihermitian(A, alg) - @test isantihermitian(Ba) - @test Ba ≈ Aa - @test A == Ac - Ba_approx = Ba + noisefactor * Ah - @test !isantihermitian(Ba_approx) - @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) - - Bh = project_hermitian!(Ac, alg) - @test Bh === Ac - @test ishermitian(Bh) - @test Bh ≈ Ah - - copy!(Ac, A) - Ba = project_antihermitian!(Ac, alg) - @test Ba === Ac - @test isantihermitian(Ba) - @test Ba ≈ Aa +BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) +GenericFloats = (Float16, BigFloat, Complex{BigFloat}) + +@isdefined(TestSuite) || include("testsuite/TestSuite.jl") +using .TestSuite + +m = 54 +for T in BLASFloats, n in (37, m, 63) + TestSuite.seed_rng!(123) + TestSuite.test_projections(T, (m, n)) + if CUDA.functional() + TestSuite.test_projections(CuMatrix{T}, (m, n); test_pivoted = false, test_blocksize = false) + TestSuite.test_projections(Diagonal{T, CuVector{T}}, m; test_pivoted = false, test_blocksize = false) end - - # test approximate error calculation - A = normalize!(randn(rng, T, m, m)) - Ah = project_hermitian(A) - Aa = project_antihermitian(A) - - Ah_approx = Ah + noisefactor * Aa - ϵ = norm(project_antihermitian(Ah_approx)) - @test !ishermitian(Ah_approx; atol = (999 // 1000) * ϵ) - @test ishermitian(Ah_approx; atol = (1001 // 1000) * ϵ) - - Aa_approx = Aa + noisefactor * Ah - ϵ = norm(project_hermitian(Aa_approx)) - @test !isantihermitian(Aa_approx; atol = (999 // 1000) * ϵ) - @test isantihermitian(Aa_approx; atol = (1001 // 1000) * ϵ) -end - -@testset "project_isometric! for T = $T" for T in BLASFloats - rng = StableRNG(123) - m = 54 - @testset "size ($m, $n)" for n in (37, m) - k = min(m, n) - if LinearAlgebra.LAPACK.version() < v"3.12.0" - svdalgs = (LAPACK_DivideAndConquer(), LAPACK_QRIteration(), LAPACK_Bisection()) - else - svdalgs = (LAPACK_DivideAndConquer(), LAPACK_QRIteration(), LAPACK_Bisection(), LAPACK_Jacobi()) - end - algs = (PolarViaSVD.(svdalgs)..., PolarNewton()) - @testset "algorithm $alg" for alg in algs - A = randn(rng, T, m, n) - W = project_isometric(A, alg) - @test isisometric(W) - W2 = project_isometric(W, alg) - @test W2 ≈ W # stability of the projection - @test W * (W' * A) ≈ A - - Ac = similar(A) - W2 = @constinferred project_isometric!(copy!(Ac, A), W, alg) - @test W2 === W - @test isisometric(W) - - # test that W is closer to A then any other isometry - for k in 1:10 - δA = randn(rng, T, m, n) - W = project_isometric(A, alg) - W2 = project_isometric(A + δA / 100, alg) - @test norm(A - W2) > norm(A - W) - end - end + if AMDGPU.functional() + TestSuite.test_projections(ROCMatrix{T}, (m, n); test_pivoted = false, test_blocksize = false) + TestSuite.test_projections(Diagonal{T, ROCVector{T}}, m; test_pivoted = false, test_blocksize = false) end end +for T in (BLASFloats..., GenericFloats...) + AT = Diagonal{T, Vector{T}} + TestSuite.test_projections(AT, m; test_pivoted = false, test_blocksize = false) +end diff --git a/test/qr.jl b/test/qr.jl index 44af61f0..3229df05 100644 --- a/test/qr.jl +++ b/test/qr.jl @@ -1,6 +1,5 @@ using MatrixAlgebraKit using Test -using TestExtras using StableRNGs using LinearAlgebra: diag, I, Diagonal using CUDA, AMDGPU diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 765d7d25..1563517c 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -8,7 +8,7 @@ Suite of tests that may be used for all packages inheriting from MatrixAlgebraKi """ module TestSuite -using Test, TestExtras +using Test using MatrixAlgebraKit using MatrixAlgebraKit: diagview using LinearAlgebra: Diagonal, norm, istriu, istril @@ -66,5 +66,6 @@ end include("qr.jl") include("lq.jl") include("polar.jl") +include("projections.jl") end diff --git a/test/testsuite/lq.jl b/test/testsuite/lq.jl index 5ad41c3b..575d299f 100644 --- a/test/testsuite/lq.jl +++ b/test/testsuite/lq.jl @@ -1,3 +1,5 @@ +using TestExtras + function test_lq(T::Type, sz; kwargs...) summary_str = testargs_summary(T, sz) return @testset "lq $summary_str" begin diff --git a/test/testsuite/polar.jl b/test/testsuite/polar.jl index 67eece8b..af747aa2 100644 --- a/test/testsuite/polar.jl +++ b/test/testsuite/polar.jl @@ -1,3 +1,5 @@ +using TestExtras + function test_polar(T::Type, sz; kwargs...) summary_str = testargs_summary(T, sz) return @testset "polar $summary_str" begin diff --git a/test/testsuite/projections.jl b/test/testsuite/projections.jl new file mode 100644 index 00000000..02cef217 --- /dev/null +++ b/test/testsuite/projections.jl @@ -0,0 +1,136 @@ +using TestExtras + +function test_projections(T::Type, sz; kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "projections $summary_str" begin + test_project_antihermitian(T, sz; kwargs...) + test_project_hermitian(T, sz; kwargs...) + test_project_isometric(T, sz; kwargs...) + end +end + +function test_project_antihermitian( + T::Type, sz; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "project_antihermitian! $summary_str" begin + noisefactor = eps(real(T))^(3 / 4) + algs = (NativeBlocked(blocksize = 16), NativeBlocked(blocksize = 32), NativeBlocked(blocksize = 64)) + @testset "algorithm $alg" for alg in algs + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + Ah = (A + A') / 2 + Aa = (A - A') / 2 + + Ba = project_antihermitian(A, alg) + @test isantihermitian(Ba) + @test Ba ≈ Aa + @test A == Ac + Ba_approx = Ba + noisefactor * Ah + @test !isantihermitian(Ba_approx) + @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) + + copy!(Ac, A) + Ba = project_antihermitian!(Ac, alg) + @test Ba === Ac + @test isantihermitian(Ba) + @test Ba ≈ Aa + end + + # test approximate error calculation + A = normalize!(randn(rng, T, m, m)) + Ah = project_hermitian(A) + Aa = project_antihermitian(A) + + Ah_approx = Ah + noisefactor * Aa + ϵ = norm(project_antihermitian(Ah_approx)) + @test !ishermitian(Ah_approx; atol = (999 // 1000) * ϵ) + @test ishermitian(Ah_approx; atol = (1001 // 1000) * ϵ) + + Aa_approx = Aa + noisefactor * Ah + ϵ = norm(project_hermitian(Aa_approx)) + @test !isantihermitian(Aa_approx; atol = (999 // 1000) * ϵ) + @test isantihermitian(Aa_approx; atol = (1001 // 1000) * ϵ) + end +end + +function test_project_hermitian( + T::Type, sz; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "project_hermitian! $summary_str" begin + noisefactor = eps(real(T))^(3 / 4) + algs = (NativeBlocked(blocksize = 16), NativeBlocked(blocksize = 32), NativeBlocked(blocksize = 64)) + @testset "algorithm $alg" for alg in algs + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + Ah = (A + A') / 2 + Aa = (A - A') / 2 + + Bh = project_hermitian(A, alg) + @test ishermitian(Bh) + @test Bh ≈ Ah + @test A == Ac + Bh_approx = Bh + noisefactor * Aa + @test !ishermitian(Bh_approx) + @test ishermitian(Bh_approx; rtol = 10 * noisefactor) + + Bh = project_hermitian!(Ac, alg) + @test Bh === Ac + @test ishermitian(Bh) + @test Bh ≈ Ah + end + + # test approximate error calculation + A = normalize!(randn(rng, T, m, m)) + Ah = project_hermitian(A) + Aa = project_antihermitian(A) + + Ah_approx = Ah + noisefactor * Aa + ϵ = norm(project_antihermitian(Ah_approx)) + @test !ishermitian(Ah_approx; atol = (999 // 1000) * ϵ) + @test ishermitian(Ah_approx; atol = (1001 // 1000) * ϵ) + + Aa_approx = Aa + noisefactor * Ah + ϵ = norm(project_hermitian(Aa_approx)) + @test !isantihermitian(Aa_approx; atol = (999 // 1000) * ϵ) + @test isantihermitian(Aa_approx; atol = (1001 // 1000) * ϵ) + end +end + +function test_project_isometric( + T::Type, sz; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "project_isometric! $summary_str" begin + algs = (PolarViaSVD(), PolarNewton()) + @testset "algorithm $alg" for alg in algs + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + k = min(size(A)...) + W = project_isometric(A, alg) + @test isisometric(W) + W2 = project_isometric(W, alg) + @test W2 ≈ W # stability of the projection + @test W * (W' * A) ≈ A + + W2 = @constinferred project_isometric!(Ac, W, alg) + @test W2 === W + @test isisometric(W) + + # test that W is closer to A then any other isometry + for k in 1:10 + δA = randn(rng, T, m, n) + W = project_isometric(A, alg) + W2 = project_isometric(A + δA / 100, alg) + @test norm(A - W2) > norm(A - W) + end + end + end +end diff --git a/test/testsuite/qr.jl b/test/testsuite/qr.jl index c354c2e4..09e2a637 100644 --- a/test/testsuite/qr.jl +++ b/test/testsuite/qr.jl @@ -1,3 +1,5 @@ +using TestExtras + function test_qr(T::Type, sz; kwargs...) summary_str = testargs_summary(T, sz) return @testset "qr $summary_str" begin