From 4becc5a87baaf27a9943c407f94f8e0273c06e1f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matteo=20Secl=C3=AC?= Date: Mon, 26 May 2025 17:28:36 +0200 Subject: [PATCH 01/11] Expand eigen() and add eig[vals,vecs]() - Expand LinearAlgebra.eigen() to handle non-symmetric matrices via recent `Xgeev` - Add LinearAlgebra.eigvals() and LinearAlgebra.eigvecs() --- lib/cusolver/linalg.jl | 53 +++++++++++++++++++++++++++++++++++++++--- 1 file changed, 50 insertions(+), 3 deletions(-) diff --git a/lib/cusolver/linalg.jl b/lib/cusolver/linalg.jl index 0f2d9fac49..d3a16f5c91 100644 --- a/lib/cusolver/linalg.jl +++ b/lib/cusolver/linalg.jl @@ -110,7 +110,7 @@ function Base.:\(F::Union{LinearAlgebra.LAPACKFactorizations{<:Any,<:CuArray}, return LinearAlgebra._cut_B(BB, 1:n) end -# eigenvalues +# eigen function LinearAlgebra.eigen(A::Symmetric{T,<:CuMatrix}) where {T<:BlasReal} A2 = copy(A.data) @@ -126,13 +126,60 @@ end function LinearAlgebra.eigen(A::CuMatrix{T}) where {T<:BlasReal} A2 = copy(A) - issymmetric(A) ? Eigen(syevd!('V', 'U', A2)...) : error("GPU eigensolver supports only Hermitian or Symmetric matrices.") + r = Xgeev!('N', 'V', A2) + Eigen(r[1], r[3]) end function LinearAlgebra.eigen(A::CuMatrix{T}) where {T<:BlasComplex} A2 = copy(A) - ishermitian(A) ? Eigen(heevd!('V', 'U', A2)...) : error("GPU eigensolver supports only Hermitian or Symmetric matrices.") + r = Xgeev!('N', 'V', A2) + Eigen(r[1], r[3]) end +# eigvals + +function LinearAlgebra.eigvals(A::Symmetric{T,<:CuMatrix}) where {T<:BlasReal} + A2 = copy(A.data) + syevd!('N', 'U', A2)[1] +end +function LinearAlgebra.eigvals(A::Hermitian{T,<:CuMatrix}) where {T<:BlasComplex} + A2 = copy(A.data) + heevd!('N', 'U', A2)[1] +end +function LinearAlgebra.eigvals(A::Hermitian{T,<:CuMatrix}) where {T<:BlasReal} + eigvals(Symmetric(A)) +end + +function LinearAlgebra.eigvals(A::CuMatrix{T}) where {T<:BlasReal} + A2 = copy(A) + Xgeev!('N', 'N', A2)[1] +end +function LinearAlgebra.eigvals(A::CuMatrix{T}) where {T<:BlasComplex} + A2 = copy(A) + Xgeev!('N', 'N', A2)[1] +end + +# eigvecs + +function LinearAlgebra.eigvecs(A::Symmetric{T,<:CuMatrix}) where {T<:BlasReal} + A2 = copy(A.data) + syevd!('V', 'U', A2)[2] +end +function LinearAlgebra.eigvecs(A::Hermitian{T,<:CuMatrix}) where {T<:BlasComplex} + A2 = copy(A.data) + heevd!('V', 'U', A2)[2] +end +function LinearAlgebra.eigvecs(A::Hermitian{T,<:CuMatrix}) where {T<:BlasReal} + eigvals(Symmetric(A)) +end + +function LinearAlgebra.eigvecs(A::CuMatrix{T}) where {T<:BlasReal} + A2 = copy(A) + Xgeev!('N', 'V', A2)[3] +end +function LinearAlgebra.eigvecs(A::CuMatrix{T}) where {T<:BlasComplex} + A2 = copy(A) + Xgeev!('N', 'V', A2)[3] +end # factorizations From 291ea99943db304b84de9432906b014b34fd3d7d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matteo=20Secl=C3=AC?= Date: Mon, 26 May 2025 17:38:14 +0200 Subject: [PATCH 02/11] Add returns --- lib/cusolver/linalg.jl | 50 +++++++++++++++++++++--------------------- 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/lib/cusolver/linalg.jl b/lib/cusolver/linalg.jl index d3a16f5c91..a5234bd499 100644 --- a/lib/cusolver/linalg.jl +++ b/lib/cusolver/linalg.jl @@ -114,71 +114,71 @@ end function LinearAlgebra.eigen(A::Symmetric{T,<:CuMatrix}) where {T<:BlasReal} A2 = copy(A.data) - Eigen(syevd!('V', 'U', A2)...) + return Eigen(syevd!('V', 'U', A2)...) end function LinearAlgebra.eigen(A::Hermitian{T,<:CuMatrix}) where {T<:BlasComplex} A2 = copy(A.data) - Eigen(heevd!('V', 'U', A2)...) + return Eigen(heevd!('V', 'U', A2)...) end function LinearAlgebra.eigen(A::Hermitian{T,<:CuMatrix}) where {T<:BlasReal} - eigen(Symmetric(A)) + return eigen(Symmetric(A)) end function LinearAlgebra.eigen(A::CuMatrix{T}) where {T<:BlasReal} A2 = copy(A) r = Xgeev!('N', 'V', A2) - Eigen(r[1], r[3]) + return Eigen(r[1], r[3]) end function LinearAlgebra.eigen(A::CuMatrix{T}) where {T<:BlasComplex} A2 = copy(A) r = Xgeev!('N', 'V', A2) - Eigen(r[1], r[3]) + return Eigen(r[1], r[3]) end # eigvals -function LinearAlgebra.eigvals(A::Symmetric{T,<:CuMatrix}) where {T<:BlasReal} +function LinearAlgebra.eigvals(A::Symmetric{T, <:CuMatrix}) where {T <: BlasReal} A2 = copy(A.data) - syevd!('N', 'U', A2)[1] + return syevd!('N', 'U', A2)[1] end -function LinearAlgebra.eigvals(A::Hermitian{T,<:CuMatrix}) where {T<:BlasComplex} +function LinearAlgebra.eigvals(A::Hermitian{T, <:CuMatrix}) where {T <: BlasComplex} A2 = copy(A.data) - heevd!('N', 'U', A2)[1] + return heevd!('N', 'U', A2)[1] end -function LinearAlgebra.eigvals(A::Hermitian{T,<:CuMatrix}) where {T<:BlasReal} - eigvals(Symmetric(A)) +function LinearAlgebra.eigvals(A::Hermitian{T, <:CuMatrix}) where {T <: BlasReal} + return eigvals(Symmetric(A)) end -function LinearAlgebra.eigvals(A::CuMatrix{T}) where {T<:BlasReal} +function LinearAlgebra.eigvals(A::CuMatrix{T}) where {T <: BlasReal} A2 = copy(A) - Xgeev!('N', 'N', A2)[1] + return Xgeev!('N', 'N', A2)[1] end -function LinearAlgebra.eigvals(A::CuMatrix{T}) where {T<:BlasComplex} +function LinearAlgebra.eigvals(A::CuMatrix{T}) where {T <: BlasComplex} A2 = copy(A) - Xgeev!('N', 'N', A2)[1] + return Xgeev!('N', 'N', A2)[1] end # eigvecs -function LinearAlgebra.eigvecs(A::Symmetric{T,<:CuMatrix}) where {T<:BlasReal} +function LinearAlgebra.eigvecs(A::Symmetric{T, <:CuMatrix}) where {T <: BlasReal} A2 = copy(A.data) - syevd!('V', 'U', A2)[2] + return syevd!('V', 'U', A2)[2] end -function LinearAlgebra.eigvecs(A::Hermitian{T,<:CuMatrix}) where {T<:BlasComplex} +function LinearAlgebra.eigvecs(A::Hermitian{T, <:CuMatrix}) where {T <: BlasComplex} A2 = copy(A.data) - heevd!('V', 'U', A2)[2] + return heevd!('V', 'U', A2)[2] end -function LinearAlgebra.eigvecs(A::Hermitian{T,<:CuMatrix}) where {T<:BlasReal} - eigvals(Symmetric(A)) +function LinearAlgebra.eigvecs(A::Hermitian{T, <:CuMatrix}) where {T <: BlasReal} + return eigvals(Symmetric(A)) end -function LinearAlgebra.eigvecs(A::CuMatrix{T}) where {T<:BlasReal} +function LinearAlgebra.eigvecs(A::CuMatrix{T}) where {T <: BlasReal} A2 = copy(A) - Xgeev!('N', 'V', A2)[3] + return Xgeev!('N', 'V', A2)[3] end -function LinearAlgebra.eigvecs(A::CuMatrix{T}) where {T<:BlasComplex} +function LinearAlgebra.eigvecs(A::CuMatrix{T}) where {T <: BlasComplex} A2 = copy(A) - Xgeev!('N', 'V', A2)[3] + return Xgeev!('N', 'V', A2)[3] end # factorizations From 11bf501e11250352afb477897e775a849aec240e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matteo=20Secl=C3=AC?= Date: Tue, 27 May 2025 14:49:18 +0200 Subject: [PATCH 03/11] Add `geev` tests for `eigen` --- test/libraries/cusolver/dense.jl | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/test/libraries/cusolver/dense.jl b/test/libraries/cusolver/dense.jl index 4df323226c..010b5306eb 100644 --- a/test/libraries/cusolver/dense.jl +++ b/test/libraries/cusolver/dense.jl @@ -315,6 +315,35 @@ k = 1 end end + @testset "geev!" begin + A = rand(elty,m,m) + d_A = CuArray(A) + local d_W, d_V + d_W, d_V = CUSOLVER.Xgeev!('N','V', d_A) + d_W_b, d_V_b = LAPACK.geev!('N','V', CuArray(A)) + @test d_W ≈ d_W_b + @test d_V ≈ d_V_b + h_W = collect(d_W) + h_V = collect(d_V) + h_V⁻¹ = inv(h_V) + Eig = eigen(A) + @test Eig.values ≈ h_W + @test abs.(Eig.vectors*h_V⁻¹) ≈ I + d_A = CuArray(A) + d_W = CUSOLVER.Xgeev!('N','N', d_A) + h_W = collect(d_W) + @test Eig.values ≈ h_W + + A = rand(elty,m,m) + d_A = CuArray(A) + Eig = eigen(A) + d_eig = eigen(d_A) + @test Eig.values ≈ collect(d_eig.values) + h_V = collect(d_eig.vectors) + h_V⁻¹ = inv(h_V) + @test abs.(Eig.vectors* h_V⁻¹) ≈ I + end + @testset "syevd!" begin A = rand(elty,m,m) A += A' From 05ccc59dc0aa73d7f17eb64c2841673473199c69 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matteo=20Secl=C3=AC?= Date: Tue, 27 May 2025 15:07:10 +0200 Subject: [PATCH 04/11] Add tests for `eigvals` and `eigvecs` --- test/libraries/cusolver/dense.jl | 44 +++++++++++++++++++++++++++++++- 1 file changed, 43 insertions(+), 1 deletion(-) diff --git a/test/libraries/cusolver/dense.jl b/test/libraries/cusolver/dense.jl index 010b5306eb..6428cc8065 100644 --- a/test/libraries/cusolver/dense.jl +++ b/test/libraries/cusolver/dense.jl @@ -341,7 +341,20 @@ k = 1 @test Eig.values ≈ collect(d_eig.values) h_V = collect(d_eig.vectors) h_V⁻¹ = inv(h_V) - @test abs.(Eig.vectors* h_V⁻¹) ≈ I + @test abs.(Eig.vectors*h_V⁻¹) ≈ I + + A = rand(elty,m,m) + d_A = CuArray(A) + W = eigvals(A) + d_W = eigvals(d_A) + @test W ≈ collect(d_W) + + A = rand(elty,m,m) + d_A = CuArray(A) + V = eigvecs(A) + d_V = eigvecs(d_A) + V⁻¹ = inv(V) + @test abs.(collect(d_V)*V⁻¹) ≈ I end @testset "syevd!" begin @@ -398,6 +411,35 @@ k = 1 @test abs.(Eig.vectors'*h_V) ≈ I end + A = rand(elty,m,m) + A += A' + d_A = CuArray(A) + W = eigvals(LinearAlgebra.Hermitian(A)) + d_W = eigvals(d_A) + @test W ≈ collect(d_W) + d_W = eigvals(LinearAlgebra.Hermitian(d_A)) + @test W ≈ collect(d_W) + if elty <: Real + W = eigvals(LinearAlgebra.Symmetric(A)) + d_W = eigvals(LinearAlgebra.Symmetric(d_A)) + @test W ≈ collect(d_W) + end + + A = rand(elty,m,m) + A += A' + d_A = CuArray(A) + V = eigvecs(LinearAlgebra.Hermitian(A)) + V⁻¹ = inv(V) + d_V = eigvecs(d_A) + @test abs.(collect(d_V)*V⁻¹) ≈ I + d_V = eigvecs(LinearAlgebra.Hermitian(d_A)) + @test abs.(collect(d_V)*V⁻¹) ≈ I + if elty <: Real + V = eigvecs(LinearAlgebra.Symmetric(A)) + V⁻¹ = inv(V) + d_V = eigvecs(LinearAlgebra.Symmetric(d_A)) + @test abs.(collect(d_V)*V⁻¹) ≈ I + end end @testset "sygvd!" begin From 757cba39a3e8e20db7df697d9eb76fcc79b8d54e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matteo=20Secl=C3=AC?= Date: Tue, 27 May 2025 15:07:51 +0200 Subject: [PATCH 05/11] Typo --- lib/cusolver/linalg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/cusolver/linalg.jl b/lib/cusolver/linalg.jl index a5234bd499..b46c3caf39 100644 --- a/lib/cusolver/linalg.jl +++ b/lib/cusolver/linalg.jl @@ -169,7 +169,7 @@ function LinearAlgebra.eigvecs(A::Hermitian{T, <:CuMatrix}) where {T <: BlasComp return heevd!('V', 'U', A2)[2] end function LinearAlgebra.eigvecs(A::Hermitian{T, <:CuMatrix}) where {T <: BlasReal} - return eigvals(Symmetric(A)) + return eigvecs(Symmetric(A)) end function LinearAlgebra.eigvecs(A::CuMatrix{T}) where {T <: BlasReal} From 01f436c24ae7f8e4b57ab2534bc531259f58f0e0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matteo=20Secl=C3=AC?= Date: Tue, 27 May 2025 21:38:56 +0200 Subject: [PATCH 06/11] Fix most tests --- lib/cusolver/linalg.jl | 42 ++++++++++++++++++---------- test/libraries/cusolver/dense.jl | 47 ++++++++++++++++++-------------- 2 files changed, 54 insertions(+), 35 deletions(-) diff --git a/lib/cusolver/linalg.jl b/lib/cusolver/linalg.jl index b46c3caf39..336c9859c3 100644 --- a/lib/cusolver/linalg.jl +++ b/lib/cusolver/linalg.jl @@ -110,6 +110,20 @@ function Base.:\(F::Union{LinearAlgebra.LAPACKFactorizations{<:Any,<:CuArray}, return LinearAlgebra._cut_B(BB, 1:n) end +# Adapted from LinearAlgebra.sorteig!(). +# Warning: not very efficient, but works. +eigsortby(λ::Real) = λ +eigsortby(λ::Complex) = (real(λ),imag(λ)) +function sorteig!(λ::AbstractVector, X::AbstractMatrix, sortby::Union{Function,Nothing}=eigsortby) + if sortby !== nothing # && !issorted(λ, by=sortby) + p = sortperm(λ; by=sortby) + λ .= λ[p] # permute!(λ, p) + X .= X[:, p] # Base.permutecols!!(X, p) + end + return λ, X +end +sorteig!(λ::AbstractVector, sortby::Union{Function,Nothing}=eigsortby) = sortby === nothing ? λ : sort!(λ, by=sortby) + # eigen function LinearAlgebra.eigen(A::Symmetric{T,<:CuMatrix}) where {T<:BlasReal} @@ -127,23 +141,23 @@ end function LinearAlgebra.eigen(A::CuMatrix{T}) where {T<:BlasReal} A2 = copy(A) r = Xgeev!('N', 'V', A2) - return Eigen(r[1], r[3]) + return Eigen(sorteig!(r[1], r[3])...) end function LinearAlgebra.eigen(A::CuMatrix{T}) where {T<:BlasComplex} A2 = copy(A) r = Xgeev!('N', 'V', A2) - return Eigen(r[1], r[3]) + return Eigen(sorteig!(r[1], r[3])...) end # eigvals function LinearAlgebra.eigvals(A::Symmetric{T, <:CuMatrix}) where {T <: BlasReal} A2 = copy(A.data) - return syevd!('N', 'U', A2)[1] + return syevd!('N', 'U', A2) end function LinearAlgebra.eigvals(A::Hermitian{T, <:CuMatrix}) where {T <: BlasComplex} A2 = copy(A.data) - return heevd!('N', 'U', A2)[1] + return heevd!('N', 'U', A2) end function LinearAlgebra.eigvals(A::Hermitian{T, <:CuMatrix}) where {T <: BlasReal} return eigvals(Symmetric(A)) @@ -151,34 +165,34 @@ end function LinearAlgebra.eigvals(A::CuMatrix{T}) where {T <: BlasReal} A2 = copy(A) - return Xgeev!('N', 'N', A2)[1] + return sorteig!(Xgeev!('N', 'N', A2)[1]) end function LinearAlgebra.eigvals(A::CuMatrix{T}) where {T <: BlasComplex} A2 = copy(A) - return Xgeev!('N', 'N', A2)[1] + return sorteig!(Xgeev!('N', 'N', A2)[1]) end # eigvecs function LinearAlgebra.eigvecs(A::Symmetric{T, <:CuMatrix}) where {T <: BlasReal} - A2 = copy(A.data) - return syevd!('V', 'U', A2)[2] + E = eigen(A) + return E.vectors end function LinearAlgebra.eigvecs(A::Hermitian{T, <:CuMatrix}) where {T <: BlasComplex} - A2 = copy(A.data) - return heevd!('V', 'U', A2)[2] + E = eigen(A) + return E.vectors end function LinearAlgebra.eigvecs(A::Hermitian{T, <:CuMatrix}) where {T <: BlasReal} return eigvecs(Symmetric(A)) end function LinearAlgebra.eigvecs(A::CuMatrix{T}) where {T <: BlasReal} - A2 = copy(A) - return Xgeev!('N', 'V', A2)[3] + E = eigen(A) + return E.vectors end function LinearAlgebra.eigvecs(A::CuMatrix{T}) where {T <: BlasComplex} - A2 = copy(A) - return Xgeev!('N', 'V', A2)[3] + E = eigen(A) + return E.vectors end # factorizations diff --git a/test/libraries/cusolver/dense.jl b/test/libraries/cusolver/dense.jl index 6428cc8065..20939e5a3e 100644 --- a/test/libraries/cusolver/dense.jl +++ b/test/libraries/cusolver/dense.jl @@ -316,23 +316,27 @@ k = 1 end @testset "geev!" begin - A = rand(elty,m,m) - d_A = CuArray(A) + ## Note: we have Xgeev in dense_generic.jl, but no geev in dense.jl. + # A = rand(elty,m,m) + # d_A = CuArray(A) local d_W, d_V - d_W, d_V = CUSOLVER.Xgeev!('N','V', d_A) - d_W_b, d_V_b = LAPACK.geev!('N','V', CuArray(A)) - @test d_W ≈ d_W_b - @test d_V ≈ d_V_b - h_W = collect(d_W) - h_V = collect(d_V) - h_V⁻¹ = inv(h_V) - Eig = eigen(A) - @test Eig.values ≈ h_W - @test abs.(Eig.vectors*h_V⁻¹) ≈ I - d_A = CuArray(A) - d_W = CUSOLVER.Xgeev!('N','N', d_A) - h_W = collect(d_W) - @test Eig.values ≈ h_W + # d_W, _, d_V = CUSOLVER.Xgeev!('N','V', d_A) + # # d_W_b, _, d_V_b = LAPACK.geev!('N','V', CuArray(A)) + # # @test d_W ≈ d_W_b + # # @test d_V ≈ d_V_b + # W_b, _, V_b = LAPACK.geev!('N','V', A) + # @test collect(d_W) ≈ W_b + # @test collect(d_V) ≈ V_b + # h_W = collect(d_W) + # h_V = collect(d_V) + # h_V⁻¹ = inv(h_V) + # Eig = eigen(A) + # @test Eig.values ≈ h_W + # @test abs.(Eig.vectors*h_V⁻¹) ≈ I + # d_A = CuArray(A) + # d_W = CUSOLVER.Xgeev!('N','N', d_A) + # h_W = collect(d_W) + # @test Eig.values ≈ h_W A = rand(elty,m,m) d_A = CuArray(A) @@ -429,16 +433,17 @@ k = 1 A += A' d_A = CuArray(A) V = eigvecs(LinearAlgebra.Hermitian(A)) - V⁻¹ = inv(V) d_V = eigvecs(d_A) - @test abs.(collect(d_V)*V⁻¹) ≈ I + h_V = collect(d_V) + @test abs.(V'*h_V) ≈ I d_V = eigvecs(LinearAlgebra.Hermitian(d_A)) - @test abs.(collect(d_V)*V⁻¹) ≈ I + h_V = collect(d_V) + @test abs.(V'*h_V) ≈ I if elty <: Real V = eigvecs(LinearAlgebra.Symmetric(A)) - V⁻¹ = inv(V) d_V = eigvecs(LinearAlgebra.Symmetric(d_A)) - @test abs.(collect(d_V)*V⁻¹) ≈ I + h_V = collect(d_V) + @test abs.(V'*h_V) ≈ I end end From 804fe1e3af2279c17697df131fccf4893b2222f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matteo=20Secl=C3=AC?= Date: Wed, 28 May 2025 00:28:48 +0200 Subject: [PATCH 07/11] Fix typo --- test/libraries/cusolver/dense.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/libraries/cusolver/dense.jl b/test/libraries/cusolver/dense.jl index 20939e5a3e..d594ab1399 100644 --- a/test/libraries/cusolver/dense.jl +++ b/test/libraries/cusolver/dense.jl @@ -345,7 +345,7 @@ k = 1 @test Eig.values ≈ collect(d_eig.values) h_V = collect(d_eig.vectors) h_V⁻¹ = inv(h_V) - @test abs.(Eig.vectors*h_V⁻¹) ≈ I + @test abs.(h_V⁻¹*Eig.vectors) ≈ I A = rand(elty,m,m) d_A = CuArray(A) @@ -358,7 +358,7 @@ k = 1 V = eigvecs(A) d_V = eigvecs(d_A) V⁻¹ = inv(V) - @test abs.(collect(d_V)*V⁻¹) ≈ I + @test abs.(V⁻¹*collect(d_V)) ≈ I end @testset "syevd!" begin From ef3bd367d32cad963a2d4180a321be376966ead0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matteo=20Secl=C3=AC?= Date: Wed, 28 May 2025 01:36:24 +0200 Subject: [PATCH 08/11] Fix real `geev` eigenvectors --- lib/cusolver/linalg.jl | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/lib/cusolver/linalg.jl b/lib/cusolver/linalg.jl index 336c9859c3..c6a2f8da31 100644 --- a/lib/cusolver/linalg.jl +++ b/lib/cusolver/linalg.jl @@ -140,8 +140,22 @@ end function LinearAlgebra.eigen(A::CuMatrix{T}) where {T<:BlasReal} A2 = copy(A) - r = Xgeev!('N', 'V', A2) - return Eigen(sorteig!(r[1], r[3])...) + W, _, VR = Xgeev!('N', 'V', A2) + C = Complex{T} + U = CuMatrix{C}([1. 1.; im -im]) + VR = CuMatrix{C}(VR) + h_W = collect(W) + n = length(W) + j = 1 + while j <= n + if imag(h_W[j]) == 0 + j += 1 + else + VR[:,j:j+1] .= VR[:,j:j+1] * U + j += 2 + end + end + return Eigen(sorteig!(W, VR)...) end function LinearAlgebra.eigen(A::CuMatrix{T}) where {T<:BlasComplex} A2 = copy(A) From 1d9e95faccfc6f01efb9e55589e706cd6b30deb6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matteo=20Secl=C3=AC?= Date: Wed, 28 May 2025 11:27:26 +0200 Subject: [PATCH 09/11] Do eigenpairs sorting in tests only --- lib/cusolver/linalg.jl | 26 ++++++-------------------- test/libraries/cusolver/dense.jl | 22 ++++++++++++++++++++++ 2 files changed, 28 insertions(+), 20 deletions(-) diff --git a/lib/cusolver/linalg.jl b/lib/cusolver/linalg.jl index c6a2f8da31..72e2113f29 100644 --- a/lib/cusolver/linalg.jl +++ b/lib/cusolver/linalg.jl @@ -110,20 +110,6 @@ function Base.:\(F::Union{LinearAlgebra.LAPACKFactorizations{<:Any,<:CuArray}, return LinearAlgebra._cut_B(BB, 1:n) end -# Adapted from LinearAlgebra.sorteig!(). -# Warning: not very efficient, but works. -eigsortby(λ::Real) = λ -eigsortby(λ::Complex) = (real(λ),imag(λ)) -function sorteig!(λ::AbstractVector, X::AbstractMatrix, sortby::Union{Function,Nothing}=eigsortby) - if sortby !== nothing # && !issorted(λ, by=sortby) - p = sortperm(λ; by=sortby) - λ .= λ[p] # permute!(λ, p) - X .= X[:, p] # Base.permutecols!!(X, p) - end - return λ, X -end -sorteig!(λ::AbstractVector, sortby::Union{Function,Nothing}=eigsortby) = sortby === nothing ? λ : sort!(λ, by=sortby) - # eigen function LinearAlgebra.eigen(A::Symmetric{T,<:CuMatrix}) where {T<:BlasReal} @@ -142,7 +128,7 @@ function LinearAlgebra.eigen(A::CuMatrix{T}) where {T<:BlasReal} A2 = copy(A) W, _, VR = Xgeev!('N', 'V', A2) C = Complex{T} - U = CuMatrix{C}([1. 1.; im -im]) + U = CuMatrix{C}([1.0 1.0; im -im]) VR = CuMatrix{C}(VR) h_W = collect(W) n = length(W) @@ -151,16 +137,16 @@ function LinearAlgebra.eigen(A::CuMatrix{T}) where {T<:BlasReal} if imag(h_W[j]) == 0 j += 1 else - VR[:,j:j+1] .= VR[:,j:j+1] * U + VR[:, j:(j + 1)] .= VR[:, j:(j + 1)] * U j += 2 end end - return Eigen(sorteig!(W, VR)...) + return Eigen(W, VR) end function LinearAlgebra.eigen(A::CuMatrix{T}) where {T<:BlasComplex} A2 = copy(A) r = Xgeev!('N', 'V', A2) - return Eigen(sorteig!(r[1], r[3])...) + return Eigen(r[1], r[3]) end # eigvals @@ -179,11 +165,11 @@ end function LinearAlgebra.eigvals(A::CuMatrix{T}) where {T <: BlasReal} A2 = copy(A) - return sorteig!(Xgeev!('N', 'N', A2)[1]) + return Xgeev!('N', 'N', A2)[1] end function LinearAlgebra.eigvals(A::CuMatrix{T}) where {T <: BlasComplex} A2 = copy(A) - return sorteig!(Xgeev!('N', 'N', A2)[1]) + return Xgeev!('N', 'N', A2)[1] end # eigvecs diff --git a/test/libraries/cusolver/dense.jl b/test/libraries/cusolver/dense.jl index d594ab1399..a9e4415477 100644 --- a/test/libraries/cusolver/dense.jl +++ b/test/libraries/cusolver/dense.jl @@ -9,6 +9,20 @@ p = 5 l = 13 k = 1 +# Adapted from LinearAlgebra.sorteig!(). +# Warning: not very efficient, but works. +eigsortby(λ::Real) = λ +eigsortby(λ::Complex) = (real(λ), imag(λ)) +function sorteig!(λ::AbstractVector, X::AbstractMatrix, sortby::Union{Function, Nothing} = eigsortby) + if sortby !== nothing # && !issorted(λ, by=sortby) + p = sortperm(λ; by = sortby) + λ .= λ[p] # permute!(λ, p) + X .= X[:, p] # Base.permutecols!!(X, p) + end + return λ, X +end +sorteig!(λ::AbstractVector, sortby::Union{Function, Nothing} = eigsortby) = sortby === nothing ? λ : sort!(λ, by = sortby) + @testset "elty = $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64] @testset "gesv!" begin @testset "irs_precision = AUTO" begin @@ -342,6 +356,7 @@ k = 1 d_A = CuArray(A) Eig = eigen(A) d_eig = eigen(d_A) + sorteig!(d_eig.values, d_eig.vectors) @test Eig.values ≈ collect(d_eig.values) h_V = collect(d_eig.vectors) h_V⁻¹ = inv(h_V) @@ -351,12 +366,15 @@ k = 1 d_A = CuArray(A) W = eigvals(A) d_W = eigvals(d_A) + sorteig!(d_W) @test W ≈ collect(d_W) A = rand(elty,m,m) d_A = CuArray(A) V = eigvecs(A) + d_W = eigvals(d_A) d_V = eigvecs(d_A) + sorteig!(d_W, d_V) V⁻¹ = inv(V) @test abs.(V⁻¹*collect(d_V)) ≈ I end @@ -402,6 +420,7 @@ k = 1 d_A = CuArray(A) Eig = eigen(LinearAlgebra.Hermitian(A)) d_eig = eigen(d_A) + sorteig!(d_eig.values, d_eig.vectors) @test Eig.values ≈ collect(d_eig.values) d_eig = eigen(LinearAlgebra.Hermitian(d_A)) @test Eig.values ≈ collect(d_eig.values) @@ -420,6 +439,7 @@ k = 1 d_A = CuArray(A) W = eigvals(LinearAlgebra.Hermitian(A)) d_W = eigvals(d_A) + sorteig!(d_W) @test W ≈ collect(d_W) d_W = eigvals(LinearAlgebra.Hermitian(d_A)) @test W ≈ collect(d_W) @@ -433,7 +453,9 @@ k = 1 A += A' d_A = CuArray(A) V = eigvecs(LinearAlgebra.Hermitian(A)) + d_W = eigvals(d_A) d_V = eigvecs(d_A) + sorteig!(d_W, d_V) h_V = collect(d_V) @test abs.(V'*h_V) ≈ I d_V = eigvecs(LinearAlgebra.Hermitian(d_A)) From 8152b97b72337f37f57cfb60578e2dce11c0f5c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matteo=20Secl=C3=AC?= Date: Wed, 28 May 2025 12:01:01 +0200 Subject: [PATCH 10/11] Test `geev` for supported CUSOLVER versions only --- test/libraries/cusolver/dense.jl | 124 ++++++++++++++++--------------- 1 file changed, 66 insertions(+), 58 deletions(-) diff --git a/test/libraries/cusolver/dense.jl b/test/libraries/cusolver/dense.jl index a9e4415477..3ff1cb5e8d 100644 --- a/test/libraries/cusolver/dense.jl +++ b/test/libraries/cusolver/dense.jl @@ -329,54 +329,56 @@ sorteig!(λ::AbstractVector, sortby::Union{Function, Nothing} = eigsortby) = sor end end - @testset "geev!" begin - ## Note: we have Xgeev in dense_generic.jl, but no geev in dense.jl. - # A = rand(elty,m,m) - # d_A = CuArray(A) - local d_W, d_V - # d_W, _, d_V = CUSOLVER.Xgeev!('N','V', d_A) - # # d_W_b, _, d_V_b = LAPACK.geev!('N','V', CuArray(A)) - # # @test d_W ≈ d_W_b - # # @test d_V ≈ d_V_b - # W_b, _, V_b = LAPACK.geev!('N','V', A) - # @test collect(d_W) ≈ W_b - # @test collect(d_V) ≈ V_b - # h_W = collect(d_W) - # h_V = collect(d_V) - # h_V⁻¹ = inv(h_V) - # Eig = eigen(A) - # @test Eig.values ≈ h_W - # @test abs.(Eig.vectors*h_V⁻¹) ≈ I - # d_A = CuArray(A) - # d_W = CUSOLVER.Xgeev!('N','N', d_A) - # h_W = collect(d_W) - # @test Eig.values ≈ h_W - - A = rand(elty,m,m) - d_A = CuArray(A) - Eig = eigen(A) - d_eig = eigen(d_A) - sorteig!(d_eig.values, d_eig.vectors) - @test Eig.values ≈ collect(d_eig.values) - h_V = collect(d_eig.vectors) - h_V⁻¹ = inv(h_V) - @test abs.(h_V⁻¹*Eig.vectors) ≈ I - - A = rand(elty,m,m) - d_A = CuArray(A) - W = eigvals(A) - d_W = eigvals(d_A) - sorteig!(d_W) - @test W ≈ collect(d_W) + if CUSOLVER.version() >= v"11.7.1" + @testset "geev!" begin + ## Note: we have Xgeev in dense_generic.jl, but no geev in dense.jl. + # A = rand(elty,m,m) + # d_A = CuArray(A) + local d_W, d_V + # d_W, _, d_V = CUSOLVER.Xgeev!('N','V', d_A) + # # d_W_b, _, d_V_b = LAPACK.geev!('N','V', CuArray(A)) + # # @test d_W ≈ d_W_b + # # @test d_V ≈ d_V_b + # W_b, _, V_b = LAPACK.geev!('N','V', A) + # @test collect(d_W) ≈ W_b + # @test collect(d_V) ≈ V_b + # h_W = collect(d_W) + # h_V = collect(d_V) + # h_V⁻¹ = inv(h_V) + # Eig = eigen(A) + # @test Eig.values ≈ h_W + # @test abs.(Eig.vectors*h_V⁻¹) ≈ I + # d_A = CuArray(A) + # d_W = CUSOLVER.Xgeev!('N','N', d_A) + # h_W = collect(d_W) + # @test Eig.values ≈ h_W + + A = rand(elty,m,m) + d_A = CuArray(A) + Eig = eigen(A) + d_eig = eigen(d_A) + sorteig!(d_eig.values, d_eig.vectors) + @test Eig.values ≈ collect(d_eig.values) + h_V = collect(d_eig.vectors) + h_V⁻¹ = inv(h_V) + @test abs.(h_V⁻¹*Eig.vectors) ≈ I + + A = rand(elty,m,m) + d_A = CuArray(A) + W = eigvals(A) + d_W = eigvals(d_A) + sorteig!(d_W) + @test W ≈ collect(d_W) - A = rand(elty,m,m) - d_A = CuArray(A) - V = eigvecs(A) - d_W = eigvals(d_A) - d_V = eigvecs(d_A) - sorteig!(d_W, d_V) - V⁻¹ = inv(V) - @test abs.(V⁻¹*collect(d_V)) ≈ I + A = rand(elty,m,m) + d_A = CuArray(A) + V = eigvecs(A) + d_W = eigvals(d_A) + d_V = eigvecs(d_A) + sorteig!(d_W, d_V) + V⁻¹ = inv(V) + @test abs.(V⁻¹*collect(d_V)) ≈ I + end end @testset "syevd!" begin @@ -419,9 +421,11 @@ sorteig!(λ::AbstractVector, sortby::Union{Function, Nothing} = eigsortby) = sor A += A' d_A = CuArray(A) Eig = eigen(LinearAlgebra.Hermitian(A)) - d_eig = eigen(d_A) - sorteig!(d_eig.values, d_eig.vectors) - @test Eig.values ≈ collect(d_eig.values) + if CUSOLVER.version() >= v"11.7.1" + d_eig = eigen(d_A) + sorteig!(d_eig.values, d_eig.vectors) + @test Eig.values ≈ collect(d_eig.values) + end d_eig = eigen(LinearAlgebra.Hermitian(d_A)) @test Eig.values ≈ collect(d_eig.values) h_V = collect(d_eig.vectors) @@ -438,9 +442,11 @@ sorteig!(λ::AbstractVector, sortby::Union{Function, Nothing} = eigsortby) = sor A += A' d_A = CuArray(A) W = eigvals(LinearAlgebra.Hermitian(A)) - d_W = eigvals(d_A) - sorteig!(d_W) - @test W ≈ collect(d_W) + if CUSOLVER.version() >= v"11.7.1" + d_W = eigvals(d_A) + sorteig!(d_W) + @test W ≈ collect(d_W) + end d_W = eigvals(LinearAlgebra.Hermitian(d_A)) @test W ≈ collect(d_W) if elty <: Real @@ -453,11 +459,13 @@ sorteig!(λ::AbstractVector, sortby::Union{Function, Nothing} = eigsortby) = sor A += A' d_A = CuArray(A) V = eigvecs(LinearAlgebra.Hermitian(A)) - d_W = eigvals(d_A) - d_V = eigvecs(d_A) - sorteig!(d_W, d_V) - h_V = collect(d_V) - @test abs.(V'*h_V) ≈ I + if CUSOLVER.version() >= v"11.7.1" + d_W = eigvals(d_A) + d_V = eigvecs(d_A) + sorteig!(d_W, d_V) + h_V = collect(d_V) + @test abs.(V'*h_V) ≈ I + end d_V = eigvecs(LinearAlgebra.Hermitian(d_A)) h_V = collect(d_V) @test abs.(V'*h_V) ≈ I From 56be9b36b6de26e9b3b66131bffc1aa6ff6edfbc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matteo=20Secl=C3=AC?= Date: Wed, 28 May 2025 14:16:24 +0200 Subject: [PATCH 11/11] Remove irrelevant tests --- test/libraries/cusolver/dense.jl | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/test/libraries/cusolver/dense.jl b/test/libraries/cusolver/dense.jl index 3ff1cb5e8d..ecc9078edc 100644 --- a/test/libraries/cusolver/dense.jl +++ b/test/libraries/cusolver/dense.jl @@ -331,27 +331,7 @@ sorteig!(λ::AbstractVector, sortby::Union{Function, Nothing} = eigsortby) = sor if CUSOLVER.version() >= v"11.7.1" @testset "geev!" begin - ## Note: we have Xgeev in dense_generic.jl, but no geev in dense.jl. - # A = rand(elty,m,m) - # d_A = CuArray(A) local d_W, d_V - # d_W, _, d_V = CUSOLVER.Xgeev!('N','V', d_A) - # # d_W_b, _, d_V_b = LAPACK.geev!('N','V', CuArray(A)) - # # @test d_W ≈ d_W_b - # # @test d_V ≈ d_V_b - # W_b, _, V_b = LAPACK.geev!('N','V', A) - # @test collect(d_W) ≈ W_b - # @test collect(d_V) ≈ V_b - # h_W = collect(d_W) - # h_V = collect(d_V) - # h_V⁻¹ = inv(h_V) - # Eig = eigen(A) - # @test Eig.values ≈ h_W - # @test abs.(Eig.vectors*h_V⁻¹) ≈ I - # d_A = CuArray(A) - # d_W = CUSOLVER.Xgeev!('N','N', d_A) - # h_W = collect(d_W) - # @test Eig.values ≈ h_W A = rand(elty,m,m) d_A = CuArray(A)