Skip to content

Expand eigen() and add eig[vals,vecs]() #2787

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 11 commits into
base: master
Choose a base branch
from
Open

Expand eigen() and add eig[vals,vecs]() #2787

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
4becc5a
Expand eigen() and add eig[vals,vecs]()
matteosecli May 26, 2025
291ea99
Add returns
matteosecli May 26, 2025
11bf501
Add `geev` tests for `eigen`
matteosecli May 27, 2025
05ccc59
Add tests for `eigvals` and `eigvecs`
matteosecli May 27, 2025
757cba3
Typo
matteosecli May 27, 2025
01f436c
Fix most tests
matteosecli May 27, 2025
804fe1e
Fix typo
matteosecli May 27, 2025
ef3bd36
Fix real `geev` eigenvectors
matteosecli May 27, 2025
1d9e95f
Do eigenpairs sorting in tests only
matteosecli May 28, 2025
8152b97
Test `geev` for supported CUSOLVER versions only
matteosecli May 28, 2025
56be9b3
Remove irrelevant tests
matteosecli May 28, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
73 changes: 67 additions & 6 deletions lib/cusolver/linalg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,29 +110,90 @@ 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)
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)
issymmetric(A) ? Eigen(syevd!('V', 'U', A2)...) : error("GPU eigensolver supports only Hermitian or Symmetric matrices.")
W, _, VR = Xgeev!('N', 'V', A2)
C = Complex{T}
U = CuMatrix{C}([1.0 1.0; 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(W, VR)
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)
return Eigen(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)
end
function LinearAlgebra.eigvals(A::Hermitian{T, <:CuMatrix}) where {T <: BlasComplex}
A2 = copy(A.data)
return heevd!('N', 'U', A2)
end
function LinearAlgebra.eigvals(A::Hermitian{T, <:CuMatrix}) where {T <: BlasReal}
return eigvals(Symmetric(A))
end

function LinearAlgebra.eigvals(A::CuMatrix{T}) where {T <: BlasReal}
A2 = copy(A)
return Xgeev!('N', 'N', A2)[1]
end
function LinearAlgebra.eigvals(A::CuMatrix{T}) where {T <: BlasComplex}
A2 = copy(A)
return Xgeev!('N', 'N', A2)[1]
end

# eigvecs

function LinearAlgebra.eigvecs(A::Symmetric{T, <:CuMatrix}) where {T <: BlasReal}
E = eigen(A)
return E.vectors
end
function LinearAlgebra.eigvecs(A::Hermitian{T, <:CuMatrix}) where {T <: BlasComplex}
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}
E = eigen(A)
return E.vectors
end
function LinearAlgebra.eigvecs(A::CuMatrix{T}) where {T <: BlasComplex}
E = eigen(A)
return E.vectors
end

# factorizations

Expand Down
90 changes: 88 additions & 2 deletions test/libraries/cusolver/dense.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -315,6 +329,38 @@ k = 1
end
end

if CUSOLVER.version() >= v"11.7.1"
@testset "geev!" begin
local d_W, d_V

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
end
end

@testset "syevd!" begin
A = rand(elty,m,m)
A += A'
Expand Down Expand Up @@ -355,8 +401,11 @@ k = 1
A += A'
d_A = CuArray(A)
Eig = eigen(LinearAlgebra.Hermitian(A))
d_eig = eigen(d_A)
@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)
Expand All @@ -369,6 +418,43 @@ 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))
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
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))
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
if elty <: Real
V = eigvecs(LinearAlgebra.Symmetric(A))
d_V = eigvecs(LinearAlgebra.Symmetric(d_A))
h_V = collect(d_V)
@test abs.(V'*h_V) ≈ I
end
end

@testset "sygvd!" begin
Expand Down