CODEWITHSUNDEEP
c++openmplapack
SirGuy
SirGuy

Reputation: 10770

Program slower with OpenMP due to LAPACK

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?

Upvotes: 1

Views: 1130

Answers (1)

Stephen Canon
Stephen Canon

Reputation: 106167

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.

Upvotes: 4

Related Questions