search

JuliaLinearAlgebra / Arpack.jl

Julia Wrappers for the arpack-ng Fortran library
https://arpack.julialinearalgebra.org/stable
MIT License
69 stars 29 forks source link

Generalized eigenvalue problem with symmetric matrices: complex solution!? #140

Open PetrKryslUCSD opened 2 years ago

PetrKryslUCSD commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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