I was trying to numerically compute eigenvalues and eigenvectors of a matrix in Python with Numpy. The problem was that the inverse similarity transformation of the eigenvalue matrix does not accurately reproduce the original matrix in the sense that the element-wise error was at the same order of the entries of the original matrix; instead, for other more "well-behave" matrices, I found that the error was of the order of, 1e-15, the floating point precision. All elements of the original matrix is of the same order. Concluding that the matrix is non-diagonalizable, I tried instead to compute the Jordan form of the same matrix in Mathematica, but instead what I got was a diagonal matrix with the same "inaccurate eigenvalues" as its diagonal, and exactly the same set of eigenvectors.
Both linear algebra routines of Numpy and Mathematica give accurate results when applied to ordinary matrices. I am looking for a numerical routine or treatment that I can try to apply to this problematic matrix to accurately diagonalize it or accurately compute its Jordan form.