Arpack.jl icon indicating copy to clipboard operation
Arpack.jl copied to clipboard

Generalized eigenvalue problem with symmetric matrices: complex solution!?

Open PetrKryslUCSD opened this issue 3 years ago • 2 comments

Arpack v0.5.3: The generalized eigenvalue problem with symmetric matrices should return real eigenvalues and real eigenvectors. The eigs function fails to do that:

using LinearAlgebra
using SparseArrays 
using Arpack
neigvs = 2
K = sparse([1, 4, 5, 2, 3, 1, 4, 8, 1, 5, 8, 6, 7, 4, 5, 8], [1, 1, 1, 2, 3, 4, 4, 4, 5, 5, 5, 6, 7, 8, 8, 8], [2.37230e+05, -2.80193e+05, 1.18615e+05, 8.82496e+05, 6.80939e+08, -2.80193e+05, 8.82496e+05, 2.80193e+05, 1.18615e+05, 4.74459e+05, 1.18615e+05, 1.00348e+05, 4.74459e+05, 2.80193e+05, 1.18615e+05, 2.37230e+05], 8, 8)                                      
M = sparse([1, 4, 5, 2, 3, 1, 4, 8, 1, 5, 8, 6, 7, 4, 5, 8], [1, 1, 1, 2, 3, 4, 4, 4, 5, 5, 5, 6, 7, 8, 8, 8], [3.16529e-01, 8.08509e-01, -2.37158e-01, 1.52896e+01, 1.37166e+01, 8.08509e-01, 1.52896e+01, -8.08509e-01, -2.37158e-01, 6.33057e-01, -2.37158e-01, 4.77869e-03, 6.33057e-01, -8.08509e-01, -2.37158e-01, 3.16529e-01], 8, 8) 
evals, evecs, nconv = eigs(Symmetric(K), Symmetric(M); nev=neigvs, which=:SM);

# First  we should check that the requested eigenvalues actually converged:
@show nconv == neigvs

fs = sqrt.(evals) / (2 * pi);
println("Approximate frequencies: $fs [Hz]")
println("Target frequencies: $([1.66013e+01, 3.76332e+01]) [Hz]")

gives

nconv == neigvs = true                                                                                                                                                                 
Approximate frequencies: ComplexF64[1.66645e+01 + 0.00000e+00im, 3.82365e+01 + 0.00000e+00im] [Hz]                                                                                     
Target frequencies: [1.66013e+01, 3.76332e+01] [Hz]   

This works as expected with Arpack v0.4.

PetrKryslUCSD avatar Jan 13 '22 15:01 PetrKryslUCSD

I found that setting explicittransform = :none in the argument list fixes this issue. I believe the explicit transform was going to be set to default to :none?

PetrKryslUCSD avatar Jan 13 '22 16:01 PetrKryslUCSD

https://github.com/JuliaLinearAlgebra/Arpack.jl/pull/120?

PetrKryslUCSD avatar Jan 13 '22 16:01 PetrKryslUCSD