Problem summary

I want to calculate the smallest eigenvalue (algebraic value) of a matrix. The matrix comes from an op4 file that I read using the pyNastran library. Following these instructions, I am trying to calculate the smallest eigevanlue using the scipy.sparse.linalg.eigsh function with shift-invert mode. In order to check that the computation is correct, I compare the result of eigsh with the result of numpy.linalg.eigvals. What I observe is quite puzzling: if I simply apply eigsh to the matrix, the calculated eigenvalue is wrong, if I save the matrix to a csv file and then I load it back into a numpy array the eigenvalue is correct. What is even more baffling is that numpy.array_equal returns True when I compare the two matrices. How is it possible that eigsh returns two different results for the same matrix?

Code

from pyNastran.op4.op4 import read_op4
from scipy.sparse.linalg import eigsh
import numpy as np
op4 = read_op4('kllrh.op4')
matrix_name = 'KLLRH'
kllrh_matrix = op4[matrix_name][1][-1]
reference_max_eigenvalue = np.max(np.linalg.eigvals(kllrh_matrix))
reference_min_eigenvalue = np.min(np.linalg.eigvals(kllrh_matrix))
eigsh_min_eigenvalue = eigsh(kllrh_matrix, 1, sigma=0, which='LM', return_eigenvectors=False)
np.savetxt('kllrh.csv', kllrh_matrix, delimiter=',')
kllrh_matrix_reloaded = np.loadtxt('kllrh.csv', delimiter=",")
reloaded_eigsh_min_eigenvalue = eigsh(kllrh_matrix_reloaded, 1, sigma=0, which='LM', return_eigenvectors=False)
print(reference_min_eigenvalue)
print(eigsh_min_eigenvalue)
print(reloaded_eigsh_min_eigenvalue)
print(np.array_equal(kllrh_matrix, kllrh_matrix_reloaded))
print(type(kllrh_matrix))
print(type(kllrh_matrix_reloaded))
print(reference_max_eigenvalue)

This returns the following:

-0.0028387385
[0.05363945]
[-0.00283876]
True
<class 'numpy.ndarray'>
<class 'numpy.ndarray'>
6502976000.0

Please find the kllrh.op4 file here.

Although np.array_equal says that kllrh_matrix and kllrh_matrix_reloaded are equal, they are of different dtypes (float32 vs float64).

If you do

kllrh_matrix = op4[matrix_name][1][-1].astype('float64')

everything is correct:

-0.0028387384680708
[-0.00283876]
[-0.00283876]
True
<class 'numpy.ndarray'>
<class 'numpy.ndarray'>
6502976138.583383

As a better alternative, you can specify the precision when reading the op4 file:

op4 = read_op4('kllrh.op4', precision='double')

which will give you slightly other results:

0.004395871268066287
[0.00439589]
[0.00439589]
True
<class 'numpy.ndarray'>
<class 'numpy.ndarray'>
6502976180.582107

The default precision comes from the op4 file: in the example the matrix type is 1 (fourth number on first line of the file) which results in single precision (float32).