Python scipy eigh(Arpack) giving wrong eigenvalues for generalized eigenvalue problem

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)
```