Jed Brown emailed the following on 27/02/12 00:39:

On Sun, Feb 26, 2012 at 04:07, Gerard Gorman <[email protected] <mailto:[email protected]>> wrote:

> Did you post a repository yet? I'd like to have a look at the code.

It's on Barry's favourite collaborative software development site of course ;-)

https://bitbucket.org/wence/petsc-dev-omp/overview

I looked through the code and I'm concerned that all the OpenMP code is inlined into vec/impls/seq and mat/impls/aij/seq with, as far as I can tell, no way to use OpenMP for some objects in a simulation, but not others. I think that all the pragmas should have num_threads(x->nthreads) clauses. We can compute the correct number of threads based on sizes when memory is allocated (or specified through command line options, inherited from related objects, etc).

The num_threads(x->nthreads) is worth investigating. However, in the benchmarks I have done so far it would seem that the only two sensible values for nthreads would be 1 or the total number of threads. When you set num_threads(1) it seems that OpenMP (for some implementations at least) is sensible enough not to do anything silly that would introduce scheduling/synchronisation overheads. Once you use more than one thread then you incur the scheduling/synchronisation overheads. As you increase the number of threads you may for small arrays see parallel efficiency decreasing, but I have not seen the actual time increasing. If this is a general result then it might not be a big deal that OpenMP is inlined into vec/impls/seq and mat/impls/aij/seq so long as the num_threads was used as suggested. For determining where the cut off array size for using all or one thread - as interesting option would be to determine this at run time, i.e. learn at run time what the appropriate is.

I don't think we can get away from OpenMP schedule overhead (several hundred cycles) even for those objects that we choose not to use threads for, but (at least with my gcc tests), that overhead is only about a third of the cost of actually starting a parallel region.

I had a quick go at trying to get some sensible benchmarks for this but I was getting too much system noise. I am particularly interested in seeing if the overhead goes to zero if num_threads(1) is used. The next step is to take a look at the EPCC OpenMP microbenchmarks to see if they have tied down these issues.

It's really not acceptable to insert unguarded

#pragma omp ...

into the code because this will generate tons of warnings or errors with compilers that don't know about OpenMP. It would be better to test for _Pragma and use

I'm surprised by this. I not aware of any compiler that doesn't have OpenMP support - and then you do not actually enable OpenMP compilers generally just ignore the pragma. Do you know of any compiler that does not have OpenMP support which will complain?

#define PetscPragmatize(x) _Pragma(#x) #if defined(PETSC_USE_OPENMP) # define PetscPragmaOMP(x) PetscPragmatize(omp x) #else # define PetscPragmaOMP(x) #endif

then use

PetscPragmaOMP(parallel for ...)

We should probably use a variant for object-based threading

#define PetscPragmaOMPObject(obj,x) PetscPragmaOMP(x num_threads((obj)->nthreads))

This may not be flexible enough. You frequently want to have a parallel region, and then have multiple omp for's within that one region.

In the case of multiple objects, I think you usually want the object being modified to control the number of threads.

I take this point.

In many cases, I would prefer more control over the partition of the loop. For example, in many cases, I'd be willing to tolerate a slight computational imbalance between threads in exchange for working exclusively within my page. Note that the arithmetic to compute such things is orders of magnitude less expensive than the schedule/distribution to threads. I don't know how to do that except to

PragmaOMP(parallel) { int nthreads = omp_get_num_threads(); int tnum = omp_get_thread_num(); int start,end; // compute start and end for (int i=start; i<end; i++) { // the work } }

We could perhaps capture some of this common logic in a macro:

#define VecOMPParallelBegin(X,args) do { \ PragmaOMPObject(X,parallel args) { \ PetscInt _start, _end; \ VecOMPGetThreadLocalPart(X,&_start,&_end); \ { do {} while(0)

#define VecOMPParallelEnd() }}} while(0)

VecOMPParallelBegin(X, shared/private ...); { PetscInt i; for (i=_start; i<_end; i++) { // the work } } VecOMPParallelEnd();

That should reasonably give us complete run-time control of the number of parallel threads per object and their distribution, within the constraints of contiguous thread partition.

I think what you describe is close to Fig 3 of this paper written by your neighbours: http://greg.bronevetsky.com/papers/2008IWOMP.pdf However, before making the implementation more complex, it would be good to benchmark the current approach and use a tool like likwid to measure the NUMA traffic so we can get a good handle on the costs.

That also leaves open the possibility of using libnuma to query and migrate pages. (For example, a short vector that needs to be accessed from multiple NUMA nodes might intentionally be faulted with pages spread apart even though other vectors of similar size might be accessed from within one NUMA nodes and thus not use threads at all. (One 4 KiB page is only 512 doubles, but if the memory is local to a single NUMA node, we wouldn't use threads until the vector length was 4 to 8 times larger.)

Well this is where the implementation details get richer and there are many options - they also become less portable. For example, what does all this mean for the sparc64 processors which are UMA. Not to mention Intel MIC which also supports OpenMP. I guess I am cautious about getting too bogged down with very invasive optimisations until we have benchmarked the basic approach which in a wide range of use cases will achieve good thread/page locality as illustrated previously. Cheers Gerard