# Computing eigen-decomposition of several matrices in parallel in C++

I am writing a program in C++ in which I am trying to reduce the run-time by computing eigen-decomposition of several matrices in parallel. This might be a programming question but since many Physics people might use it, I am asking the question here. This might be a simple example of how to perform computational tasks in parallel. Can anybody tell through a simple example, how to do this ? (Assume I have a multi-core processor )

• It would help to hear (i) how big the matrices are, (ii) how many matrices you have, (iii) what you have been using to compute eigendecompositions on a single core, etc. Commented Jul 12, 2013 at 14:22
• @WolfgangBangerth : matrices are of the order 10^5 * 10^5. I am using the armadillo library to compute eigendecompositions on a single core. Commented Jul 19, 2013 at 16:04

## 2 Answers

You would use something like OpenMP to exploit thread level parallelism. An example in C++ would be like

#pragma omp parallel for
for(int i = 0; i < num_matrices; i++){
DoEigendecomposition(matrix[i]);
}


The pragma will automatically cause each loop iteration to be performed in parallel (assuming you compile with the right flags; e.g. -fopenmp for GCC). Within DoEigendecomposition, you would allocate the necessary memory for temporary work arrays (must be distinct for each thread, so don't make them static or shared), and call the proper routine to do the eigen-decomposition as you would in the single-core case (most likely in Lapack for dense matrices, or some Krylov-like method for sparse). You can also use pthreads, but that is much lower level, and you have to manually create, synchronize, and destroy the threads.

Note: You need to make sure that whatever code does the eigendecomposition is thread-safe (does not use static variables in C, no COMMON or SAVE blocks in Fortran). Lapack's non-symmetric eigensolvers (double and zomplex precision) typically are not thread-safe due to the way the library is usually compiled. For details, see this forum post. Disclaimer: I am the original bug reporter.

• Note that COMMON is outdated (only F77) and SAVE was "fixed" for the Fortran 90 updates. See here for more details on SAVE. Commented Jul 13, 2013 at 0:28

Parallel solution for eigenvalues and eigenvectors of a given matrix can be performed with the package SLEPc available online (which is built on top of a large package PETSc). There are usage examples in the SLEPc distribution.