I am trying to diagonalize some dense, ill-conditioned matrices. In machine precision, results are inaccurate (returning negative eigenvalues, eigenvectors do not have the expected symmetries). I switched over to Mathematica's Eigensystem function to take advantage of arbitrary precision, but computations are extremely slow. I am open to any number of solutions. Are there packages/algorithms that are well suited to ill-conditioned problems? I am not an expert on preconditioning, so I am not sure how much this could help. Otherwise, all I can think of are parallelized arbitrary precision eigenvalue solvers, but I am not familiar with anything beyond Mathematica, MATLAB and C++.
To give some background on the problem, the matrices are large, but not huge (4096x4096 to 32768x32768 at the most). They are real, symmetric, and the eigenvalues are bounded between 0 and 1 (exclusive), with many eigenvalues being very close to 0 and none close to 1. The matrix is essentially a convolution operator. I do not need to diagonalize all of my matrices, but the larger I can go, the better. I have access to computing clusters with many processors and distributed computing capabilities.