diff --git a/Project.toml b/Project.toml index f09ad557..0e9a26cc 100644 --- a/Project.toml +++ b/Project.toml @@ -33,10 +33,11 @@ GenericSchur = "0.5.6" JET = "0.9, 0.10" LinearAlgebra = "1" Mooncake = "0.4.183" +Random = "1" SafeTestsets = "0.1" StableRNGs = "1" Test = "1" -TestExtras = "0.2,0.3" +TestExtras = "0.3.2" Zygote = "0.7" julia = "1.10" @@ -47,6 +48,7 @@ ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -54,4 +56,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", "Mooncake"] +test = ["Aqua", "JET", "SafeTestsets", "Test", "TestExtras", "ChainRulesCore", "ChainRulesTestUtils", "Random", "StableRNGs", "Zygote", "CUDA", "AMDGPU", "GenericLinearAlgebra", "GenericSchur", "Mooncake"] diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index 9b558b11..e4637875 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -56,7 +56,7 @@ function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix, D, ::GLA_QRIteration) return eigvals!(Hermitian(A); sortby = real) end -function MatrixAlgebraKit.default_qr_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} +function MatrixAlgebraKit.default_qr_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{Float16, ComplexF16, BigFloat, Complex{BigFloat}}}} return GLA_HouseholderQR(; kwargs...) end @@ -109,7 +109,7 @@ function _gla_householder_qr!(A::AbstractMatrix, Q, R; positive = false, blocksi return Q, R end -function MatrixAlgebraKit.default_lq_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} +function MatrixAlgebraKit.default_lq_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{Float16, ComplexF16, BigFloat, Complex{BigFloat}}}} return MatrixAlgebraKit.LQViaTransposedQR(GLA_HouseholderQR(; kwargs...)) end diff --git a/src/implementations/lq.jl b/src/implementations/lq.jl index a1648d3e..758bf53e 100644 --- a/src/implementations/lq.jl +++ b/src/implementations/lq.jl @@ -153,8 +153,8 @@ function _lapack_lq!( computeL = length(L) > 0 inplaceQ = Q === A - if pivoted - throw(ArgumentError("LAPACK does not provide an implementation for a pivoted LQ decomposition")) + if pivoted && (blocksize > 1) + throw(ArgumentError("LAPACK does not provide a blocked implementation for a pivoted LQ decomposition")) end if inplaceQ && (computeL || positive || blocksize > 1 || n < m) throw(ArgumentError("inplace Q only supported if matrix is wide (`m <= n`), L is not required, and using the unblocked algorithm (`blocksize=1`) with `positive=false`")) diff --git a/src/implementations/qr.jl b/src/implementations/qr.jl index a3c9d1f5..54951a45 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; 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) 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/lq.jl b/test/amd/lq.jl deleted file mode 100644 index ef46bbe9..00000000 --- a/test/amd/lq.jl +++ /dev/null @@ -1,163 +0,0 @@ -using MatrixAlgebraKit -using MatrixAlgebraKit: diagview -using Test -using TestExtras -using StableRNGs -using AMDGPU -using LinearAlgebra - -include(joinpath("..", "utilities.jl")) - -BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) - -@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 = 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 BLASFloats - 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 - -@testset "lq_compact, lq_full and lq_null for Diagonal{$T}" for T in BLASFloats - rng = StableRNG(123) - atol = eps(real(T))^(3 / 4) - for m in (54, 0) - Ad = ROCArray(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(isposdef.(diagview(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(isposdef.(diagview(Lp))) - - # null - N = @constinferred lq_null(A) - @test N isa AbstractMatrix{T} && size(N) == (0, m) - end -end diff --git a/test/amd/qr.jl b/test/amd/qr.jl deleted file mode 100644 index b0708ae1..00000000 --- a/test/amd/qr.jl +++ /dev/null @@ -1,167 +0,0 @@ -using MatrixAlgebraKit -using MatrixAlgebraKit: diagview -using Test -using TestExtras -using StableRNGs -using AMDGPU -using LinearAlgebra - -include(joinpath("..", "utilities.jl")) - -BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) - -@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 = 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 BLASFloats - 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 - -@testset "qr_compact, qr_full and qr_null for Diagonal{$T}" for T in BLASFloats - rng = StableRNG(123) - atol = eps(real(T))^(3 / 4) - for m in (54, 0) - Ad = ROCArray(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(isposdef.(diagview(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(isposdef.(diagview(Rp))) - - # null - N = @constinferred qr_null(A) - @test N isa AbstractMatrix{T} && size(N) == (m, 0) - end -end diff --git a/test/cuda/lq.jl b/test/cuda/lq.jl deleted file mode 100644 index de904c88..00000000 --- a/test/cuda/lq.jl +++ /dev/null @@ -1,163 +0,0 @@ -using MatrixAlgebraKit -using MatrixAlgebraKit: diagview -using Test -using TestExtras -using StableRNGs -using CUDA -using LinearAlgebra - -include(joinpath("..", "utilities.jl")) - -BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) - -@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 = 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 BLASFloats - 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 - -@testset "lq_compact, lq_full and lq_null for Diagonal{$T}" for T in BLASFloats - rng = StableRNG(123) - atol = eps(real(T))^(3 / 4) - for m in (54, 0) - Ad = CuArray(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(isposdef.(diagview(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(isposdef.(diagview(Lp))) - - # null - N = @constinferred lq_null(A) - @test N isa AbstractMatrix{T} && size(N) == (0, m) - end -end diff --git a/test/cuda/qr.jl b/test/cuda/qr.jl deleted file mode 100644 index 73fa3985..00000000 --- a/test/cuda/qr.jl +++ /dev/null @@ -1,168 +0,0 @@ -using MatrixAlgebraKit -using MatrixAlgebraKit: diagview -using Test -using TestExtras -using StableRNGs -using CUDA -using LinearAlgebra -using LinearAlgebra: isposdef - -include(joinpath("..", "utilities.jl")) - -BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) - -@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 = 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 BLASFloats - 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 - -@testset "qr_compact, qr_full and qr_null for Diagonal{$T}" for T in BLASFloats - rng = StableRNG(123) - atol = eps(real(T))^(3 / 4) - for m in (54, 0) - Ad = CuArray(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(isposdef.(diagview(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(isposdef.(diagview(Rp))) - - # null - N = @constinferred qr_null(A) - @test N isa AbstractMatrix{T} && size(N) == (m, 0) - end -end diff --git a/test/genericlinearalgebra/lq.jl b/test/genericlinearalgebra/lq.jl deleted file mode 100644 index dc186dcb..00000000 --- a/test/genericlinearalgebra/lq.jl +++ /dev/null @@ -1,124 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using LinearAlgebra: diag, I, Diagonal -using GenericLinearAlgebra - -eltypes = (BigFloat, Complex{BigFloat}) - -@testset "qr_compact! for T = $T" for T in eltypes - - rng = StableRNG(123) - m = 54 - for n in (37, m, 63) - minmn = min(m, n) - m = 54 - 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) - - Ac = similar(A) - L2, Q2 = @constinferred lq_compact!(copy!(Ac, A), (L, Q)) - @test L2 === L - @test Q2 === Q - - noL = similar(A, 0, minmn) - Q2 = similar(Q) - lq_compact!(copy!(Ac, A), (noL, Q2)) - @test Q == Q2 - - # Transposed QR algorithm - qr_alg = GLA_HouseholderQR() - lq_alg = LQViaTransposedQR(qr_alg) - L2, Q2 = @constinferred lq_compact!(copy!(Ac, A), (L, Q), lq_alg) - @test L2 === L - @test Q2 === Q - noL = similar(A, 0, minmn) - Q2 = similar(Q) - lq_compact!(copy!(Ac, A), (noL, Q2), lq_alg) - @test Q == Q2 - - @test_throws ArgumentError lq_compact(A; blocksize = 2) - @test_throws ArgumentError lq_compact(A; 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 - end -end - -@testset "lq_full! for T = $T" for T in eltypes - 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[1:minmn, n] ≈ Q2[1:minmn, n] - - # Transposed QR algorithm - qr_alg = GLA_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[1:minmn, n] ≈ Q2[1:minmn, n] - - # Argument errors for unsupported options - @test_throws ArgumentError lq_full(A; blocksize = 2) - @test_throws ArgumentError lq_full(A; 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[1:minmn, n] ≈ Q2[1:minmn, n] - - qr_alg = GLA_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[1:minmn, n] ≈ Q2[1:minmn, n] - - # 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[1:minmn, n] ≈ Q2[1:minmn, n] - end -end diff --git a/test/genericlinearalgebra/qr.jl b/test/genericlinearalgebra/qr.jl deleted file mode 100644 index 3ce530bb..00000000 --- a/test/genericlinearalgebra/qr.jl +++ /dev/null @@ -1,109 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using LinearAlgebra: diag, I, Diagonal -using GenericLinearAlgebra - -eltypes = (BigFloat, Complex{BigFloat}) - -@testset "qr_compact! for T = $T" for T in eltypes - - rng = StableRNG(123) - m = 54 - for n in (37, m, 63) - minmn = min(m, n) - m = 54 - 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 - - Ac = similar(A) - Q2, R2 = @constinferred qr_compact!(copy!(Ac, A), (Q, R)) - @test Q2 === Q - @test R2 === R - - Q2 = similar(Q) - noR = similar(A, minmn, 0) - qr_compact!(copy!(Ac, A), (Q2, noR)) - @test Q == Q2 - - @test_throws ArgumentError qr_compact(A; blocksize = 2) - @test_throws ArgumentError qr_compact(A; pivoted = true) - - # 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 - end -end - -@testset "qr_full! for T = $T" for T in eltypes - 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) - Qc, Rc = qr_compact(A) - @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 - - # Argument errors for unsupported options - @test_throws ArgumentError qr_full(A; blocksize = 2) - @test_throws ArgumentError qr_compact(A; pivoted = true) - - # 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 - 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 - end -end diff --git a/test/lq.jl b/test/lq.jl index 8de5a582..6829e935 100644 --- a/test/lq.jl +++ b/test/lq.jl @@ -1,255 +1,63 @@ using MatrixAlgebraKit using Test -using TestExtras using StableRNGs using LinearAlgebra: diag, I, Diagonal -using MatrixAlgebraKit: LQViaTransposedQR, LAPACK_HouseholderQR +using MatrixAlgebraKit: LQViaTransposedQR, LAPACK_HouseholderLQ +using CUDA, AMDGPU, GenericLinearAlgebra BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) GenericFloats = (Float16, BigFloat, Complex{BigFloat}) -@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) +@isdefined(TestSuite) || include("testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +m = 54 +for T in (BLASFloats..., GenericFloats...), n in (37, m, 63) + TestSuite.seed_rng!(123) + if T ∈ BLASFloats + if CUDA.functional() + CUDA_LQ_ALGS = LQViaTransposedQR.(CUSOLVER_HouseholderLQ(; positive = false), CUSOLVER_HouseholderLQ(; positive = true)) + TestSuite.test_lq(CuMatrix{T}, (m, n); test_pivoted = false, test_blocksize = false) + TestSuite.test_lq_algs(CuMatrix{T}, (m, n), CUDA_LQ_ALGS) + if n == m + TestSuite.test_lq(Diagonal{T, CuVector{T}}, m; test_pivoted = false, test_blocksize = false) + TestSuite.test_lq_algs(Diagonal{T, CuVector{T}}, m, (DiagonalAlgorithm(),)) + end 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 + if AMDGPU.functional() + ROC_LQ_ALGS = LQViaTransposedQR.(ROCSOLVER_HouseholderLQ(; positive = false), ROCSOLVER_HouseholderLQ(; positive = true)) + TestSuite.test_lq(ROCMatrix{T}, (m, n); test_pivoted = false, test_blocksize = false) + TestSuite.test_lq_algs(ROCMatrix{T}, (m, n), CUDA_LQ_ALGS) + if n == m + TestSuite.test_lq(Diagonal{T, ROCVector{T}}, m; test_pivoted = false, test_blocksize = false) + TestSuite.test_lq_algs(Diagonal{T, ROCVector{T}}, m, (DiagonalAlgorithm(),)) + end 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 + elseif !is_buildkite + if T ∈ BLASFloats + TestSuite.test_lq(T, (m, n)) + LAPACK_LQ_ALGS = ( + LAPACK_HouseholderLQ(; positive = false, pivoted = false, blocksize = 1), + LAPACK_HouseholderLQ(; positive = false, pivoted = false, blocksize = 8), + LAPACK_HouseholderLQ(; positive = false, pivoted = true, blocksize = 1), + #LAPACK_HouseholderLQ(; positive=false, pivoted=true, blocksize=8), # not supported + LAPACK_HouseholderLQ(; positive = true, pivoted = false, blocksize = 1), + LAPACK_HouseholderLQ(; positive = true, pivoted = false, blocksize = 8), + LAPACK_HouseholderLQ(; positive = true, pivoted = true, blocksize = 1), + #LAPACK_HouseholderLQ(; positive=true, pivoted=true, blocksize=8), # not supported + ) + TestSuite.test_lq_algs(T, (m, n), LAPACK_LQ_ALGS) + elseif T ∈ GenericFloats + TestSuite.test_lq(T, (m, n); test_null = false, test_pivoted = false, test_blocksize = false) + GLA_LQ_ALGS = (LQViaTransposedQR(GLA_HouseholderQR()),) + TestSuite.test_lq_algs(T, (m, n), GLA_LQ_ALGS; test_null = false) + end + if m == n + AT = Diagonal{T, Vector{T}} + TestSuite.test_lq(AT, m; test_pivoted = false, test_blocksize = false) + TestSuite.test_lq_algs(AT, m, (DiagonalAlgorithm(),)) 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 c4f0c9d6..19d1cc00 100644 --- a/test/qr.jl +++ b/test/qr.jl @@ -1,223 +1,62 @@ using MatrixAlgebraKit using Test -using TestExtras using StableRNGs using LinearAlgebra: diag, I, Diagonal +using CUDA, AMDGPU, GenericLinearAlgebra 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) - - 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) +@isdefined(TestSuite) || include("testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +m = 54 +for T in (BLASFloats..., GenericFloats...), n in (37, m, 63) + TestSuite.seed_rng!(123) + if T ∈ BLASFloats + if CUDA.functional() + CUDA_QR_ALGS = (CUSOLVER_HouseholderQR(; positive = false), CUSOLVER_HouseholderQR(; positive = true)) + TestSuite.test_qr(CuMatrix{T}, (m, n); test_pivoted = false, test_blocksize = false) + TestSuite.test_qr_algs(CuMatrix{T}, (m, n), CUDA_QR_ALGS) + if n == m + TestSuite.test_qr(Diagonal{T, CuVector{T}}, m; test_pivoted = false, test_blocksize = false) + TestSuite.test_qr_algs(Diagonal{T, CuVector{T}}, m, (DiagonalAlgorithm(),)) + end 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)) + if AMDGPU.functional() + ROC_QR_ALGS = (ROCSOLVER_HouseholderQR(; positive = false), ROCSOLVER_HouseholderQR(; positive = true)) + TestSuite.test_qr(ROCMatrix{T}, (m, n); test_pivoted = false, test_blocksize = false) + TestSuite.test_qr_algs(ROCMatrix{T}, (m, n), ROC_QR_ALGS) + if n == m + TestSuite.test_qr(Diagonal{T, ROCVector{T}}, m; test_pivoted = false, test_blocksize = false) + TestSuite.test_qr_algs(Diagonal{T, ROCVector{T}}, m, (DiagonalAlgorithm(),)) 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 + elseif !is_buildkite + if T ∈ BLASFloats + TestSuite.test_qr(T, (m, n)) + LAPACK_QR_ALGS = ( + LAPACK_HouseholderQR(; positive = false, pivoted = false, blocksize = 1), + LAPACK_HouseholderQR(; positive = false, pivoted = false, blocksize = 8), + LAPACK_HouseholderQR(; positive = false, pivoted = true, blocksize = 1), + #LAPACK_HouseholderQR(; positive=false, pivoted=true, blocksize=8), # not supported + LAPACK_HouseholderQR(; positive = true, pivoted = false, blocksize = 1), + LAPACK_HouseholderQR(; positive = true, pivoted = false, blocksize = 8), + LAPACK_HouseholderQR(; positive = true, pivoted = true, blocksize = 1), + #LAPACK_HouseholderQR(; positive=true, pivoted=true, blocksize=8), # not supported + ) + TestSuite.test_qr_algs(T, (m, n), LAPACK_QR_ALGS) + elseif T ∈ GenericFloats + TestSuite.test_qr(T, (m, n); test_null = false, test_pivoted = false, test_blocksize = false) + GLA_QR_ALGS = (GLA_HouseholderQR(),) + TestSuite.test_qr_algs(T, (m, n), GLA_QR_ALGS; test_null = false) 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 + if m == n + AT = Diagonal{T, Vector{T}} + TestSuite.test_qr(AT, m; test_pivoted = false, test_blocksize = false) + TestSuite.test_qr_algs(AT, m, (DiagonalAlgorithm(),)) 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 1ed1f456..392bd490 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -58,10 +58,6 @@ if !is_buildkite end using GenericLinearAlgebra - @safetestset "QR / LQ Decomposition" begin - include("genericlinearalgebra/qr.jl") - include("genericlinearalgebra/lq.jl") - end @safetestset "Singular Value Decomposition" begin include("genericlinearalgebra/svd.jl") end @@ -77,12 +73,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 @safetestset "CUDA Projections" begin include("cuda/projections.jl") end @@ -105,12 +95,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 @safetestset "AMDGPU Projections" begin include("amd/projections.jl") end diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl new file mode 100644 index 00000000..83119b71 --- /dev/null +++ b/test/testsuite/TestSuite.jl @@ -0,0 +1,75 @@ +# Based on the design of GPUArrays.jl + +""" + TestSuite + +Suite of tests that may be used for all packages inheriting from MatrixAlgebraKit. + +""" +module TestSuite + +using Test +using MatrixAlgebraKit +using MatrixAlgebraKit: diagview +using LinearAlgebra: Diagonal, norm, istriu, istril +using Random, StableRNGs +using AMDGPU, CUDA + +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 + +testargs_summary(args...) = string(args) + +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) +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)) + +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; 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)) + +is_positive(::MatrixAlgebraKit.AbstractAlgorithm) = false +is_pivoted(::MatrixAlgebraKit.AbstractAlgorithm) = false +is_positive(alg::MatrixAlgebraKit.LAPACK_HouseholderQR) = alg.positive +is_pivoted(alg::MatrixAlgebraKit.LAPACK_HouseholderQR) = alg.pivoted +is_positive(alg::MatrixAlgebraKit.LAPACK_HouseholderLQ) = alg.positive +is_pivoted(alg::MatrixAlgebraKit.LAPACK_HouseholderLQ) = alg.pivoted +is_positive(alg::MatrixAlgebraKit.CUSOLVER_HouseholderQR) = alg.positive +is_positive(alg::MatrixAlgebraKit.ROCSOLVER_HouseholderQR) = alg.positive +is_positive(alg::MatrixAlgebraKit.LQViaTransposedQR) = is_positive(alg.qr_alg) +is_pivoted(alg::MatrixAlgebraKit.LQViaTransposedQR) = is_pivoted(alg.qr_alg) + +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..d53ff443 --- /dev/null +++ b/test/testsuite/lq.jl @@ -0,0 +1,268 @@ +using TestExtras + +function test_lq(T::Type, sz; test_null = true, 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_null && test_lq_null(T, sz; kwargs...) + end +end + +function test_lq_algs(T::Type, sz, algs; test_null = true, kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "lq algorithms $summary_str" begin + test_lq_compact_algs(T, sz, algs; kwargs...) + test_lq_full_algs(T, sz, algs; kwargs...) + test_null && test_lq_null_algs(T, sz, algs; 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(Qpiv; side = :right, atol, rtol) + # pivoted AND blocksize not yet supported + @test_throws ArgumentError lq_compact(A; pivoted = true, blocksize = 2) + 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_compact_algs( + T::Type, sz, algs; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "lq_compact! algorithm $alg $summary_str" for alg in algs + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + + # does the elementary functionality work + L, Q = @testinferred lq_compact(A; alg) + @test L * Q ≈ A + @test isisometric(Q; side = :right, atol, rtol) + @test A == Ac + if !is_pivoted(alg) + @test istril(L) + if is_positive(alg) + @test has_positive_diagonal(L) + end + end + + # can I pass in outputs? + L2, Q2 = @testinferred lq_compact!(Ac, (L, Q); alg) + @test L2 * Q2 ≈ A + @test isisometric(Q2; side = :right, atol, rtol) + if !is_pivoted(alg) + @test istril(L2) + if is_positive(alg) + @test has_positive_diagonal(L2) + end + 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) + # pivoted AND blocksize not yet supported + @test_throws ArgumentError lq_full(A; pivoted = true, blocksize = 2) + 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_full_algs( + T::Type, sz, algs; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "lq_full! algorithm $alg $summary_str" for alg in algs + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + + # does the elementary functionality work + L, Q = @testinferred lq_full(A; alg) + @test L * Q ≈ A + @test isisometric(Q; side = :right, atol, rtol) + @test A == Ac + if !is_pivoted(alg) + @test istril(L) + if is_positive(alg) + @test has_positive_diagonal(L) + end + end + + # can I pass in outputs? + L2, Q2 = @testinferred lq_full!(Ac, (L, Q); alg) + @test L2 * Q2 ≈ A + @test isisometric(Q2; side = :right, atol, rtol) + if !is_pivoted(alg) + @test istril(L2) + if is_positive(alg) + @test has_positive_diagonal(L2) + end + 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!(Ac, 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 + @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 + +function test_lq_null_algs( + T::Type, sz, algs; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "lq_null! algorithm $alg $summary_str" for alg in algs + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + + # does the elementary functionality work + Nᴴ = @testinferred lq_null(A; alg) + @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!(Ac, Nᴴ; alg) + @test isrightnull(Nᴴ2, A; atol, rtol) + @test isisometric(Nᴴ2; side = :right, atol, rtol) + end +end diff --git a/test/testsuite/qr.jl b/test/testsuite/qr.jl new file mode 100644 index 00000000..8bba9e62 --- /dev/null +++ b/test/testsuite/qr.jl @@ -0,0 +1,271 @@ +using TestExtras + +function test_qr(T::Type, sz; test_null = true, 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_null && test_qr_null(T, sz; kwargs...) + end +end + +function test_qr_algs(T::Type, sz, algs; test_null = true, kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "qr algorithms $summary_str" begin + test_qr_compact_algs(T, sz, algs; kwargs...) + test_qr_full_algs(T, sz, algs; kwargs...) + test_null && test_qr_null_algs(T, sz, algs; kwargs...) + end +end + +# test correctness and interface for QR regardless of algorithm +function test_qr_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 "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; 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; 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; atol, rtol) + @test istriu(Rpos) + @test has_positive_diagonal(Rpos) + else + @test_throws Exception 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; atol, rtol) + # pivoted AND blocksize not yet supported + @test_throws ArgumentError qr_compact(A; pivoted = true, blocksize = 2) + else + @test_throws Exception 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; atol, rtol) + else + @test_throws Exception qr_compact(A; blocksize = 2) + end + end +end + +# test various algorithms +function test_qr_compact_algs( + T::Type, sz, algs; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "qr_compact! algorithm $alg $summary_str" for alg in algs + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + + # does the elementary functionality work + Q, R = @testinferred qr_compact(A; alg) + @test Q * R ≈ A + @test isisometric(Q; atol, rtol) + @test A == Ac + if !is_pivoted(alg) + @test istriu(R) + if is_positive(alg) + @test has_positive_diagonal(R) + end + end + + # can I pass in outputs? + Q2, R2 = @testinferred qr_compact!(Ac, (Q, R); alg) + @test Q2 * R2 ≈ A + @test isisometric(Q2; atol, rtol) + if !is_pivoted(alg) + @test istriu(R2) + if is_positive(alg) + @test has_positive_diagonal(R2) + end + end + end +end + +function test_qr_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 "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; 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; 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; atol, rtol) + @test istriu(Rpos) + @test has_positive_diagonal(Rpos) + else + @test_throws Exception 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; atol, rtol) + # pivoted AND blocksize not yet supported + @test_throws ArgumentError qr_full(A; pivoted = true, blocksize = 2) + else + @test_throws Exception 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; atol, rtol) + else + @test_throws Exception qr_full(A; blocksize = 2) + end + end +end + +function test_qr_full_algs( + T::Type, sz, algs; + test_blocksize = true, + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "qr_full! algorithm $alg $summary_str" for alg in algs + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + + # does the elementary functionality work + Q, R = @testinferred qr_full(A; alg) + @test Q * R ≈ A + @test isunitary(Q; atol, rtol) + @test A == Ac + if !is_pivoted(alg) + @test istriu(R) + if is_positive(alg) + @test has_positive_diagonal(R) + end + end + + # can I pass in outputs? + Q2, R2 = @testinferred qr_full!(Ac, (Q, R); alg) + @test Q2 * R2 ≈ A + @test isunitary(Q2; atol, rtol) + if !is_pivoted(alg) + @test istriu(R2) + if is_positive(alg) + @test has_positive_diagonal(R2) + end + end + end +end + +function test_qr_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 "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; 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; 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; atol, rtol) + @test isisometric(Npiv; atol, rtol) + else + @test_throws Exception 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; atol, rtol) + @test isisometric(Nblocked; atol, rtol) + else + @test_throws Exception qr_null(A; blocksize = 2) + end + end +end + +function test_qr_null_algs( + T::Type, sz, algs; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "qr_null! algorithm $alg $summary_str" for alg in algs + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + + # does the elementary functionality work + N = @testinferred qr_null(A; alg) + @test isleftnull(N, A; atol, rtol) + @test isisometric(N; atol, rtol) + @test A == Ac + + # can I pass in outputs? + N2 = @testinferred qr_null!(Ac, N; alg) + @test isleftnull(N2, A; atol, rtol) + @test isisometric(N2; atol, rtol) + end +end