OpenMP nested parallelism bug in BLIS

[jkd2016]

There appears to be an threading issue with BLIS compiled with OpenMP and run inside a parallel nested loop.

It crashes the Elk code with a segmentation fault.

I can't reproduce the seg. fault with a small example, but the following program never terminates on our AMD machine:

program test
implicit none
integer, parameter :: n=20
integer i
real(8), allocatable :: a(:,:),b(:,:),c(:,:)

!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(a,b,c)
!$OMP DO
do i=1,100
  allocate(a(n,n),b(n,n),c(n,n))
  a=1.d0
  b=1.d0
  call dgemm('N','N',n,n,n,1.d0,a,n,b,n,0.d0,c,n)
  deallocate(a,b,c)
end do
!$OMP END DO
!$OMP END PARALLEL

end program

The above code has to be compiled with

gfortran -fopenmp test.f90 lapack.a libblis.a

Both the number of OpenMP and BLIS threads have to be larger than 1:

export OMP_NUM_THREADS=2
export BLIS_NUM_THREADS=2

If either OMP_NUM_THREADS=1 or BLIS_NUM_THREADS=1 then the code runs fine.

If BLIS is compiled with pthreads:

./configure -t pthreads zen

then the code also runs fine.

The official AMD release (https://developer.amd.com/amd-cpu-libraries/blas-library/) was compiled with OpenMP and therefore also has the bug.

[devinamatthews]

@fgvanzee So far I've been able to trace this as far as frame/3/bli_l3_packm.c where cntl_mem_p->pblk.buf_sys = 0x1, so almost definitely a memory pool problem. I'll keep digging.

[devinamatthews]

Interestingly the program runs just fine in valgrind, which is quite annoying.

[jeffhammond]

@devinamatthews Doesn't Valgrind pad allocations to detect overruns and thus may inadvertently fix the bug if the buffer overrun is by a small amount? I've found that Valgrind missed off-by-1 errors before.

[devinamatthews]

NVM, I recompiled and that went away and I got the behavior described by @jkd2016. Setting OMP_NESTED=TRUE makes it go away, so I think the fix is to use omp_get_max_threads() to limit how many threads are used (expected) in the OpenMP thread decorator.

[jeffhammond]

I'm not looking at the spec right now, but I'm not sure omp_get_max_threads() solves a problem caused by nested parallel regions. This function only tells you the maximum number of threads in a given parallel region GCC docs.

[devinamatthews]

Wouldn't it have to return 1 if nesting (at this level) is disabled, though? I get there is always omp_get_nested() but some implementations have a max nesting level/max number of total threads etc.

[devinamatthews]

Or, we could test omp_num_threads() after the fact and adjust accordingly.

[devinamatthews]

@jkd2016 will you try the nested-omp-patch branch and verify that it fixes your problem?

[devinamatthews]

Well, it doesn't quite work yet: I forgot that the number of threads is also saved in the rntm_t struct. But, the threads are partitioned amongst the various loops there so not easy to just change willy nilly. I will make the patch only work for getting a number of threads exactly equal to what we wanted, or one.

[devinamatthews]

It should work now.

[fgvanzee]

@devinamatthews I'm having a hard time understanding your patch. Would you mind walking me through what it does?

[devinamatthews]

Basically it checks if the actual number of threads spawned is equal to the number of threads requested. If not, then it reinitialized the number of threads (expressed by and in the n_threads, g_comm, and rntm variables). I specifically check if the new number of threads is one because if it is not, then it would be somewhat more complicated to reinitialize the rntm srtructure properly.

[devinamatthews]

The problem was that BLIS expected to get N threads, but the omp parallel region only spawned one (because nesting was disabled). Then, any barrier will get stuck while that one thread waits on the N-1 other nonexistent ones.