Multithreading symmetric eigenvalue problem

[bluetakamaka]

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.

[rodney]

There are no global variables used in LAPACK, so if all of your variables are indeed local, it should work correctly. However, on some platforms, you may need to add an extra FORTRAN compiler flag to produce a thread-safe library. What type of computer and compiler are you using, and what are the compiler flags in your make.inc? Also, are you using the latest release of LAPACK, 3.4.1?

Rodney

[bluetakamaka]

Thank you for the quick reply. I just talked to our system administrator inquiring the answers to your questions:
version of LAPACK: 3.4.1
type of compiler: gcc 4.5.1
type of system: Intel Xeon CPU 64bit, operating system openSUSE 11.4, 64bit

compiler flags in make.inc (taken from the make.inc.example file):
SHELL = /bin/sh
FORTRAN = gfortran
OPTS = -O2
DRVOPTS = $(OPTS)
NOOPT = -O0
LOADER = gfortran
LOADOPTS =
TIMER = INT_ETIME
CC = gcc
CFLAGS = -O3
ARCH = ar
ARCHFLAGS= cr
RANLIB = ranlib
XBLASLIB =
BLASLIB = ../../librefblas.a
LAPACKLIB = liblapack.a
TMGLIB = libtmglib.a
LAPACKELIB = liblapacke.a

What would be an extra FORTRAN compiler flag that one has to set to produce a thread-safe library?

Joerg

[rodney]

Try adding the -frecursive option to OPTS and NOOPT.

Rodney

[bluetakamaka]

Thanks for the suggestion. Our system administrator compiled LAPACK with these additional flags, however, the same problem still persists.

Joerg

[CyLith]

You should verify that the flags were actually applied since the compilation process can be tricky. The best way to do this would be to write a simple test program that calls dsyevx in more than one thread at the same time, repeatedly, and check that all results are correct and consistent. I encountered a similarly frustrating bug (now in the errata #61) due to Lapack's use of large stack arrays (their statement that Lapack is thread safe is only a half truth). The eigenvector computation routines call the random number generators dlarnv and dlaruv, both of which allocate stack arrays (though not large ones). If your test code cannot reproduce the problem, it might be somewhere else in your code where you generate the matrices (I mistakenly thought I encountered a Lapack bug once when in reality I was accidentally calling the routine in multiple threads on the same input data (which is bad since the input matrix is modified)).