From f8c97f144efc4a70ceaf1f89c76e2633ca7a328a Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Mon, 17 Nov 2025 13:29:43 -0500 Subject: [PATCH 01/14] A few more updates for GPU compatibility for TensorKit --- Project.toml | 3 + .../MatrixAlgebraKitAMDGPUExt.jl | 14 ++- .../MatrixAlgebraKitCUDAExt.jl | 14 ++- src/MatrixAlgebraKit.jl | 2 +- src/implementations/projections.jl | 15 ++- src/implementations/svd.jl | 10 ++ test/amd/projections.jl | 91 ++++++++++--------- test/cuda/projections.jl | 91 ++++++++++--------- 8 files changed, 143 insertions(+), 97 deletions(-) diff --git a/Project.toml b/Project.toml index a7851a45..ec1efe36 100644 --- a/Project.toml +++ b/Project.toml @@ -55,3 +55,6 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] test = ["Aqua", "JET", "SafeTestsets", "Test", "TestExtras", "ChainRulesCore", "ChainRulesTestUtils", "StableRNGs", "Zygote", "CUDA", "AMDGPU", "GenericLinearAlgebra", "GenericSchur", "Mooncake"] + +[sources] +CUDA = {url="https://github.com/juliagpu/cuda.jl", rev="master"} diff --git a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl index 04bee40e..a7830914 100644 --- a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl +++ b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl @@ -128,17 +128,23 @@ function MatrixAlgebraKit._project_hermitian_diag!(A::StridedROCMatrix, B::Strid end MatrixAlgebraKit.ishermitian_exact(A::StridedROCMatrix) = all(A .== adjoint(A)) -MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T} = - all(A.diag .== adjoint(A.diag)) +MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T <: Real} = true +function MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T <: Complex} + return all(isreal.(A.diag)) +end MatrixAlgebraKit.ishermitian_approx(A::StridedROCMatrix; kwargs...) = @invoke MatrixAlgebraKit.ishermitian_approx(A::Any; kwargs...) MatrixAlgebraKit.isantihermitian_exact(A::StridedROCMatrix) = all(A .== -adjoint(A)) -MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T} = - all(A.diag .== -adjoint(A.diag)) +MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T <: Real} = all(iszero(A.diag)) +MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T <: Complex} = + all(iszero.(real.(A.diag))) MatrixAlgebraKit.isantihermitian_approx(A::StridedROCMatrix; kwargs...) = @invoke MatrixAlgebraKit.isantihermitian_approx(A::Any; kwargs...) +function MatrixAlgebraKit.isantihermitian_approx(A::Diagonal{T, <:StridedROCVector{T}}; atol, rtol, kwargs...) where {T <: Real} + return sum(abs2, A.diag) ≤ max(atol, rtol * norm(A)) +end function MatrixAlgebraKit._avgdiff!(A::StridedROCMatrix, B::StridedROCMatrix) axes(A) == axes(B) || throw(DimensionMismatch()) diff --git a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl index e4c77453..c681d3eb 100644 --- a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl +++ b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl @@ -153,17 +153,23 @@ end MatrixAlgebraKit.ishermitian_exact(A::StridedCuMatrix) = all(A .== adjoint(A)) -MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T} = - all(A.diag .== adjoint(A.diag)) +MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T <: Real} = true +function MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T <: Complex} + return all(isreal.(A.diag)) +end MatrixAlgebraKit.ishermitian_approx(A::StridedCuMatrix; kwargs...) = @invoke MatrixAlgebraKit.ishermitian_approx(A::Any; kwargs...) MatrixAlgebraKit.isantihermitian_exact(A::StridedCuMatrix) = all(A .== -adjoint(A)) -MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T} = - all(A.diag .== -adjoint(A.diag)) +MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T <: Real} = all(iszero(A.diag)) +MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T <: Complex} = + all(iszero.(real.(A.diag))) MatrixAlgebraKit.isantihermitian_approx(A::StridedCuMatrix; kwargs...) = @invoke MatrixAlgebraKit.isantihermitian_approx(A::Any; kwargs...) +function MatrixAlgebraKit.isantihermitian_approx(A::Diagonal{T, <:StridedCuVector{T}}; atol, rtol, kwargs...) where {T <: Real} + return sum(abs2, A.diag) ≤ max(atol, rtol * norm(A)) +end function MatrixAlgebraKit._avgdiff!(A::StridedCuMatrix, B::StridedCuMatrix) axes(A) == axes(B) || throw(DimensionMismatch()) diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl index d3a04741..493f3a91 100644 --- a/src/MatrixAlgebraKit.jl +++ b/src/MatrixAlgebraKit.jl @@ -3,7 +3,7 @@ module MatrixAlgebraKit using LinearAlgebra: LinearAlgebra using LinearAlgebra: norm # TODO: eleminate if we use VectorInterface.jl? using LinearAlgebra: mul!, rmul!, lmul!, adjoint!, rdiv!, ldiv! -using LinearAlgebra: sylvester, lu! +using LinearAlgebra: sylvester, lu!, diagm using LinearAlgebra: isposdef, issymmetric using LinearAlgebra: Diagonal, diag, diagind, isdiag using LinearAlgebra: UpperTriangular, LowerTriangular diff --git a/src/implementations/projections.jl b/src/implementations/projections.jl index 8baf09a5..054b5b77 100644 --- a/src/implementations/projections.jl +++ b/src/implementations/projections.jl @@ -9,13 +9,15 @@ copy_input(::typeof(project_isometric), A) = copy_input(left_polar, A) function check_input(::typeof(project_hermitian!), A::AbstractMatrix, B::AbstractMatrix, ::AbstractAlgorithm) LinearAlgebra.checksquare(A) - n = Base.require_one_based_indexing(A) + Base.require_one_based_indexing(A) + n = size(A, 1) B === A || @check_size(B, (n, n)) return nothing end function check_input(::typeof(project_antihermitian!), A::AbstractMatrix, B::AbstractMatrix, ::AbstractAlgorithm) LinearAlgebra.checksquare(A) - n = Base.require_one_based_indexing(A) + Base.require_one_based_indexing(A) + n = size(A, 1) B === A || @check_size(B, (n, n)) return nothing end @@ -61,6 +63,15 @@ function project_isometric!(A::AbstractMatrix, W, alg::AbstractAlgorithm) return W end +function project_hermitian_native!(A::Diagonal, B::Diagonal, ::Val{anti}; kwargs...) where {anti} + if anti + diagview(A) .= imag.(diagview(B)) .* im + else + diagview(A) .= real.(diagview(B)) + end + return A +end + function project_hermitian_native!(A::AbstractMatrix, B::AbstractMatrix, anti::Val; blocksize = 32) n = size(A, 1) for j in 1:blocksize:n diff --git a/src/implementations/svd.jl b/src/implementations/svd.jl index c8b27f3f..ef0c2bc3 100644 --- a/src/implementations/svd.jl +++ b/src/implementations/svd.jl @@ -356,6 +356,7 @@ function svd_full!(A::AbstractMatrix, USVᴴ, alg::GPU_SVDAlgorithm) return USVᴴ end +svd_full!(A::Diagonal, USVᴴ, alg::GPU_SVDAlgorithm) = svd_full!(diagm(A.diag), USVᴴ, alg) function svd_trunc!(A::AbstractMatrix, USVᴴ, alg::TruncatedAlgorithm{<:GPU_Randomized}) check_input(svd_trunc!, A, USVᴴ, alg.alg) @@ -373,6 +374,7 @@ function svd_trunc!(A::AbstractMatrix, USVᴴ, alg::TruncatedAlgorithm{<:GPU_Ran return Utr, Str, Vᴴtr, ϵ end +svd_trunc!(A::Diagonal, USVᴴ, alg::GPU_SVDAlgorithm) = svd_trunc!(diagm(A.diag), USVᴴ, alg) function svd_compact!(A::AbstractMatrix, USVᴴ, alg::GPU_SVDAlgorithm) check_input(svd_compact!, A, USVᴴ, alg) @@ -396,6 +398,7 @@ function svd_compact!(A::AbstractMatrix, USVᴴ, alg::GPU_SVDAlgorithm) return USVᴴ end +svd_compact!(A::Diagonal, USVᴴ, alg::GPU_SVDAlgorithm) = svd_compact!(diagm(A.diag), USVᴴ, alg) _argmaxabs(x) = reduce(_largest, x; init = zero(eltype(x))) _largest(x, y) = abs(x) < abs(y) ? y : x @@ -418,3 +421,10 @@ function svd_vals!(A::AbstractMatrix, S, alg::GPU_SVDAlgorithm) return S end +function svd_vals!(A::Diagonal, S, alg::GPU_SVDAlgorithm) + check_input(svd_vals!, A, S, alg) + Ad = diagview(A) + S .= abs.(Ad) + sort!(S; rev = true) + return S +end diff --git a/test/amd/projections.jl b/test/amd/projections.jl index b22148e2..cab09206 100644 --- a/test/amd/projections.jl +++ b/test/amd/projections.jl @@ -12,37 +12,41 @@ const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) m = 54 noisefactor = eps(real(T))^(3 / 4) for alg in (NativeBlocked(blocksize = 16), NativeBlocked(blocksize = 32), NativeBlocked(blocksize = 64)) - A = ROCArray(randn(rng, T, m, m)) - Ah = (A + A') / 2 - Aa = (A - A') / 2 - Ac = copy(A) + for A in (ROCArray(randn(rng, T, m, m)), Diagonal(ROCArray(randn(rng, T, m)))) + Ah = (A + A') / 2 + Aa = (A - A') / 2 + Ac = copy(A) - Bh = project_hermitian(A, alg) - @test ishermitian(Bh) - @test Bh ≈ Ah - @test A == Ac - Bh_approx = Bh + noisefactor * Aa - @test !ishermitian(Bh_approx) - @test ishermitian(Bh_approx; rtol = 10 * noisefactor) + Bh = project_hermitian(A, alg) + @test ishermitian(Bh) + @test Bh ≈ Ah + @test A == Ac + # this is still hermitian for real Diagonals! + Bh_approx = Bh + noisefactor * Aa + if !isa(A, Diagonal) && !(T <: Real) + @test !ishermitian(Bh_approx) + end + @test ishermitian(Bh_approx; rtol = 10 * noisefactor) - Ba = project_antihermitian(A, alg) - @test isantihermitian(Ba) - @test Ba ≈ Aa - @test A == Ac - Ba_approx = Ba + noisefactor * Ah - @test !isantihermitian(Ba_approx) - @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) + Ba = project_antihermitian(A, alg) + @test isantihermitian(Ba) + @test Ba ≈ Aa + @test A == Ac + Ba_approx = Ba + noisefactor * Ah + @test !isantihermitian(Ba_approx) + @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) - Bh = project_hermitian!(Ac, alg) - @test Bh === Ac - @test ishermitian(Bh) - @test Bh ≈ Ah + Bh = project_hermitian!(Ac, alg) + @test Bh === Ac + @test ishermitian(Bh) + @test Bh ≈ Ah - copy!(Ac, A) - Ba = project_antihermitian!(Ac, alg) - @test Ba === Ac - @test isantihermitian(Ba) - @test Ba ≈ Aa + copy!(Ac, A) + Ba = project_antihermitian!(Ac, alg) + @test Ba === Ac + @test isantihermitian(Ba) + @test Ba ≈ Aa + end end end @@ -54,24 +58,25 @@ end svdalgs = (ROCSOLVER_QRIteration(), ROCSOLVER_Jacobi()) algs = (PolarViaSVD.(svdalgs)...,) # PolarNewton()) # TODO @testset "algorithm $alg" for alg in algs - A = ROCArray(randn(rng, T, m, n)) - W = project_isometric(A, alg) - @test isisometric(W) - W2 = project_isometric(W, alg) - @test W2 ≈ W # stability of the projection - @test W * (W' * A) ≈ A + for A in (ROCArray(randn(rng, T, m, n)), Diagonal(ROCArray(randn(rng, T, m)))) + W = project_isometric(A, alg) + @test isisometric(W) + W2 = project_isometric(W, alg) + @test W2 ≈ W # stability of the projection + @test W * (W' * A) ≈ A - Ac = similar(A) - W2 = @constinferred project_isometric!(copy!(Ac, A), W, alg) - @test W2 === W - @test isisometric(W) + Ac = similar(A) + W2 = @constinferred project_isometric!(copy!(Ac, A), W, alg) + @test W2 === W + @test isisometric(W) - # test that W is closer to A then any other isometry - for k in 1:10 - δA = ROCArray(randn(rng, T, m, n)) - W = project_isometric(A, alg) - W2 = project_isometric(A + δA / 100, alg) - @test norm(A - W2) > norm(A - W) + # test that W is closer to A then any other isometry + for k in 1:10 + δA = ROCArray(randn(rng, T, size(A)...)) + W = project_isometric(A, alg) + W2 = project_isometric(A + δA / 100, alg) + @test norm(A - W2) > norm(A - W) + end end end end diff --git a/test/cuda/projections.jl b/test/cuda/projections.jl index d9f7bd7f..9eb3d478 100644 --- a/test/cuda/projections.jl +++ b/test/cuda/projections.jl @@ -12,37 +12,41 @@ const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) m = 54 noisefactor = eps(real(T))^(3 / 4) for alg in (NativeBlocked(blocksize = 16), NativeBlocked(blocksize = 32), NativeBlocked(blocksize = 64)) - A = CuArray(randn(rng, T, m, m)) - Ah = (A + A') / 2 - Aa = (A - A') / 2 - Ac = copy(A) + for A in (CuArray(randn(rng, T, m, m)), Diagonal(CuArray(randn(rng, T, m)))) + Ah = (A + A') / 2 + Aa = (A - A') / 2 + Ac = copy(A) - Bh = project_hermitian(A, alg) - @test ishermitian(Bh) - @test Bh ≈ Ah - @test A == Ac - Bh_approx = Bh + noisefactor * Aa - @test !ishermitian(Bh_approx) - @test ishermitian(Bh_approx; rtol = 10 * noisefactor) + Bh = project_hermitian(A, alg) + @test ishermitian(Bh) + @test Bh ≈ Ah + @test A == Ac + # this is still hermitian for real Diagonals! + Bh_approx = Bh + noisefactor * Aa + if !isa(A, Diagonal) && !(T <: Real) + @test !ishermitian(Bh_approx) + end + @test ishermitian(Bh_approx; rtol = 10 * noisefactor) - Ba = project_antihermitian(A, alg) - @test isantihermitian(Ba) - @test Ba ≈ Aa - @test A == Ac - Ba_approx = Ba + noisefactor * Ah - @test !isantihermitian(Ba_approx) - @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) + Ba = project_antihermitian(A, alg) + @test isantihermitian(Ba) + @test Ba ≈ Aa + @test A == Ac + Ba_approx = Ba + noisefactor * Ah + @test !isantihermitian(Ba_approx) + @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) - Bh = project_hermitian!(Ac, alg) - @test Bh === Ac - @test ishermitian(Bh) - @test Bh ≈ Ah + Bh = project_hermitian!(Ac, alg) + @test Bh === Ac + @test ishermitian(Bh) + @test Bh ≈ Ah - copy!(Ac, A) - Ba = project_antihermitian!(Ac, alg) - @test Ba === Ac - @test isantihermitian(Ba) - @test Ba ≈ Aa + copy!(Ac, A) + Ba = project_antihermitian!(Ac, alg) + @test Ba === Ac + @test isantihermitian(Ba) + @test Ba ≈ Aa + end end end @@ -54,24 +58,25 @@ end svdalgs = (CUSOLVER_SVDPolar(), CUSOLVER_QRIteration(), CUSOLVER_Jacobi()) algs = (PolarViaSVD.(svdalgs)...,) # PolarNewton()) # TODO @testset "algorithm $alg" for alg in algs - A = CuArray(randn(rng, T, m, n)) - W = project_isometric(A, alg) - @test isisometric(W) - W2 = project_isometric(W, alg) - @test W2 ≈ W # stability of the projection - @test W * (W' * A) ≈ A + for A in (CuArray(randn(rng, T, m, n)), Diagonal(CuArray(randn(rng, T, m)))) + W = project_isometric(A, alg) + @test isisometric(W) + W2 = project_isometric(W, alg) + @test W2 ≈ W # stability of the projection + @test W * (W' * A) ≈ A - Ac = similar(A) - W2 = @constinferred project_isometric!(copy!(Ac, A), W, alg) - @test W2 === W - @test isisometric(W) + Ac = similar(A) + W2 = @constinferred project_isometric!(copy!(Ac, A), W, alg) + @test W2 === W + @test isisometric(W) - # test that W is closer to A then any other isometry - for k in 1:10 - δA = CuArray(randn(rng, T, m, n)) - W = project_isometric(A, alg) - W2 = project_isometric(A + δA / 100, alg) - @test norm(A - W2) > norm(A - W) + # test that W is closer to A then any other isometry + for k in 1:10 + δA = CuArray(randn(rng, T, size(A)...)) + W = project_isometric(A, alg) + W2 = project_isometric(A + δA / 100, alg) + @test norm(A - W2) > norm(A - W) + end end end end From 4ba611a6acfc2fd607c751a3ffd926055e9daefa Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Wed, 26 Nov 2025 08:44:50 +0100 Subject: [PATCH 02/14] Update src/implementations/projections.jl Co-authored-by: Jutho --- src/implementations/projections.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/implementations/projections.jl b/src/implementations/projections.jl index 054b5b77..857366ed 100644 --- a/src/implementations/projections.jl +++ b/src/implementations/projections.jl @@ -65,7 +65,7 @@ end function project_hermitian_native!(A::Diagonal, B::Diagonal, ::Val{anti}; kwargs...) where {anti} if anti - diagview(A) .= imag.(diagview(B)) .* im + diagview(A) .= _imimag.(diagview(B)) else diagview(A) .= real.(diagview(B)) end From 4af44a3e85d51c623251b5dfef7e8959288436c2 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Wed, 26 Nov 2025 03:16:02 -0500 Subject: [PATCH 03/14] Update antihermiticity check --- ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl | 2 +- ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl index a7830914..59ad1dd9 100644 --- a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl +++ b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl @@ -143,7 +143,7 @@ MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) wh MatrixAlgebraKit.isantihermitian_approx(A::StridedROCMatrix; kwargs...) = @invoke MatrixAlgebraKit.isantihermitian_approx(A::Any; kwargs...) function MatrixAlgebraKit.isantihermitian_approx(A::Diagonal{T, <:StridedROCVector{T}}; atol, rtol, kwargs...) where {T <: Real} - return sum(abs2, A.diag) ≤ max(atol, rtol * norm(A)) + return norm(A) ≤ max(atol, rtol) end function MatrixAlgebraKit._avgdiff!(A::StridedROCMatrix, B::StridedROCMatrix) diff --git a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl index c681d3eb..d8e54ef8 100644 --- a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl +++ b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl @@ -168,7 +168,7 @@ MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) whe MatrixAlgebraKit.isantihermitian_approx(A::StridedCuMatrix; kwargs...) = @invoke MatrixAlgebraKit.isantihermitian_approx(A::Any; kwargs...) function MatrixAlgebraKit.isantihermitian_approx(A::Diagonal{T, <:StridedCuVector{T}}; atol, rtol, kwargs...) where {T <: Real} - return sum(abs2, A.diag) ≤ max(atol, rtol * norm(A)) + return norm(A) ≤ max(atol, rtol) end function MatrixAlgebraKit._avgdiff!(A::StridedCuMatrix, B::StridedCuMatrix) From 8c574ec0ae24a028ceaa2466c4649b22c28db3fe Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Mon, 1 Dec 2025 06:34:31 -0500 Subject: [PATCH 04/14] Respond to comments, update tests --- Project.toml | 3 --- ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl | 2 +- ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl | 2 +- test/amd/projections.jl | 5 ++++- test/cuda/projections.jl | 5 ++++- 5 files changed, 10 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index ec1efe36..a7851a45 100644 --- a/Project.toml +++ b/Project.toml @@ -55,6 +55,3 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] test = ["Aqua", "JET", "SafeTestsets", "Test", "TestExtras", "ChainRulesCore", "ChainRulesTestUtils", "StableRNGs", "Zygote", "CUDA", "AMDGPU", "GenericLinearAlgebra", "GenericSchur", "Mooncake"] - -[sources] -CUDA = {url="https://github.com/juliagpu/cuda.jl", rev="master"} diff --git a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl index 59ad1dd9..bd9d31f0 100644 --- a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl +++ b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl @@ -143,7 +143,7 @@ MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) wh MatrixAlgebraKit.isantihermitian_approx(A::StridedROCMatrix; kwargs...) = @invoke MatrixAlgebraKit.isantihermitian_approx(A::Any; kwargs...) function MatrixAlgebraKit.isantihermitian_approx(A::Diagonal{T, <:StridedROCVector{T}}; atol, rtol, kwargs...) where {T <: Real} - return norm(A) ≤ max(atol, rtol) + return norm(A) ≤ atol end function MatrixAlgebraKit._avgdiff!(A::StridedROCMatrix, B::StridedROCMatrix) diff --git a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl index d8e54ef8..d2c6b3f3 100644 --- a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl +++ b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl @@ -168,7 +168,7 @@ MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) whe MatrixAlgebraKit.isantihermitian_approx(A::StridedCuMatrix; kwargs...) = @invoke MatrixAlgebraKit.isantihermitian_approx(A::Any; kwargs...) function MatrixAlgebraKit.isantihermitian_approx(A::Diagonal{T, <:StridedCuVector{T}}; atol, rtol, kwargs...) where {T <: Real} - return norm(A) ≤ max(atol, rtol) + return norm(A) ≤ atol end function MatrixAlgebraKit._avgdiff!(A::StridedCuMatrix, B::StridedCuMatrix) diff --git a/test/amd/projections.jl b/test/amd/projections.jl index cab09206..2e449be3 100644 --- a/test/amd/projections.jl +++ b/test/amd/projections.jl @@ -32,9 +32,12 @@ const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) @test isantihermitian(Ba) @test Ba ≈ Aa @test A == Ac + # this is still hermitian for real Diagonals! Ba_approx = Ba + noisefactor * Ah @test !isantihermitian(Ba_approx) - @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) + if !isa(A, Diagonal) && !(T <: Real) + @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) + end Bh = project_hermitian!(Ac, alg) @test Bh === Ac diff --git a/test/cuda/projections.jl b/test/cuda/projections.jl index 9eb3d478..b7fc42b1 100644 --- a/test/cuda/projections.jl +++ b/test/cuda/projections.jl @@ -34,7 +34,10 @@ const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) @test A == Ac Ba_approx = Ba + noisefactor * Ah @test !isantihermitian(Ba_approx) - @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) + # this is still hermitian for real Diagonals! + if !isa(A, Diagonal) && !(T <: Real) + @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) + end Bh = project_hermitian!(Ac, alg) @test Bh === Ac From 1b4399f52a8e132a7aab2a685e6726734d5a8298 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Mon, 1 Dec 2025 11:02:27 -0500 Subject: [PATCH 05/14] Add tests for SVD --- test/cuda/svd.jl | 176 +++++++++++++++++++++++++---------------------- 1 file changed, 92 insertions(+), 84 deletions(-) diff --git a/test/cuda/svd.jl b/test/cuda/svd.jl index e2173d03..6a2f6e1e 100644 --- a/test/cuda/svd.jl +++ b/test/cuda/svd.jl @@ -15,33 +15,34 @@ include(joinpath("..", "utilities.jl")) k = min(m, n) algs = (CUSOLVER_QRIteration(), CUSOLVER_SVDPolar(), CUSOLVER_Jacobi()) @testset "algorithm $alg" for alg in algs - minmn = min(m, n) - A = CuArray(randn(rng, T, m, n)) + As = m == n ? (CuArray(randn(rng, T, m, n)), Diagonal(CuArray(randn(rng, T, m)))) : (CuArray(randn(rng, T, m, n)),) + for A in As + minmn = min(m, n) + U, S, Vᴴ = svd_compact(A; alg) + @test U isa CuMatrix{T} && size(U) == (m, minmn) + @test S isa Diagonal{real(T), <:CuVector} && size(S) == (minmn, minmn) + @test Vᴴ isa CuMatrix{T} && size(Vᴴ) == (minmn, n) + @test U * S * Vᴴ ≈ A + @test isapproxone(U' * U) + @test isapproxone(Vᴴ * Vᴴ') + @test isposdef(S) - U, S, Vᴴ = svd_compact(A; alg) - @test U isa CuMatrix{T} && size(U) == (m, minmn) - @test S isa Diagonal{real(T), <:CuVector} && size(S) == (minmn, minmn) - @test Vᴴ isa CuMatrix{T} && size(Vᴴ) == (minmn, n) - @test U * S * Vᴴ ≈ A - @test isapproxone(U' * U) - @test isapproxone(Vᴴ * Vᴴ') - @test isposdef(S) + Ac = similar(A) + U2, S2, V2ᴴ = @constinferred svd_compact!(copy!(Ac, A), (U, S, Vᴴ), alg) + @test U2 === U + @test S2 === S + @test V2ᴴ === Vᴴ + @test U * S * Vᴴ ≈ A + @test isapproxone(U' * U) + @test isapproxone(Vᴴ * Vᴴ') + @test isposdef(S) - Ac = similar(A) - U2, S2, V2ᴴ = @constinferred svd_compact!(copy!(Ac, A), (U, S, Vᴴ), alg) - @test U2 === U - @test S2 === S - @test V2ᴴ === Vᴴ - @test U * S * Vᴴ ≈ A - @test isapproxone(U' * U) - @test isapproxone(Vᴴ * Vᴴ') - @test isposdef(S) - - Sd = svd_vals(A, alg) - @test CuArray(diagview(S)) ≈ Sd - # CuArray is necessary because norm of CuArray view with non-unit step is broken - if alg isa CUSOLVER_QRIteration - @test_warn "invalid keyword arguments for GPU_QRIteration" svd_compact!(copy!(Ac, A), (U, S, Vᴴ), CUSOLVER_QRIteration(; bad = "bad")) + Sd = svd_vals(A, alg) + @test CuArray(diagview(S)) ≈ Sd + # CuArray is necessary because norm of CuArray view with non-unit step is broken + if alg isa CUSOLVER_QRIteration + @test_warn "invalid keyword arguments for GPU_QRIteration" svd_compact!(copy!(Ac, A), (U, S, Vᴴ), CUSOLVER_QRIteration(; bad = "bad")) + end end end end @@ -53,56 +54,61 @@ end @testset "size ($m, $n)" for n in (37, m, 63) algs = (CUSOLVER_QRIteration(), CUSOLVER_SVDPolar(), CUSOLVER_Jacobi()) @testset "algorithm $alg" for alg in algs - A = CuArray(randn(rng, T, m, n)) - U, S, Vᴴ = svd_full(A; alg) - @test U isa CuMatrix{T} && size(U) == (m, m) - @test S isa CuMatrix{real(T)} && size(S) == (m, n) - @test Vᴴ isa CuMatrix{T} && size(Vᴴ) == (n, n) - @test U * S * Vᴴ ≈ A - @test isapproxone(U' * U) - @test isapproxone(U * U') - @test isapproxone(Vᴴ * Vᴴ') - @test isapproxone(Vᴴ' * Vᴴ) - @test all(isposdef, diagview(S)) + As = m == n ? (CuArray(randn(rng, T, m, n)), Diagonal(CuArray(randn(rng, T, m)))) : (CuArray(randn(rng, T, m, n)),) + for A in As + minmn = min(m, n) + U, S, Vᴴ = svd_full(A; alg) + @test U isa CuMatrix{T} && size(U) == (m, m) + @test S isa CuMatrix{real(T)} && size(S) == (m, n) + @test Vᴴ isa CuMatrix{T} && size(Vᴴ) == (n, n) + @test U * S * Vᴴ ≈ A + @test isapproxone(U' * U) + @test isapproxone(U * U') + @test isapproxone(Vᴴ * Vᴴ') + @test isapproxone(Vᴴ' * Vᴴ) + @test all(isposdef, diagview(S)) - Ac = similar(A) - U2, S2, V2ᴴ = @constinferred svd_full!(copy!(Ac, A), (U, S, Vᴴ), alg) - @test U2 === U - @test S2 === S - @test V2ᴴ === Vᴴ - @test U * S * Vᴴ ≈ A - @test isapproxone(U' * U) - @test isapproxone(U * U') - @test isapproxone(Vᴴ * Vᴴ') - @test isapproxone(Vᴴ' * Vᴴ) - @test all(isposdef, diagview(S)) + Ac = similar(A) + U2, S2, V2ᴴ = @constinferred svd_full!(copy!(Ac, A), (U, S, Vᴴ), alg) + @test U2 === U + @test S2 === S + @test V2ᴴ === Vᴴ + @test U * S * Vᴴ ≈ A + @test isapproxone(U' * U) + @test isapproxone(U * U') + @test isapproxone(Vᴴ * Vᴴ') + @test isapproxone(Vᴴ' * Vᴴ) + @test all(isposdef, diagview(S)) - minmn = min(m, n) - Sc = similar(A, real(T), minmn) - Sc2 = svd_vals!(copy!(Ac, A), Sc, alg) - @test Sc === Sc2 - @test CuArray(diagview(S)) ≈ Sc - # CuArray is necessary because norm of CuArray view with non-unit step is broken - if alg isa CUSOLVER_QRIteration - @test_warn "invalid keyword arguments for GPU_QRIteration" svd_full!(copy!(Ac, A), (U, S, Vᴴ), CUSOLVER_QRIteration(; bad = "bad")) - @test_warn "invalid keyword arguments for GPU_QRIteration" svd_vals!(copy!(Ac, A), Sc, CUSOLVER_QRIteration(; bad = "bad")) + Sc = similar(A, real(T), minmn) + Sc2 = svd_vals!(copy!(Ac, A), Sc, alg) + @test Sc === Sc2 + @test CuArray(diagview(S)) ≈ Sc + # CuArray is necessary because norm of CuArray view with non-unit step is broken + if alg isa CUSOLVER_QRIteration + @test_warn "invalid keyword arguments for GPU_QRIteration" svd_full!(copy!(Ac, A), (U, S, Vᴴ), CUSOLVER_QRIteration(; bad = "bad")) + if !isa(A, Diagonal) + @test_warn "invalid keyword arguments for GPU_QRIteration" svd_vals!(copy!(Ac, A), Sc, CUSOLVER_QRIteration(; bad = "bad")) + end + end end end end @testset "size (0, 0)" begin algs = (CUSOLVER_QRIteration(), CUSOLVER_SVDPolar(), CUSOLVER_Jacobi()) @testset "algorithm $alg" for alg in algs - A = CuArray(randn(rng, T, 0, 0)) - U, S, Vᴴ = svd_full(A; alg) - @test U isa CuMatrix{T} && size(U) == (0, 0) - @test S isa CuMatrix{real(T)} && size(S) == (0, 0) - @test Vᴴ isa CuMatrix{T} && size(Vᴴ) == (0, 0) - @test U * S * Vᴴ ≈ A - @test isapproxone(U' * U) - @test isapproxone(U * U') - @test isapproxone(Vᴴ * Vᴴ') - @test isapproxone(Vᴴ' * Vᴴ) - @test all(isposdef, diagview(S)) + for A in (CuArray(randn(rng, T, 0, 0)), Diagonal(CuArray(randn(rng, T, 0)))) + U, S, Vᴴ = svd_full(A; alg) + @test U isa CuMatrix{T} && size(U) == (0, 0) + @test S isa CuMatrix{real(T)} && size(S) == (0, 0) + @test Vᴴ isa CuMatrix{T} && size(Vᴴ) == (0, 0) + @test U * S * Vᴴ ≈ A + @test isapproxone(U' * U) + @test isapproxone(U * U') + @test isapproxone(Vᴴ * Vᴴ') + @test isapproxone(Vᴴ' * Vᴴ) + @test all(isposdef, diagview(S)) + end end end end @@ -115,26 +121,28 @@ end p = min(m, n) - k - 1 algs = (CUSOLVER_QRIteration(), CUSOLVER_SVDPolar(), CUSOLVER_Jacobi(), CUSOLVER_Randomized(; k = k, p = p, niters = 100)) @testset "algorithm $alg" for alg in algs - hA = randn(rng, T, m, n) - S₀ = svd_vals(hA) - A = CuArray(hA) - minmn = min(m, n) - r = k + hAs = m == n ? (randn(rng, T, m, n), Diagonal(randn(rng, T, m))) : (randn(rng, T, m, n),) + for hA in hAs + S₀ = svd_vals(hA) + A = CuArray(hA) + minmn = min(m, n) + r = k - U1, S1, V1ᴴ, ϵ1 = @constinferred svd_trunc(A; alg, trunc = truncrank(r)) - @test length(S1.diag) == r - @test opnorm(A - U1 * S1 * V1ᴴ) ≈ S₀[r + 1] - @test norm(A - U1 * S1 * V1ᴴ) ≈ ϵ1 + U1, S1, V1ᴴ, ϵ1 = @constinferred svd_trunc(A; alg, trunc = truncrank(r)) + @test length(S1.diag) == r + @test opnorm(A - U1 * S1 * V1ᴴ) ≈ S₀[r + 1] + @test norm(A - U1 * S1 * V1ᴴ) ≈ ϵ1 - if !(alg isa CUSOLVER_Randomized) - s = 1 + sqrt(eps(real(T))) - trunc2 = trunctol(; atol = s * S₀[r + 1]) + if !(alg isa CUSOLVER_Randomized) + s = 1 + sqrt(eps(real(T))) + trunc2 = trunctol(; atol = s * S₀[r + 1]) - U2, S2, V2ᴴ, ϵ2 = @constinferred svd_trunc(A; alg, trunc = trunctol(; atol = s * S₀[r + 1])) - @test length(S2.diag) == r - @test U1 ≈ U2 - @test parent(S1) ≈ parent(S2) - @test V1ᴴ ≈ V2ᴴ + U2, S2, V2ᴴ, ϵ2 = @constinferred svd_trunc(A; alg, trunc = trunctol(; atol = s * S₀[r + 1])) + @test length(S2.diag) == r + @test U1 ≈ U2 + @test parent(S1) ≈ parent(S2) + @test V1ᴴ ≈ V2ᴴ + end end end end From e83f6cc1200b22c1fe36d9bb51d25ca9bb8f752e Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 1 Dec 2025 16:21:49 -0500 Subject: [PATCH 06/14] remove `@invoke` calls --- .../MatrixAlgebraKitAMDGPUExt.jl | 11 +++++++---- .../MatrixAlgebraKitCUDAExt.jl | 17 +++++++++-------- 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl index bd9d31f0..c1933c61 100644 --- a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl +++ b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl @@ -132,20 +132,23 @@ MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where function MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T <: Complex} return all(isreal.(A.diag)) end -MatrixAlgebraKit.ishermitian_approx(A::StridedROCMatrix; kwargs...) = - @invoke MatrixAlgebraKit.ishermitian_approx(A::Any; kwargs...) MatrixAlgebraKit.isantihermitian_exact(A::StridedROCMatrix) = all(A .== -adjoint(A)) MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T <: Real} = all(iszero(A.diag)) MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T <: Complex} = all(iszero.(real.(A.diag))) -MatrixAlgebraKit.isantihermitian_approx(A::StridedROCMatrix; kwargs...) = - @invoke MatrixAlgebraKit.isantihermitian_approx(A::Any; kwargs...) function MatrixAlgebraKit.isantihermitian_approx(A::Diagonal{T, <:StridedROCVector{T}}; atol, rtol, kwargs...) where {T <: Real} return norm(A) ≤ atol end +# avoids calling the `StridedMatrix` specialization to avoid scalar indexing, +# use (allocating) fallback instead until we write a dedicated kernel +MatrixAlgebraKit.ishermitian_approx(A::StridedROCMatrix; kwargs...) = + norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) +MatrixAlgebraKit.isantihermitian_approx(A::StridedROCMatrix; kwargs...) = + norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) + function MatrixAlgebraKit._avgdiff!(A::StridedROCMatrix, B::StridedROCMatrix) axes(A) == axes(B) || throw(DimensionMismatch()) # COV_EXCL_START diff --git a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl index d2c6b3f3..f7fe5b84 100644 --- a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl +++ b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl @@ -151,26 +151,27 @@ function MatrixAlgebraKit._project_hermitian_diag!(A::StridedCuMatrix, B::Stride return nothing end -MatrixAlgebraKit.ishermitian_exact(A::StridedCuMatrix) = - all(A .== adjoint(A)) +MatrixAlgebraKit.ishermitian_exact(A::StridedCuMatrix) = all(A .== adjoint(A)) MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T <: Real} = true function MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T <: Complex} return all(isreal.(A.diag)) end -MatrixAlgebraKit.ishermitian_approx(A::StridedCuMatrix; kwargs...) = - @invoke MatrixAlgebraKit.ishermitian_approx(A::Any; kwargs...) -MatrixAlgebraKit.isantihermitian_exact(A::StridedCuMatrix) = - all(A .== -adjoint(A)) +MatrixAlgebraKit.isantihermitian_exact(A::StridedCuMatrix) = all(A .== -adjoint(A)) MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T <: Real} = all(iszero(A.diag)) MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T <: Complex} = all(iszero.(real.(A.diag))) -MatrixAlgebraKit.isantihermitian_approx(A::StridedCuMatrix; kwargs...) = - @invoke MatrixAlgebraKit.isantihermitian_approx(A::Any; kwargs...) function MatrixAlgebraKit.isantihermitian_approx(A::Diagonal{T, <:StridedCuVector{T}}; atol, rtol, kwargs...) where {T <: Real} return norm(A) ≤ atol end +# avoids calling the `StridedMatrix` specialization to avoid scalar indexing, +# use (allocating) fallback instead until we write a dedicated kernel +MatrixAlgebraKit.ishermitian_approx(A::StridedCuMatrix; kwargs...) = + norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) +MatrixAlgebraKit.isantihermitian_approx(A::StridedCuMatrix; kwargs...) = + norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) + function MatrixAlgebraKit._avgdiff!(A::StridedCuMatrix, B::StridedCuMatrix) axes(A) == axes(B) || throw(DimensionMismatch()) # COV_EXCL_START From 73828badebefebe2fd8983682ae062034325731a Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 1 Dec 2025 16:38:23 -0500 Subject: [PATCH 07/14] fix kwargs --- ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl | 4 ++-- ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl index c1933c61..673fdc12 100644 --- a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl +++ b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl @@ -144,9 +144,9 @@ end # avoids calling the `StridedMatrix` specialization to avoid scalar indexing, # use (allocating) fallback instead until we write a dedicated kernel -MatrixAlgebraKit.ishermitian_approx(A::StridedROCMatrix; kwargs...) = +MatrixAlgebraKit.ishermitian_approx(A::StridedROCMatrix; atol, rtol, kwargs...) = norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) -MatrixAlgebraKit.isantihermitian_approx(A::StridedROCMatrix; kwargs...) = +MatrixAlgebraKit.isantihermitian_approx(A::StridedROCMatrix; atol, rtol, kwargs...) = norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) function MatrixAlgebraKit._avgdiff!(A::StridedROCMatrix, B::StridedROCMatrix) diff --git a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl index f7fe5b84..4cc30d67 100644 --- a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl +++ b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl @@ -167,9 +167,9 @@ end # avoids calling the `StridedMatrix` specialization to avoid scalar indexing, # use (allocating) fallback instead until we write a dedicated kernel -MatrixAlgebraKit.ishermitian_approx(A::StridedCuMatrix; kwargs...) = +MatrixAlgebraKit.ishermitian_approx(A::StridedCuMatrix; atol, rtol, kwargs...) = norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) -MatrixAlgebraKit.isantihermitian_approx(A::StridedCuMatrix; kwargs...) = +MatrixAlgebraKit.isantihermitian_approx(A::StridedCuMatrix; atol, rtol, kwargs...) = norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) function MatrixAlgebraKit._avgdiff!(A::StridedCuMatrix, B::StridedCuMatrix) From 3a3f05b0c4cee7615403687545fe0341a7ea39a5 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 1 Dec 2025 16:53:47 -0500 Subject: [PATCH 08/14] further cleanup of `ishermitian` and friends --- .../MatrixAlgebraKitAMDGPUExt.jl | 17 ++--------- .../MatrixAlgebraKitCUDAExt.jl | 16 ++--------- src/common/matrixproperties.jl | 28 ++++++++++++++----- 3 files changed, 25 insertions(+), 36 deletions(-) diff --git a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl index 673fdc12..258b46d6 100644 --- a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl +++ b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl @@ -127,25 +127,12 @@ function MatrixAlgebraKit._project_hermitian_diag!(A::StridedROCMatrix, B::Strid return nothing end -MatrixAlgebraKit.ishermitian_exact(A::StridedROCMatrix) = all(A .== adjoint(A)) -MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T <: Real} = true -function MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T <: Complex} - return all(isreal.(A.diag)) -end - -MatrixAlgebraKit.isantihermitian_exact(A::StridedROCMatrix) = - all(A .== -adjoint(A)) -MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T <: Real} = all(iszero(A.diag)) -MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T <: Complex} = - all(iszero.(real.(A.diag))) -function MatrixAlgebraKit.isantihermitian_approx(A::Diagonal{T, <:StridedROCVector{T}}; atol, rtol, kwargs...) where {T <: Real} - return norm(A) ≤ atol -end - # avoids calling the `StridedMatrix` specialization to avoid scalar indexing, # use (allocating) fallback instead until we write a dedicated kernel +MatrixAlgebraKit.ishermitian_exact(A::StridedROCMatrix) = A == A' MatrixAlgebraKit.ishermitian_approx(A::StridedROCMatrix; atol, rtol, kwargs...) = norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) +MatrixAlgebraKit.isantihermitian_exact(A::StridedROCMatrix) = A == -A' MatrixAlgebraKit.isantihermitian_approx(A::StridedROCMatrix; atol, rtol, kwargs...) = norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) diff --git a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl index 4cc30d67..189f5825 100644 --- a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl +++ b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl @@ -151,24 +151,12 @@ function MatrixAlgebraKit._project_hermitian_diag!(A::StridedCuMatrix, B::Stride return nothing end -MatrixAlgebraKit.ishermitian_exact(A::StridedCuMatrix) = all(A .== adjoint(A)) -MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T <: Real} = true -function MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T <: Complex} - return all(isreal.(A.diag)) -end - -MatrixAlgebraKit.isantihermitian_exact(A::StridedCuMatrix) = all(A .== -adjoint(A)) -MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T <: Real} = all(iszero(A.diag)) -MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T <: Complex} = - all(iszero.(real.(A.diag))) -function MatrixAlgebraKit.isantihermitian_approx(A::Diagonal{T, <:StridedCuVector{T}}; atol, rtol, kwargs...) where {T <: Real} - return norm(A) ≤ atol -end - # avoids calling the `StridedMatrix` specialization to avoid scalar indexing, # use (allocating) fallback instead until we write a dedicated kernel +MatrixAlgebraKit.ishermitian_exact(A::StridedCuMatrix) = A == A' MatrixAlgebraKit.ishermitian_approx(A::StridedCuMatrix; atol, rtol, kwargs...) = norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) +MatrixAlgebraKit.isantihermitian_exact(A::StridedCuMatrix) = A == -A' MatrixAlgebraKit.isantihermitian_approx(A::StridedCuMatrix; atol, rtol, kwargs...) = norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) diff --git a/src/common/matrixproperties.jl b/src/common/matrixproperties.jl index 925f92e6..6a384437 100644 --- a/src/common/matrixproperties.jl +++ b/src/common/matrixproperties.jl @@ -79,10 +79,13 @@ end ishermitian_exact(A) = A == A' ishermitian_exact(A::StridedMatrix; kwargs...) = strided_ishermitian_exact(A, Val(false); kwargs...) +ishermitian_exact(A::Diagonal; kwargs...) = diagonal_ishermitian_exact(A, Val(false); kwargs...) + function ishermitian_approx(A; atol, rtol, kwargs...) return norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) end ishermitian_approx(A::StridedMatrix; kwargs...) = strided_ishermitian_approx(A, Val(false); kwargs...) +ishermitian_approx(A::Diagonal; kwargs...) = diagonal_ishermitian_approx(A, Val(false); kwargs...) """ isantihermitian(A; isapprox_kwargs...) @@ -97,16 +100,15 @@ function isantihermitian(A; atol::Real = 0, rtol::Real = 0, kwargs...) return isantihermitian_approx(A; atol, rtol, kwargs...) end end -function isantihermitian_exact(A) - return A == -A' -end -function isantihermitian_exact(A::StridedMatrix; kwargs...) - return strided_ishermitian_exact(A, Val(true); kwargs...) -end +isantihermitian_exact(A) = A == -A' +isantihermitian_exact(A::StridedMatrix; kwargs...) = strided_ishermitian_exact(A, Val(true); kwargs...) +isantihermitian_exact(A::Diagonal; kwargs...) = diagonal_ishermitian_exact(A, Val(true); kwargs...) + function isantihermitian_approx(A; atol, rtol, kwargs...) return norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) end isantihermitian_approx(A::StridedMatrix; kwargs...) = strided_ishermitian_approx(A, Val(true); kwargs...) +isantihermitian_approx(A::Diagonal; kwargs...) = diagonal_ishermitian_approx(A, Val(true); kwargs...) # blocked implementation of exact checks for strided matrices # ----------------------------------------------------------- @@ -145,7 +147,6 @@ function _ishermitian_exact_offdiag(Al, Au, ::Val{anti}) where {anti} return true end - function strided_ishermitian_approx( A::AbstractMatrix, anti::Val; blocksize = 32, atol::Real = default_hermitian_tol(A), rtol::Real = 0 @@ -192,3 +193,16 @@ function _ishermitian_approx_offdiag(Al, Au, ::Val{anti}) where {anti} end return ϵ² end + +diagonal_ishermitian_exact(A, ::Val{anti}) where {anti} = all(iszero ∘ (anti ? real : imag), diagview(A)) + +function diagonal_ishermitian_approx( + A, ::Val{anti}; atol::Real = default_hermitian_tol(A), rtol::Real = 0 + ) where {anti} + m, n = size(A) + m == n || throw(DimensionMismatch()) + init = abs2(zero(eltype(A))) + ϵ² = sum(abs2 ∘ (anti ? real : imag), diagview(A); init) + ϵ²max = oftype(ϵ², rtol > 0 ? max(atol, rtol * norm(A)) : atol)^2 + return ϵ² ≤ ϵ²max +end From 8a402bda7aa1809500624dcc2ea0f86b8fe3b06c Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 1 Dec 2025 17:15:59 -0500 Subject: [PATCH 09/14] include diagonal tests in projections on CPU --- test/amd/projections.jl | 12 +++------ test/cuda/projections.jl | 12 +++------ test/projections.jl | 57 +++++++++++++++++++++------------------- 3 files changed, 38 insertions(+), 43 deletions(-) diff --git a/test/amd/projections.jl b/test/amd/projections.jl index 2e449be3..0236aa66 100644 --- a/test/amd/projections.jl +++ b/test/amd/projections.jl @@ -21,23 +21,19 @@ const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) @test ishermitian(Bh) @test Bh ≈ Ah @test A == Ac - # this is still hermitian for real Diagonals! Bh_approx = Bh + noisefactor * Aa - if !isa(A, Diagonal) && !(T <: Real) - @test !ishermitian(Bh_approx) - end + # this is still hermitian for real Diagonal: |A - A'| == 0 + @test !ishermitian(Bh_approx) || norm(Aa) == 0 @test ishermitian(Bh_approx; rtol = 10 * noisefactor) Ba = project_antihermitian(A, alg) @test isantihermitian(Ba) @test Ba ≈ Aa @test A == Ac - # this is still hermitian for real Diagonals! Ba_approx = Ba + noisefactor * Ah @test !isantihermitian(Ba_approx) - if !isa(A, Diagonal) && !(T <: Real) - @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) - end + # this is never anti-hermitian for real Diagonal: |A - A'| == 0 + @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) || norm(Aa) == 0 Bh = project_hermitian!(Ac, alg) @test Bh === Ac diff --git a/test/cuda/projections.jl b/test/cuda/projections.jl index b7fc42b1..d157ff5f 100644 --- a/test/cuda/projections.jl +++ b/test/cuda/projections.jl @@ -21,11 +21,9 @@ const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) @test ishermitian(Bh) @test Bh ≈ Ah @test A == Ac - # this is still hermitian for real Diagonals! Bh_approx = Bh + noisefactor * Aa - if !isa(A, Diagonal) && !(T <: Real) - @test !ishermitian(Bh_approx) - end + # this is still hermitian for real Diagonal: |A - A'| == 0 + @test !ishermitian(Bh_approx) || norm(Aa) == 0 @test ishermitian(Bh_approx; rtol = 10 * noisefactor) Ba = project_antihermitian(A, alg) @@ -34,10 +32,8 @@ const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) @test A == Ac Ba_approx = Ba + noisefactor * Ah @test !isantihermitian(Ba_approx) - # this is still hermitian for real Diagonals! - if !isa(A, Diagonal) && !(T <: Real) - @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) - end + # this is never anti-hermitian for real Diagonal: |A - A'| == 0 + @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) || norm(Aa) == 0 Bh = project_hermitian!(Ac, alg) @test Bh === Ac diff --git a/test/projections.jl b/test/projections.jl index 56bc1a04..3923528e 100644 --- a/test/projections.jl +++ b/test/projections.jl @@ -11,37 +11,40 @@ const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) m = 54 noisefactor = eps(real(T))^(3 / 4) for alg in (NativeBlocked(blocksize = 16), NativeBlocked(blocksize = 32), NativeBlocked(blocksize = 64)) - A = randn(rng, T, m, m) - Ah = (A + A') / 2 - Aa = (A - A') / 2 - Ac = copy(A) + for A in (randn(rng, T, m, m), Diagonal(randn(rng, T, m))) + Ah = (A + A') / 2 + Aa = (A - A') / 2 + Ac = copy(A) - Bh = project_hermitian(A, alg) - @test ishermitian(Bh) - @test Bh ≈ Ah - @test A == Ac - Bh_approx = Bh + noisefactor * Aa - @test !ishermitian(Bh_approx) - @test ishermitian(Bh_approx; rtol = 10 * noisefactor) + Bh = project_hermitian(A, alg) + @test ishermitian(Bh) + @test Bh ≈ Ah + @test A == Ac + Bh_approx = Bh + noisefactor * Aa + # this is still hermitian for real Diagonal: |A - A'| == 0 + @test !ishermitian(Bh_approx) || norm(Aa) == 0 + @test ishermitian(Bh_approx; rtol = 10 * noisefactor) - Ba = project_antihermitian(A, alg) - @test isantihermitian(Ba) - @test Ba ≈ Aa - @test A == Ac - Ba_approx = Ba + noisefactor * Ah - @test !isantihermitian(Ba_approx) - @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) + Ba = project_antihermitian(A, alg) + @test isantihermitian(Ba) + @test Ba ≈ Aa + @test A == Ac + Ba_approx = Ba + noisefactor * Ah + @test !isantihermitian(Ba_approx) + # this is never anti-hermitian for real Diagonal: |A - A'| == 0 + @test isantihermitian(Ba_approx; rtol = 10 * noisefactor) || norm(Aa) == 0 - Bh = project_hermitian!(Ac, alg) - @test Bh === Ac - @test ishermitian(Bh) - @test Bh ≈ Ah + Bh = project_hermitian!(Ac, alg) + @test Bh === Ac + @test ishermitian(Bh) + @test Bh ≈ Ah - copy!(Ac, A) - Ba = project_antihermitian!(Ac, alg) - @test Ba === Ac - @test isantihermitian(Ba) - @test Ba ≈ Aa + copy!(Ac, A) + Ba = project_antihermitian!(Ac, alg) + @test Ba === Ac + @test isantihermitian(Ba) + @test Ba ≈ Aa + end end # test approximate error calculation From 59ac47575e3a84c47cc36a0b70f99c0ff4bb57b3 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 1 Dec 2025 18:33:19 -0500 Subject: [PATCH 10/14] make JET happy --- src/common/matrixproperties.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/common/matrixproperties.jl b/src/common/matrixproperties.jl index 6a384437..5a73176e 100644 --- a/src/common/matrixproperties.jl +++ b/src/common/matrixproperties.jl @@ -79,7 +79,7 @@ end ishermitian_exact(A) = A == A' ishermitian_exact(A::StridedMatrix; kwargs...) = strided_ishermitian_exact(A, Val(false); kwargs...) -ishermitian_exact(A::Diagonal; kwargs...) = diagonal_ishermitian_exact(A, Val(false); kwargs...) +ishermitian_exact(A::Diagonal) = diagonal_ishermitian_exact(A, Val(false)) function ishermitian_approx(A; atol, rtol, kwargs...) return norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) @@ -102,7 +102,7 @@ function isantihermitian(A; atol::Real = 0, rtol::Real = 0, kwargs...) end isantihermitian_exact(A) = A == -A' isantihermitian_exact(A::StridedMatrix; kwargs...) = strided_ishermitian_exact(A, Val(true); kwargs...) -isantihermitian_exact(A::Diagonal; kwargs...) = diagonal_ishermitian_exact(A, Val(true); kwargs...) +isantihermitian_exact(A::Diagonal) = diagonal_ishermitian_exact(A, Val(true)) function isantihermitian_approx(A; atol, rtol, kwargs...) return norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A)) From 978ae4946e33f92599284acb5505b4dff3e3b183 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 1 Dec 2025 18:41:55 -0500 Subject: [PATCH 11/14] Revert "Add tests for SVD" This reverts commit 1b4399f52a8e132a7aab2a685e6726734d5a8298. --- test/cuda/svd.jl | 176 ++++++++++++++++++++++------------------------- 1 file changed, 84 insertions(+), 92 deletions(-) diff --git a/test/cuda/svd.jl b/test/cuda/svd.jl index 6a2f6e1e..e2173d03 100644 --- a/test/cuda/svd.jl +++ b/test/cuda/svd.jl @@ -15,34 +15,33 @@ include(joinpath("..", "utilities.jl")) k = min(m, n) algs = (CUSOLVER_QRIteration(), CUSOLVER_SVDPolar(), CUSOLVER_Jacobi()) @testset "algorithm $alg" for alg in algs - As = m == n ? (CuArray(randn(rng, T, m, n)), Diagonal(CuArray(randn(rng, T, m)))) : (CuArray(randn(rng, T, m, n)),) - for A in As - minmn = min(m, n) - U, S, Vᴴ = svd_compact(A; alg) - @test U isa CuMatrix{T} && size(U) == (m, minmn) - @test S isa Diagonal{real(T), <:CuVector} && size(S) == (minmn, minmn) - @test Vᴴ isa CuMatrix{T} && size(Vᴴ) == (minmn, n) - @test U * S * Vᴴ ≈ A - @test isapproxone(U' * U) - @test isapproxone(Vᴴ * Vᴴ') - @test isposdef(S) + minmn = min(m, n) + A = CuArray(randn(rng, T, m, n)) - Ac = similar(A) - U2, S2, V2ᴴ = @constinferred svd_compact!(copy!(Ac, A), (U, S, Vᴴ), alg) - @test U2 === U - @test S2 === S - @test V2ᴴ === Vᴴ - @test U * S * Vᴴ ≈ A - @test isapproxone(U' * U) - @test isapproxone(Vᴴ * Vᴴ') - @test isposdef(S) + U, S, Vᴴ = svd_compact(A; alg) + @test U isa CuMatrix{T} && size(U) == (m, minmn) + @test S isa Diagonal{real(T), <:CuVector} && size(S) == (minmn, minmn) + @test Vᴴ isa CuMatrix{T} && size(Vᴴ) == (minmn, n) + @test U * S * Vᴴ ≈ A + @test isapproxone(U' * U) + @test isapproxone(Vᴴ * Vᴴ') + @test isposdef(S) - Sd = svd_vals(A, alg) - @test CuArray(diagview(S)) ≈ Sd - # CuArray is necessary because norm of CuArray view with non-unit step is broken - if alg isa CUSOLVER_QRIteration - @test_warn "invalid keyword arguments for GPU_QRIteration" svd_compact!(copy!(Ac, A), (U, S, Vᴴ), CUSOLVER_QRIteration(; bad = "bad")) - end + Ac = similar(A) + U2, S2, V2ᴴ = @constinferred svd_compact!(copy!(Ac, A), (U, S, Vᴴ), alg) + @test U2 === U + @test S2 === S + @test V2ᴴ === Vᴴ + @test U * S * Vᴴ ≈ A + @test isapproxone(U' * U) + @test isapproxone(Vᴴ * Vᴴ') + @test isposdef(S) + + Sd = svd_vals(A, alg) + @test CuArray(diagview(S)) ≈ Sd + # CuArray is necessary because norm of CuArray view with non-unit step is broken + if alg isa CUSOLVER_QRIteration + @test_warn "invalid keyword arguments for GPU_QRIteration" svd_compact!(copy!(Ac, A), (U, S, Vᴴ), CUSOLVER_QRIteration(; bad = "bad")) end end end @@ -54,61 +53,56 @@ end @testset "size ($m, $n)" for n in (37, m, 63) algs = (CUSOLVER_QRIteration(), CUSOLVER_SVDPolar(), CUSOLVER_Jacobi()) @testset "algorithm $alg" for alg in algs - As = m == n ? (CuArray(randn(rng, T, m, n)), Diagonal(CuArray(randn(rng, T, m)))) : (CuArray(randn(rng, T, m, n)),) - for A in As - minmn = min(m, n) - U, S, Vᴴ = svd_full(A; alg) - @test U isa CuMatrix{T} && size(U) == (m, m) - @test S isa CuMatrix{real(T)} && size(S) == (m, n) - @test Vᴴ isa CuMatrix{T} && size(Vᴴ) == (n, n) - @test U * S * Vᴴ ≈ A - @test isapproxone(U' * U) - @test isapproxone(U * U') - @test isapproxone(Vᴴ * Vᴴ') - @test isapproxone(Vᴴ' * Vᴴ) - @test all(isposdef, diagview(S)) + A = CuArray(randn(rng, T, m, n)) + U, S, Vᴴ = svd_full(A; alg) + @test U isa CuMatrix{T} && size(U) == (m, m) + @test S isa CuMatrix{real(T)} && size(S) == (m, n) + @test Vᴴ isa CuMatrix{T} && size(Vᴴ) == (n, n) + @test U * S * Vᴴ ≈ A + @test isapproxone(U' * U) + @test isapproxone(U * U') + @test isapproxone(Vᴴ * Vᴴ') + @test isapproxone(Vᴴ' * Vᴴ) + @test all(isposdef, diagview(S)) - Ac = similar(A) - U2, S2, V2ᴴ = @constinferred svd_full!(copy!(Ac, A), (U, S, Vᴴ), alg) - @test U2 === U - @test S2 === S - @test V2ᴴ === Vᴴ - @test U * S * Vᴴ ≈ A - @test isapproxone(U' * U) - @test isapproxone(U * U') - @test isapproxone(Vᴴ * Vᴴ') - @test isapproxone(Vᴴ' * Vᴴ) - @test all(isposdef, diagview(S)) + Ac = similar(A) + U2, S2, V2ᴴ = @constinferred svd_full!(copy!(Ac, A), (U, S, Vᴴ), alg) + @test U2 === U + @test S2 === S + @test V2ᴴ === Vᴴ + @test U * S * Vᴴ ≈ A + @test isapproxone(U' * U) + @test isapproxone(U * U') + @test isapproxone(Vᴴ * Vᴴ') + @test isapproxone(Vᴴ' * Vᴴ) + @test all(isposdef, diagview(S)) - Sc = similar(A, real(T), minmn) - Sc2 = svd_vals!(copy!(Ac, A), Sc, alg) - @test Sc === Sc2 - @test CuArray(diagview(S)) ≈ Sc - # CuArray is necessary because norm of CuArray view with non-unit step is broken - if alg isa CUSOLVER_QRIteration - @test_warn "invalid keyword arguments for GPU_QRIteration" svd_full!(copy!(Ac, A), (U, S, Vᴴ), CUSOLVER_QRIteration(; bad = "bad")) - if !isa(A, Diagonal) - @test_warn "invalid keyword arguments for GPU_QRIteration" svd_vals!(copy!(Ac, A), Sc, CUSOLVER_QRIteration(; bad = "bad")) - end - end + minmn = min(m, n) + Sc = similar(A, real(T), minmn) + Sc2 = svd_vals!(copy!(Ac, A), Sc, alg) + @test Sc === Sc2 + @test CuArray(diagview(S)) ≈ Sc + # CuArray is necessary because norm of CuArray view with non-unit step is broken + if alg isa CUSOLVER_QRIteration + @test_warn "invalid keyword arguments for GPU_QRIteration" svd_full!(copy!(Ac, A), (U, S, Vᴴ), CUSOLVER_QRIteration(; bad = "bad")) + @test_warn "invalid keyword arguments for GPU_QRIteration" svd_vals!(copy!(Ac, A), Sc, CUSOLVER_QRIteration(; bad = "bad")) end end end @testset "size (0, 0)" begin algs = (CUSOLVER_QRIteration(), CUSOLVER_SVDPolar(), CUSOLVER_Jacobi()) @testset "algorithm $alg" for alg in algs - for A in (CuArray(randn(rng, T, 0, 0)), Diagonal(CuArray(randn(rng, T, 0)))) - U, S, Vᴴ = svd_full(A; alg) - @test U isa CuMatrix{T} && size(U) == (0, 0) - @test S isa CuMatrix{real(T)} && size(S) == (0, 0) - @test Vᴴ isa CuMatrix{T} && size(Vᴴ) == (0, 0) - @test U * S * Vᴴ ≈ A - @test isapproxone(U' * U) - @test isapproxone(U * U') - @test isapproxone(Vᴴ * Vᴴ') - @test isapproxone(Vᴴ' * Vᴴ) - @test all(isposdef, diagview(S)) - end + A = CuArray(randn(rng, T, 0, 0)) + U, S, Vᴴ = svd_full(A; alg) + @test U isa CuMatrix{T} && size(U) == (0, 0) + @test S isa CuMatrix{real(T)} && size(S) == (0, 0) + @test Vᴴ isa CuMatrix{T} && size(Vᴴ) == (0, 0) + @test U * S * Vᴴ ≈ A + @test isapproxone(U' * U) + @test isapproxone(U * U') + @test isapproxone(Vᴴ * Vᴴ') + @test isapproxone(Vᴴ' * Vᴴ) + @test all(isposdef, diagview(S)) end end end @@ -121,28 +115,26 @@ end p = min(m, n) - k - 1 algs = (CUSOLVER_QRIteration(), CUSOLVER_SVDPolar(), CUSOLVER_Jacobi(), CUSOLVER_Randomized(; k = k, p = p, niters = 100)) @testset "algorithm $alg" for alg in algs - hAs = m == n ? (randn(rng, T, m, n), Diagonal(randn(rng, T, m))) : (randn(rng, T, m, n),) - for hA in hAs - S₀ = svd_vals(hA) - A = CuArray(hA) - minmn = min(m, n) - r = k + hA = randn(rng, T, m, n) + S₀ = svd_vals(hA) + A = CuArray(hA) + minmn = min(m, n) + r = k - U1, S1, V1ᴴ, ϵ1 = @constinferred svd_trunc(A; alg, trunc = truncrank(r)) - @test length(S1.diag) == r - @test opnorm(A - U1 * S1 * V1ᴴ) ≈ S₀[r + 1] - @test norm(A - U1 * S1 * V1ᴴ) ≈ ϵ1 + U1, S1, V1ᴴ, ϵ1 = @constinferred svd_trunc(A; alg, trunc = truncrank(r)) + @test length(S1.diag) == r + @test opnorm(A - U1 * S1 * V1ᴴ) ≈ S₀[r + 1] + @test norm(A - U1 * S1 * V1ᴴ) ≈ ϵ1 - if !(alg isa CUSOLVER_Randomized) - s = 1 + sqrt(eps(real(T))) - trunc2 = trunctol(; atol = s * S₀[r + 1]) + if !(alg isa CUSOLVER_Randomized) + s = 1 + sqrt(eps(real(T))) + trunc2 = trunctol(; atol = s * S₀[r + 1]) - U2, S2, V2ᴴ, ϵ2 = @constinferred svd_trunc(A; alg, trunc = trunctol(; atol = s * S₀[r + 1])) - @test length(S2.diag) == r - @test U1 ≈ U2 - @test parent(S1) ≈ parent(S2) - @test V1ᴴ ≈ V2ᴴ - end + U2, S2, V2ᴴ, ϵ2 = @constinferred svd_trunc(A; alg, trunc = trunctol(; atol = s * S₀[r + 1])) + @test length(S2.diag) == r + @test U1 ≈ U2 + @test parent(S1) ≈ parent(S2) + @test V1ᴴ ≈ V2ᴴ end end end From 6bf485497d059a4da6ea7756c19b57bff2bcce53 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 1 Dec 2025 18:42:24 -0500 Subject: [PATCH 12/14] revert Diagonal to `diagm` for svd with gpuarrays --- src/implementations/svd.jl | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/implementations/svd.jl b/src/implementations/svd.jl index ef0c2bc3..c8b27f3f 100644 --- a/src/implementations/svd.jl +++ b/src/implementations/svd.jl @@ -356,7 +356,6 @@ function svd_full!(A::AbstractMatrix, USVᴴ, alg::GPU_SVDAlgorithm) return USVᴴ end -svd_full!(A::Diagonal, USVᴴ, alg::GPU_SVDAlgorithm) = svd_full!(diagm(A.diag), USVᴴ, alg) function svd_trunc!(A::AbstractMatrix, USVᴴ, alg::TruncatedAlgorithm{<:GPU_Randomized}) check_input(svd_trunc!, A, USVᴴ, alg.alg) @@ -374,7 +373,6 @@ function svd_trunc!(A::AbstractMatrix, USVᴴ, alg::TruncatedAlgorithm{<:GPU_Ran return Utr, Str, Vᴴtr, ϵ end -svd_trunc!(A::Diagonal, USVᴴ, alg::GPU_SVDAlgorithm) = svd_trunc!(diagm(A.diag), USVᴴ, alg) function svd_compact!(A::AbstractMatrix, USVᴴ, alg::GPU_SVDAlgorithm) check_input(svd_compact!, A, USVᴴ, alg) @@ -398,7 +396,6 @@ function svd_compact!(A::AbstractMatrix, USVᴴ, alg::GPU_SVDAlgorithm) return USVᴴ end -svd_compact!(A::Diagonal, USVᴴ, alg::GPU_SVDAlgorithm) = svd_compact!(diagm(A.diag), USVᴴ, alg) _argmaxabs(x) = reduce(_largest, x; init = zero(eltype(x))) _largest(x, y) = abs(x) < abs(y) ? y : x @@ -421,10 +418,3 @@ function svd_vals!(A::AbstractMatrix, S, alg::GPU_SVDAlgorithm) return S end -function svd_vals!(A::Diagonal, S, alg::GPU_SVDAlgorithm) - check_input(svd_vals!, A, S, alg) - Ad = diagview(A) - S .= abs.(Ad) - sort!(S; rev = true) - return S -end From d154baef9c8ed398aa9dae326959d6465ae76aca Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 1 Dec 2025 19:06:01 -0500 Subject: [PATCH 13/14] dont call project_isometric with invalid alg-input types --- test/amd/projections.jl | 53 +++++++++++++++++++++++++++------------ test/cuda/projections.jl | 54 ++++++++++++++++++++++++++++------------ 2 files changed, 75 insertions(+), 32 deletions(-) diff --git a/test/amd/projections.jl b/test/amd/projections.jl index 0236aa66..5a1e8a63 100644 --- a/test/amd/projections.jl +++ b/test/amd/projections.jl @@ -57,25 +57,46 @@ end svdalgs = (ROCSOLVER_QRIteration(), ROCSOLVER_Jacobi()) algs = (PolarViaSVD.(svdalgs)...,) # PolarNewton()) # TODO @testset "algorithm $alg" for alg in algs - for A in (ROCArray(randn(rng, T, m, n)), Diagonal(ROCArray(randn(rng, T, m)))) + A = ROCArray(randn(rng, T, m, n)) + W = project_isometric(A, alg) + @test isisometric(W) + W2 = project_isometric(W, alg) + @test W2 ≈ W # stability of the projection + @test W * (W' * A) ≈ A + + Ac = similar(A) + W2 = @constinferred project_isometric!(copy!(Ac, A), W, alg) + @test W2 === W + @test isisometric(W) + + # test that W is closer to A then any other isometry + for k in 1:10 + δA = ROCArray(randn(rng, T, size(A)...)) W = project_isometric(A, alg) - @test isisometric(W) - W2 = project_isometric(W, alg) - @test W2 ≈ W # stability of the projection - @test W * (W' * A) ≈ A + W2 = project_isometric(A + δA / 100, alg) + @test norm(A - W2) >= norm(A - W) + end + end + + m == n && @testset "DiagonalAlgorithm" begin + A = Diagonal(ROCArray(randn(rng, T, m))) + W = project_isometric(A, alg) + @test isisometric(W) + W2 = project_isometric(W, alg) + @test W2 ≈ W # stability of the projection + @test W * (W' * A) ≈ A - Ac = similar(A) - W2 = @constinferred project_isometric!(copy!(Ac, A), W, alg) - @test W2 === W - @test isisometric(W) + Ac = similar(A) + W2 = @constinferred project_isometric!(copy!(Ac, A), W, alg) + @test W2 === W + @test isisometric(W) - # test that W is closer to A then any other isometry - for k in 1:10 - δA = ROCArray(randn(rng, T, size(A)...)) - W = project_isometric(A, alg) - W2 = project_isometric(A + δA / 100, alg) - @test norm(A - W2) > norm(A - W) - end + # test that W is closer to A then any other isometry + for k in 1:10 + δA = Diagonal(ROCArray(randn(rng, T, m))) + W = project_isometric(A, alg) + W2 = project_isometric(A + δA / 100, alg) + @test norm(A - W2) >= norm(A - W) end end end diff --git a/test/cuda/projections.jl b/test/cuda/projections.jl index d157ff5f..677ef520 100644 --- a/test/cuda/projections.jl +++ b/test/cuda/projections.jl @@ -57,25 +57,47 @@ end svdalgs = (CUSOLVER_SVDPolar(), CUSOLVER_QRIteration(), CUSOLVER_Jacobi()) algs = (PolarViaSVD.(svdalgs)...,) # PolarNewton()) # TODO @testset "algorithm $alg" for alg in algs - for A in (CuArray(randn(rng, T, m, n)), Diagonal(CuArray(randn(rng, T, m)))) + A = CuArray(randn(rng, T, m, n)) + W = project_isometric(A, alg) + @test isisometric(W) + W2 = project_isometric(W, alg) + @test W2 ≈ W # stability of the projection + @test W * (W' * A) ≈ A + + Ac = similar(A) + W2 = @constinferred project_isometric!(copy!(Ac, A), W, alg) + @test W2 === W + @test isisometric(W) + + # test that W is closer to A then any other isometry + for k in 1:10 + δA = CuArray(randn(rng, T, size(A)...)) W = project_isometric(A, alg) - @test isisometric(W) - W2 = project_isometric(W, alg) - @test W2 ≈ W # stability of the projection - @test W * (W' * A) ≈ A + W2 = project_isometric(A + δA / 100, alg) + @test norm(A - W2) >= norm(A - W) + end + end + + m == n && @testset "DiagonalAlgorithm" begin + A = Diagonal(CuArray(randn(rng, T, m))) + alg = PolarViaSVD(DiagonalAlgorithm()) + W = project_isometric(A, alg) + @test isisometric(W) + W2 = project_isometric(W, alg) + @test W2 ≈ W # stability of the projection + @test W * (W' * A) ≈ A - Ac = similar(A) - W2 = @constinferred project_isometric!(copy!(Ac, A), W, alg) - @test W2 === W - @test isisometric(W) + Ac = similar(A) + W2 = @constinferred project_isometric!(copy!(Ac, A), W, alg) + @test W2 === W + @test isisometric(W) - # test that W is closer to A then any other isometry - for k in 1:10 - δA = CuArray(randn(rng, T, size(A)...)) - W = project_isometric(A, alg) - W2 = project_isometric(A + δA / 100, alg) - @test norm(A - W2) > norm(A - W) - end + # test that W is closer to A then any other isometry + for k in 1:10 + δA = Diagonal(CuArray(randn(rng, T, m))) + W = project_isometric(A, alg) + W2 = project_isometric(A + δA / 100, alg) + @test norm(A - W2) >= norm(A - W) end end end From 60461ab19e92999e77ae307daca2a9686ac1f4e8 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Tue, 2 Dec 2025 03:38:00 -0500 Subject: [PATCH 14/14] Add alg for amd test --- test/amd/projections.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/amd/projections.jl b/test/amd/projections.jl index 5a1e8a63..a06b152c 100644 --- a/test/amd/projections.jl +++ b/test/amd/projections.jl @@ -80,6 +80,7 @@ end m == n && @testset "DiagonalAlgorithm" begin A = Diagonal(ROCArray(randn(rng, T, m))) + alg = PolarViaSVD(DiagonalAlgorithm()) W = project_isometric(A, alg) @test isisometric(W) W2 = project_isometric(W, alg)