In the following minimal example, I am trying to create a few Eigen::Matrix in a for loop parallelized by OMP. Each matrix is contained within the loop, so there is no data sharing or race conditions among the threads. The code works perfectly when number of threads equals one, otherwise I get a segmentation fault. What is mind boggling that I get seg faults for matrix of size 600x600, but not for e.g., 599x599 or 601x601 or 1000x1000. Any help is appreciated. Thanks :)
#include <iostream>
#define EIGEN_DONT_ALIGN_STATICALLY
#define EIGEN_STACK_ALLOCATION_LIMIT 0
#include <Eigen/Core>
#define SIZE 600
#define THREADS 2
int main(int argc, char *argv[]) {
// The following code always works for THREADS=1
// When THREADS!=1, there is a seg fault if SIZE=600.
// There is no seg fault when THREADS!=1 and SIZE=599 or SIZE=601
#pragma omp parallel for num_threads(THREADS)
for(int n=0; n<5; ++n){
Eigen::Matrix<double,SIZE,SIZE> mat = Eigen::Matrix<double,SIZE,SIZE>::Zero();
}
return 0;
}
The problem is coming from fixed sizes and the limited stack size. Indeed, Matrix<double,SIZE,SIZE> is allocated on the stack and since SIZE is quite big, it should allocate 2.7 MiB on most systems while the stack size is usually set to few MiB (set by default 2 MiB on most Linux platforms). You should use Dynamic matrix sizes instead.
See the documentation for more information:
The limitation of using fixed sizes, of course, is that this is only possible when you know the sizes at compile time. Also, for large enough sizes, say for sizes greater than (roughly) 32, the performance benefit of using fixed sizes becomes negligible. Worse, trying to create a very large matrix using fixed sizes inside a function could result in a stack overflow, since Eigen will try to allocate the array automatically as a local variable, and this is normally done on the stack.