Generalized eigensolver for Hermitian matrices fails

[abussy]

Calling the generalized eigensolver eigen(Hermitian(A), Hermitian(B)) fails when both A and B are Hermitian{T, <:ROCArray}. The error message reads:

ERROR: LoadError: PosDefException: matrix is not positive definite; Factorization failed.

Here is a reproducer:

using LinearAlgebra
using AMDGPU

### On the CPU
A = randn(50, 50)
A = A * A'
B = randn(50, 50)
B = B * B'

ref_vals, ref_vecs = eigen(A, B)
ref_hvals, ref_hvecs = eigen(Hermitian(A), Hermitian(B))

### On the GPU
A_gpu = ROCArray(A)
B_gpu = ROCArray(B)

gpu_vals, gpu_vecs = eigen(A_gpu, B_gpu) # This works, and is correct
@assert maximum(abs, ref_vals .- Array(gpu_vals)) < 1.0e-14
@assert maximum(abs, ref_vecs .- Array(gpu_vecs)) < 1.0e-14

eigen(Hermitian(A_gpu), Hermitian(B_gpu)) # This fails

Interestingly enough, eigen(A, B) without the explicit Hermitian wrapper works, despite methods(eigen) showing no specific method for eigen(A::ROCArray{T, 2}, B::ROCArray{T, 2}). I guess we are lucky that the default LinearAlgebra path works on AMD GPUs.

For completeness, I have a workaround:

function LinearAlgebra.eigen(A::Hermitian{T, <:ROCArray}, B::Hermitian{T, <:ROCArray}) where {T}     
    C = cholesky(B)                                                                                  
    D = C.L \ A / C.L'                                                                               
    F = eigen(Hermitian(D))                                                                          
    return (; values = F.values, vectors = C.L' \ F.vectors)                                         
end