You might want to try the SelfAdjointEigenSolver class
in the Eigen C++ class library
I did some numerical experiments with a 30x30 SPD matrix with a condition
number of 1000 constructed
according to the procedure described by Neumaier here
Generating Symmetric Positive Definite Matrices using indices
I used a Windows 7 computer with a 2.8 GHz AMD Phenom II X4 830 CPU.
The compiler is VS 2013.
For the experiments with Lapack ssyev/dsyev, I used OpenBlas
In previous experiments I have found the BLAS operations in this library to
have very good performance. I suspect (but did not verify) that their
implementation of ssyev/dsyev simply uses the standard Fortran code but that
the performance of that routine depends significantly on the quality of the
underlying BLAS functions.
Eigen::SelfAdjointEigenSolver was significantly faster than ssyev/dsyev in this
case-- ~0.25 ms per call compared with ~2 ms per call.
I saw only very small reduction of the time when I switched from double to
single precision. But I noticed large errors in the eigenvalues computed.
This was true for both Eigen and Lapack.