Program slower with OpenMP due to LAPACK

1
votes

I have a program that needs to find the eigenvalues and eigenvectors of 3x3 matrices millions of times. I have just switched to using LAPACK's zheev (yes they are hermitian matrices) for this, and the program runs in about 1m20s for a particular case. I have parallelized my algorithm with OpenMP (as we were doing before) and suddenly my program runs in about 9m. I comment out the call to zheev and my program runs in 9s.

I have looked around online and found (as I understand it) that you can compile your BLAS library to use OpenMP, but I don't think that is the issue here.

Unfortunately this code is from my work, I don't have lapack installed in the default location and I don't know what compiler options were used when it was compiled. This also makes it difficult for me to compile a minimum testing program to demonstrate the issue.

Any ideas on what the issue could be?

EDIT:

I just found out that with OpenMP zheev is failing, which probably ties in to it running slower. I have read that some routines in LAPACK are not thread safe (or they have thread safe variants), how can I find out if zheev is calling one of those routines and can I change that?

1
c++openmplapack
How long does the program take to run serially with the call to zheev commented out? Also: what platform are you on, and which distributions of LAPACK and BLAS are you using? - Stephen Canon
@StephenCanon good point, it takes about 23s, OpenSuSe 12 with LAPACK 3.2.1 and as far as I can tell the BLAS that came with it - SirGuy
LAPACK 3.2.1 should be thread safe, except that you may need to call dlamch and slamch before parallel execution starts. - Stephen Canon
@StephenCanon Do those functions get called if I don't call them explicitly? - SirGuy
@StephenCanon You may be pleased to know that my remembering your comment about dlamch and slamch from two years ago was the only thing that helped me figure out a recent problem I was having with just that. - SirGuy

1 Answers

4
votes

Leaving aside your OpenMP issues for the moment, if your code is performance-sensitive, you may not want to use LAPACK to find the eigenvalues and eigenvectors of a 3x3 matrix; LAPACK is targeted at "large" problems. More importantly, for the specific case of matrices of dimension smaller than 5, you can directly compute the eigenvalues, so you can use a simpler algorithm than is used for general matrices (which necessarily requires iteration).

Recall that the characteristic polynomial of a 3x3 matrix is a cubic polynomial, which means that you can directly compute its roots (which are the eigenvalues). Once you know the eigenvalues, you can directly solve (A - lambda * I)x = 0 for each eigenvalue lambda to get the corresponding eigenvector.