I am using the routine dsyevx for computing eigenpairs of symmetric matrices. The routine is called through the interface of the finite element library deal.ii
When running the program on a single thread everything works very nicely. However, when using two or more threads the eigenvektors of some matrices are not computed correctly - unlike the eigenvalues, which are always correct. I realize that according to the LAPACK 3.3 release notes "all routines in LAPACK-3.3 are now thread-safe.", which hopefully also applies to LAPACK 3.4, which is what I am using. Nevertheless, I was wondering if anyone has experienced a similar problem.
As I see it the room for error in my code is rather limited. My matrices are local variables on different threads. Right before computing the eigenpairs I output each matrix to a file. After the computation I do the same for the eigenvalues and eigenvectors. All variables are local. When I check the results by Matlab, I see that the eigenvectors for some matrices are not computed correctly.
I realize that one source of error could be the interface of LAPACK with deall.ii. On the other hand, the code of this interface is rather short and clear, and the main developer of deal.ii agrees that the source of error is not there.