diff --git a/lib/cusolver/linalg.jl b/lib/cusolver/linalg.jl index 0f2d9fac49..72e2113f29 100644 --- a/lib/cusolver/linalg.jl +++ b/lib/cusolver/linalg.jl @@ -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 diff --git a/test/libraries/cusolver/dense.jl b/test/libraries/cusolver/dense.jl index 4df323226c..ecc9078edc 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 @@ -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' @@ -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) @@ -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