I am trying to solve a generalized eigenvalue problem using Arpack, right now the code is using LAPACK but that's too slow, we only need a few eigenvalues and the matrices are sparse so using Arpack should be the way to go.
Before I start working with the original code I decided to test a simple case using scipy wrapper for Arpack (eigs) but the results that I am getting are wrong and change every time the code runs.
Minimum working example:
import numpy as np
from scipy.linalg import eig
from scipy.sparse.linalg import eigs
n = 8
A = np.diag(np.arange(1,n+1,1.0))
B = np.eye(n) # We want symmetric but a non-diagonal B. eigs gives correct answer for B=np.eye(n)
B[0][n-1] = 2
B[n-1][0] = 2
evals,_ = eigs(A,k=3,M=B,which='LM')
print("The eigenvalues obtained by eigs (uses Arpack)")
print(evals)
print("Correct eigenvalues using eig (uses Lapack):")
evals_l,_ = eig(A,b=B)
print(evals_l)
```