From 9825144347d6f86e86460b822cf8ad15968985e4 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Mon, 22 Dec 2025 13:28:15 +0100 Subject: [PATCH 1/4] Testsuite for polar --- src/implementations/polar.jl | 12 ++-- src/yalapack.jl | 16 +++--- test/amd/polar.jl | 83 ---------------------------- test/cuda/polar.jl | 83 ---------------------------- test/lq.jl | 3 +- test/polar.jl | 104 +++++++++++------------------------ test/qr.jl | 3 +- test/runtests.jl | 12 +--- test/testsuite/TestSuite.jl | 1 + test/testsuite/polar.jl | 85 ++++++++++++++++++++++++++++ 10 files changed, 140 insertions(+), 262 deletions(-) delete mode 100644 test/amd/polar.jl delete mode 100644 test/cuda/polar.jl create mode 100644 test/testsuite/polar.jl diff --git a/src/implementations/polar.jl b/src/implementations/polar.jl index 9a767441..74e86ae1 100644 --- a/src/implementations/polar.jl +++ b/src/implementations/polar.jl @@ -6,8 +6,8 @@ copy_input(::typeof(right_polar), A) = copy_input(svd_full, A) function check_input(::typeof(left_polar!), A::AbstractMatrix, WP, ::AbstractAlgorithm) m, n = size(A) W, P = WP - m >= n || - throw(ArgumentError("input matrix needs at least as many rows as columns")) + m ≥ n || + throw(ArgumentError("input matrix needs at least as many rows ($m) as columns ($n)")) @assert W isa AbstractMatrix && P isa AbstractMatrix @check_size(W, (m, n)) @check_scalar(W, A) @@ -18,8 +18,8 @@ end function check_input(::typeof(right_polar!), A::AbstractMatrix, PWᴴ, ::AbstractAlgorithm) m, n = size(A) P, Wᴴ = PWᴴ - n >= m || - throw(ArgumentError("input matrix needs at least as many columns as rows")) + n ≥ m || + throw(ArgumentError("input matrix needs at least as many columns ($n) as rows ($m)")) @assert P isa AbstractMatrix && Wᴴ isa AbstractMatrix isempty(P) || @check_size(P, (m, m)) @check_scalar(P, A) @@ -152,7 +152,7 @@ function _right_polarnewton!(A::AbstractMatrix, Wᴴ, P = similar(A, (0, 0)); to else # m == n L = A Lc = view(Wᴴ, 1:m, 1:m) - copy!(Lc, L) + Lc .= L Lᴴinv = ldiv!(lu!(Lc)', one!(Lᴴinv)) end γ = sqrt(norm(Lᴴinv) / norm(L)) # scaling factor @@ -168,7 +168,7 @@ function _right_polarnewton!(A::AbstractMatrix, Wᴴ, P = similar(A, (0, 0)); to rmul!(L, γ) rmul!(Lᴴinv, 1 / γ) L, Lᴴinv = _avgdiff!(L, Lᴴinv) - copy!(Lc, L) + Lc .= L conv = norm(Lᴴinv, Inf) i += 1 end diff --git a/src/yalapack.jl b/src/yalapack.jl index b2d4c27e..18541d41 100644 --- a/src/yalapack.jl +++ b/src/yalapack.jl @@ -2162,7 +2162,7 @@ for (gesvd, gesdd, gesvdx, gejsv, gesvj, elty, relty) in jobu = 'N' else size(U, 1) == m || - throw(DimensionMismatch("row size mismatch between A and U")) + throw(DimensionMismatch("row size mismatch between A ($m) and U ($(size(U, 1)))")) size(U, 2) >= (range == 'I' ? iu - il + 1 : minmn) || throw(DimensionMismatch("invalid column size of U")) jobu = 'V' @@ -2171,13 +2171,13 @@ for (gesvd, gesdd, gesvdx, gejsv, gesvj, elty, relty) in jobvt = 'N' else size(Vᴴ, 2) == n || - throw(DimensionMismatch("column size mismatch between A and Vᴴ")) + throw(DimensionMismatch("column size mismatch between A ($n) and Vᴴ ($(size(Vᴴ, 2)))")) size(Vᴴ, 1) >= (range == 'I' ? iu - il + 1 : minmn) || throw(DimensionMismatch("invalid row size of Vᴴ")) jobvt = 'V' end length(S) == minmn || - throw(DimensionMismatch("length mismatch between A and S")) + throw(DimensionMismatch("length mismatch between A ($minmn) and S ($(length(S)))")) lda = max(1, stride(A, 2)) ldu = max(1, stride(U, 2)) @@ -2247,15 +2247,15 @@ for (gesvd, gesdd, gesvdx, gejsv, gesvj, elty, relty) in require_one_based_indexing(A, U, Vᴴ, S) chkstride1(A, U, Vᴴ, S) m, n = size(A) - m >= n || - throw(ArgumentError("gejsv! requires a matrix with at least as many rows as columns")) + m ≥ n || + throw(ArgumentError("gejsv! requires a matrix with at least as many rows ($m) as columns ($n)")) joba = 'G' if length(U) == 0 jobu = 'N' else size(U, 1) == m || - throw(DimensionMismatch("row size mismatch between A and U")) + throw(DimensionMismatch("row size mismatch between A ($m) and U ($(size(U, 1)))")) if size(U, 2) == n jobu = 'U' elseif size(U, 2) == m @@ -2268,7 +2268,7 @@ for (gesvd, gesdd, gesvdx, gejsv, gesvj, elty, relty) in jobv = 'N' else size(Vᴴ, 2) == n || - throw(DimensionMismatch("column size mismatch between A and Vᴴ")) + throw(DimensionMismatch("column size mismatch between A ($n) and Vᴴ ($(size(Vᴴ, 2)))")) if size(Vᴴ, 1) == n jobv = 'V' else @@ -2276,7 +2276,7 @@ for (gesvd, gesdd, gesvdx, gejsv, gesvj, elty, relty) in end end length(S) == n || - throw(DimensionMismatch("length mismatch between A and S")) + throw(DimensionMismatch("length mismatch between A ($minmn) and S ($(length(S)))")) lda = max(1, stride(A, 2)) mv = Ref{BlasInt}() # unused diff --git a/test/amd/polar.jl b/test/amd/polar.jl deleted file mode 100644 index 4040b674..00000000 --- a/test/amd/polar.jl +++ /dev/null @@ -1,83 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using LinearAlgebra: LinearAlgebra, I, isposdef, Hermitian -using MatrixAlgebraKit: PolarViaSVD -using AMDGPU - -@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) - svd_algs = (ROCSOLVER_QRIteration(), ROCSOLVER_Jacobi()) - @testset "algorithm $svd_alg" for svd_alg in svd_algs - A = ROCArray(randn(rng, T, m, n)) - alg = PolarViaSVD(svd_alg) - W, P = left_polar(A; alg) - @test W isa ROCMatrix{T} && size(W) == (m, n) - @test P isa ROCMatrix{T} && size(P) == (n, n) - @test W * P ≈ A - @test isisometric(W) - # work around extremely strict Julia criteria for Hermiticity - @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) && isposdef(Hermitian(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) - # work around extremely strict Julia criteria for Hermiticity - @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) && isposdef(Hermitian(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(Hermitian(project_hermitian!(P))) - end - 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) - svd_algs = (ROCSOLVER_QRIteration(), ROCSOLVER_Jacobi()) - @testset "algorithm $svd_alg" for svd_alg in svd_algs - A = ROCArray(randn(rng, T, m, n)) - alg = PolarViaSVD(svd_alg) - P, Wᴴ = right_polar(A; alg) - @test Wᴴ isa ROCMatrix{T} && size(Wᴴ) == (m, n) - @test P isa ROCMatrix{T} && size(P) == (m, m) - @test P * Wᴴ ≈ A - @test isisometric(Wᴴ; side = :right) - # work around extremely strict Julia criteria for Hermiticity - @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) && isposdef(Hermitian(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) - # work around extremely strict Julia criteria for Hermiticity - @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) && isposdef(Hermitian(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(Hermitian(project_hermitian!(P))) - end - end -end diff --git a/test/cuda/polar.jl b/test/cuda/polar.jl deleted file mode 100644 index 1f512367..00000000 --- a/test/cuda/polar.jl +++ /dev/null @@ -1,83 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using LinearAlgebra: LinearAlgebra, I, isposdef, Hermitian -using MatrixAlgebraKit: PolarViaSVD -using CUDA - -@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) - svd_algs = (CUSOLVER_QRIteration(), CUSOLVER_Jacobi()) - @testset "algorithm $svd_alg" for svd_alg in svd_algs - A = CuArray(randn(rng, T, m, n)) - alg = PolarViaSVD(svd_alg) - W, P = left_polar(A; alg) - @test W isa CuMatrix{T} && size(W) == (m, n) - @test P isa CuMatrix{T} && size(P) == (n, n) - @test W * P ≈ A - @test isisometric(W) - # work around extremely strict Julia criteria for Hermiticity - @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) && isposdef(Hermitian(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) - # work around extremely strict Julia criteria for Hermiticity - @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) && isposdef(Hermitian(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(Hermitian(project_hermitian!(P))) - end - 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) - svd_algs = (CUSOLVER_QRIteration(), CUSOLVER_Jacobi()) - @testset "algorithm $svd_alg" for svd_alg in svd_algs - A = CuArray(randn(rng, T, m, n)) - alg = PolarViaSVD(svd_alg) - P, Wᴴ = right_polar(A; alg) - @test Wᴴ isa CuMatrix{T} && size(Wᴴ) == (m, n) - @test P isa CuMatrix{T} && size(P) == (m, m) - @test P * Wᴴ ≈ A - @test isisometric(Wᴴ; side = :right) - # work around extremely strict Julia criteria for Hermiticity - @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) && isposdef(Hermitian(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) - # work around extremely strict Julia criteria for Hermiticity - @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) && isposdef(Hermitian(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(Hermitian(project_hermitian!(P))) - end - end -end diff --git a/test/lq.jl b/test/lq.jl index f9eedd5f..2f34e846 100644 --- a/test/lq.jl +++ b/test/lq.jl @@ -35,7 +35,8 @@ for T in (BLASFloats..., GenericFloats...), n in (37, m, 63) TestSuite.test_lq_algs(Diagonal{T, ROCVector{T}}, m, (DiagonalAlgorithm(),)) end end - elseif !is_buildkite + end + if !is_buildkite if T ∈ BLASFloats TestSuite.test_lq(T, (m, n)) LAPACK_LQ_ALGS = ( diff --git a/test/polar.jl b/test/polar.jl index 8087149e..21d74cd6 100644 --- a/test/polar.jl +++ b/test/polar.jl @@ -1,83 +1,45 @@ using MatrixAlgebraKit using Test -using TestExtras using StableRNGs -using LinearAlgebra: LinearAlgebra, I, isposdef +using LinearAlgebra: Diagonal +using CUDA, AMDGPU -@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) +BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) +GenericFloats = (BigFloat, Complex{BigFloat}) - 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) +@isdefined(TestSuite) || include("testsuite/TestSuite.jl") +using .TestSuite - 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) +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" - 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)) +m = 54 +for T in (BLASFloats..., GenericFloats...), n in (37, m, 63) + TestSuite.seed_rng!(123) + if T ∈ BLASFloats + if CUDA.functional() + CUDA_POLAR_ALGS = (PolarViaSVD.((CUSOLVER_QRIteration(), CUSOLVER_SVDPolar(), CUSOLVER_Jacobi()))..., PolarNewton()) + TestSuite.test_polar(CuMatrix{T}, (m, n), CUDA_POLAR_ALGS) + n == m && TestSuite.test_polar(Diagonal{T, CuVector{T}}, m, (PolarNewton(),)) + end + if AMDGPU.functional() + ROC_POLAR_ALGS = (PolarViaSVD.((ROCSOLVER_QRIteration(), ROCSOLVER_Jacobi()))..., PolarNewton()) + TestSuite.test_polar(ROCMatrix{T}, (m, n), ROC_POLAR_ALGS) + n == m && TestSuite.test_polar(Diagonal{T, ROCVector{T}}, m, (PolarNewton(),)) end 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)) + if !is_buildkite + if T ∈ BLASFloats + LAPACK_POLAR_ALGS = (PolarViaSVD.((LAPACK_QRIteration(), LAPACK_Bisection(), LAPACK_DivideAndConquer()))..., PolarNewton()) + TestSuite.test_polar(T, (m, n), LAPACK_POLAR_ALGS) + LAPACK_JACOBI = (PolarViaSVD(LAPACK_Jacobi()),) + TestSuite.test_polar(T, (m, n), LAPACK_JACOBI; test_right = false) + elseif T ∈ GenericFloats + GLA_POLAR_ALGS = (PolarViaSVD.((GLA_QRIteration(),))..., PolarNewton()) + TestSuite.test_polar(T, (m, n), GLA_POLAR_ALGS) + end + if m == n + AT = Diagonal{T, Vector{T}} + TestSuite.test_polar(AT, m, (PolarNewton(),)) end end end diff --git a/test/qr.jl b/test/qr.jl index 19d1cc00..8b420fc6 100644 --- a/test/qr.jl +++ b/test/qr.jl @@ -34,7 +34,8 @@ for T in (BLASFloats..., GenericFloats...), n in (37, m, 63) TestSuite.test_qr_algs(Diagonal{T, ROCVector{T}}, m, (DiagonalAlgorithm(),)) end end - elseif !is_buildkite + end + if !is_buildkite if T ∈ BLASFloats TestSuite.test_qr(T, (m, n)) LAPACK_QR_ALGS = ( diff --git a/test/runtests.jl b/test/runtests.jl index 37796e2a..a32f26d6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -28,9 +28,6 @@ if !is_buildkite @safetestset "Schur Decomposition" begin include("schur.jl") end - @safetestset "Polar Decomposition" begin - include("polar.jl") - end @safetestset "Image and Null Space" begin include("orthnull.jl") end @@ -71,6 +68,9 @@ end include("qr.jl") include("lq.jl") end +@safetestset "Polar Decomposition" begin + include("polar.jl") +end using CUDA if CUDA.functional() @@ -86,9 +86,6 @@ if CUDA.functional() @safetestset "CUDA Hermitian Eigenvalue Decomposition" begin include("cuda/eigh.jl") end - @safetestset "CUDA Polar Decomposition" begin - include("cuda/polar.jl") - end @safetestset "CUDA Image and Null Space" begin include("cuda/orthnull.jl") end @@ -105,9 +102,6 @@ if AMDGPU.functional() @safetestset "AMDGPU Hermitian Eigenvalue Decomposition" begin include("amd/eigh.jl") end - @safetestset "AMDGPU Polar Decomposition" begin - include("amd/polar.jl") - end @safetestset "AMDGPU Image and Null Space" begin include("amd/orthnull.jl") end diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 83119b71..737c8c68 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -71,5 +71,6 @@ is_pivoted(alg::MatrixAlgebraKit.LQViaTransposedQR) = is_pivoted(alg.qr_alg) 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..c610ba34 --- /dev/null +++ b/test/testsuite/polar.jl @@ -0,0 +1,85 @@ +using TestExtras +using LinearAlgebra: isposdef + +function test_polar(T::Type, sz, algs; test_right = true, kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "polar $summary_str" begin + (length(sz) == 1 || sz[1] ≥ sz[2]) && test_left_polar(T, sz, algs; kwargs...) + # doesn't work for Jacobi + (length(sz) == 1 || sz[2] ≥ sz[1]) && test_right && test_right_polar(T, sz, algs; kwargs...) + end +end + +function test_left_polar( + T::Type, sz, algs; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "left_polar! algorithm $alg $summary_str" for alg in algs + @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 ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) + @test isposdef(project_hermitian!(P)) + + W2, P2 = @testinferred left_polar!(Ac, (W, P), alg) + @test W2 === W + @test P2 === P + @test W * P ≈ A + @test isisometric(W) + @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) + @test isposdef(project_hermitian!(P)) + + noP = similar(P, (0, 0)) + W2, P2 = @testinferred 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, algs; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "right_polar! algorithm $alg $summary_str" 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 ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) + @test isposdef(project_hermitian!(P)) + + P2, Wᴴ2 = @testinferred right_polar!(Ac, (P, Wᴴ), alg) + @test P2 === P + @test Wᴴ2 === Wᴴ + @test P * Wᴴ ≈ A + @test isisometric(Wᴴ; side = :right) + @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P)) + @test isposdef(project_hermitian!(P)) + + noP = similar(P, (0, 0)) + P2, Wᴴ2 = @testinferred 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 From f91f42aba6da3c3c0c63f28614822f38d2ff4bea Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Mon, 22 Dec 2025 18:46:51 +0100 Subject: [PATCH 2/4] Make left polar newton more GPU friendly --- src/implementations/polar.jl | 8 ++++---- test/polar.jl | 10 ++++++---- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/src/implementations/polar.jl b/src/implementations/polar.jl index 74e86ae1..b94b2afa 100644 --- a/src/implementations/polar.jl +++ b/src/implementations/polar.jl @@ -107,19 +107,19 @@ function _left_polarnewton!(A::AbstractMatrix, W, P = similar(A, (0, 0)); tol = if m > n # initial QR Q, R = qr_compact!(A) Rc = view(A, 1:n, 1:n) - copy!(Rc, R) + Rc .= R Rᴴinv = ldiv!(UpperTriangular(Rc)', one!(Rᴴinv)) else # m == n R = A Rc = view(W, 1:n, 1:n) - copy!(Rc, R) + Rc .= R Rᴴinv = ldiv!(lu!(Rc)', one!(Rᴴinv)) end γ = sqrt(norm(Rᴴinv) / norm(R)) # scaling factor rmul!(R, γ) rmul!(Rᴴinv, 1 / γ) R, Rᴴinv = _avgdiff!(R, Rᴴinv) - copy!(Rc, R) + Rc .= R i = 1 conv = norm(Rᴴinv, Inf) while i < maxiter && conv > tol @@ -128,7 +128,7 @@ function _left_polarnewton!(A::AbstractMatrix, W, P = similar(A, (0, 0)); tol = rmul!(R, γ) rmul!(Rᴴinv, 1 / γ) R, Rᴴinv = _avgdiff!(R, Rᴴinv) - copy!(Rc, R) + Rc .= R conv = norm(Rᴴinv, Inf) i += 1 end diff --git a/test/polar.jl b/test/polar.jl index 21d74cd6..7bea5760 100644 --- a/test/polar.jl +++ b/test/polar.jl @@ -17,14 +17,16 @@ for T in (BLASFloats..., GenericFloats...), n in (37, m, 63) TestSuite.seed_rng!(123) if T ∈ BLASFloats if CUDA.functional() - CUDA_POLAR_ALGS = (PolarViaSVD.((CUSOLVER_QRIteration(), CUSOLVER_SVDPolar(), CUSOLVER_Jacobi()))..., PolarNewton()) + # PolarNewton does not work yet on GPU + CUDA_POLAR_ALGS = (PolarViaSVD.((CUSOLVER_QRIteration(), CUSOLVER_SVDPolar(), CUSOLVER_Jacobi()))...,) # PolarNewton()) TestSuite.test_polar(CuMatrix{T}, (m, n), CUDA_POLAR_ALGS) - n == m && TestSuite.test_polar(Diagonal{T, CuVector{T}}, m, (PolarNewton(),)) + #n == m && TestSuite.test_polar(Diagonal{T, CuVector{T}}, m, (PolarNewton(),)) end if AMDGPU.functional() - ROC_POLAR_ALGS = (PolarViaSVD.((ROCSOLVER_QRIteration(), ROCSOLVER_Jacobi()))..., PolarNewton()) + # PolarNewton does not work yet on GPU + ROC_POLAR_ALGS = (PolarViaSVD.((ROCSOLVER_QRIteration(), ROCSOLVER_Jacobi()))...,) # PolarNewton()) TestSuite.test_polar(ROCMatrix{T}, (m, n), ROC_POLAR_ALGS) - n == m && TestSuite.test_polar(Diagonal{T, ROCVector{T}}, m, (PolarNewton(),)) + #n == m && TestSuite.test_polar(Diagonal{T, ROCVector{T}}, m, (PolarNewton(),)) end end if !is_buildkite From 00c578cfbc946e0749070229d0079c64ece61ae0 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Mon, 22 Dec 2025 20:20:34 +0100 Subject: [PATCH 3/4] Don't test Jacobi on old Julia --- test/polar.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/polar.jl b/test/polar.jl index 7bea5760..3b5d9bbc 100644 --- a/test/polar.jl +++ b/test/polar.jl @@ -33,8 +33,10 @@ for T in (BLASFloats..., GenericFloats...), n in (37, m, 63) if T ∈ BLASFloats LAPACK_POLAR_ALGS = (PolarViaSVD.((LAPACK_QRIteration(), LAPACK_Bisection(), LAPACK_DivideAndConquer()))..., PolarNewton()) TestSuite.test_polar(T, (m, n), LAPACK_POLAR_ALGS) - LAPACK_JACOBI = (PolarViaSVD(LAPACK_Jacobi()),) - TestSuite.test_polar(T, (m, n), LAPACK_JACOBI; test_right = false) + if LinearAlgebra.LAPACK.version() ≥ v"3.12.0" + LAPACK_JACOBI = (PolarViaSVD(LAPACK_Jacobi()),) + TestSuite.test_polar(T, (m, n), LAPACK_JACOBI; test_right = false) + end elseif T ∈ GenericFloats GLA_POLAR_ALGS = (PolarViaSVD.((GLA_QRIteration(),))..., PolarNewton()) TestSuite.test_polar(T, (m, n), GLA_POLAR_ALGS) From d468a6ef656c0b1fa5e7b084948e297cfb879035 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Mon, 22 Dec 2025 15:08:37 -0500 Subject: [PATCH 4/4] Include properly --- test/polar.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/polar.jl b/test/polar.jl index 3b5d9bbc..4f1df5a9 100644 --- a/test/polar.jl +++ b/test/polar.jl @@ -1,7 +1,7 @@ using MatrixAlgebraKit using Test using StableRNGs -using LinearAlgebra: Diagonal +using LinearAlgebra: Diagonal, LAPACK using CUDA, AMDGPU BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) @@ -33,7 +33,7 @@ for T in (BLASFloats..., GenericFloats...), n in (37, m, 63) if T ∈ BLASFloats LAPACK_POLAR_ALGS = (PolarViaSVD.((LAPACK_QRIteration(), LAPACK_Bisection(), LAPACK_DivideAndConquer()))..., PolarNewton()) TestSuite.test_polar(T, (m, n), LAPACK_POLAR_ALGS) - if LinearAlgebra.LAPACK.version() ≥ v"3.12.0" + if LAPACK.version() ≥ v"3.12.0" LAPACK_JACOBI = (PolarViaSVD(LAPACK_Jacobi()),) TestSuite.test_polar(T, (m, n), LAPACK_JACOBI; test_right = false) end