<div class="gmail_quote">On Tue, Feb 28, 2012 at 10:41, Lawrence Mitchell <span dir="ltr"><<a href="mailto:lawrence.mitchell@ed.ac.uk">lawrence.mitchell@ed.ac.uk</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div id=":6cj">It's possible that if one asks the gcc folk nicely enough that the overhead might go to zero.  Alternately, if it turns out that overheads are unavoidable, and we've gone down the macro-ification route anyway (see below) we could explicitly make a runtime decision based on numthreads in PETSc.  So that:<br>

<br>
PetscOMPParallelFor(..., numthreads)<br>
for ( ... )<br>
    ;<br>
<br>
would turn into<br>
if ( numthreads == 1 )<br>
    for ( ... )<br>
        ;<br>
else<br>
#pragma omp parallel for<br>
    for ( ... )<br>
        ;<br>
<br>
The branch should hopefully be cheaper than the thread sync overhead.</div></blockquote><div><br></div><div>The branch is trivially cheap, the growth in binary size (having a threaded and non-threaded) is possibly an issue. </div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div id=":6cj"><div class="im"></div>
As far as attaching the threading information to the object goes, it would be nice to do so staying within the sequential implementation (rather than copying all the code as is done for pthreads right now). I'd envisage something like this:<br>

<br>
PetscPragmaOMPObject(vec, parallel)<br>
{<br>
   PetscPragmaOMP(for...)<br>
   for ( i = PetscVecLowerBound(vec); i < PetscVecUpperBound(vec); i++ )<br>
       ...;<br>
}<br>
<br>
Where one has:<br>
<br>
#if defined(PETSC_HAVE_OPENMP)<br>
#define PetscVecLowerBound(vec) compute_lower_bound<br>
#define PetscVecUpperBound(vec) compute_upper_bound<br>
#else<br>
#define PetscVecLowerBound(vec) 0<br>
#define PetscVecUpperBound(vec) vec->map->n<br>
#endif<br>
<br>
Or is this too ugly for words?<br></div></blockquote><div><br></div><div>The reason for my magic was simply to avoid calling omp_get_num_threads() and omp_get_thread_num() many times to make this determination. Note that OpenMP implementations create functions (run by the threads) that capture the local state and call these functions to find their work unit. Thus the generated code is quite similar to PETSc's pthreads versions.</div>
<div> </div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div id=":6cj"><div class="im"></div>
So to take a single, representative, kernel we move from something like:<br>
<br>
PetscErrorCode VecConjugate_Seq(Vec xin)<br>
{<br>
  PetscScalar    *x;<br>
  PetscInt       n = xin->map->n;<br>
  PetscInt       i;<br>
  PetscErrorCode ierr;<br>
<br>
  PetscFunctionBegin;<br>
  ierr = VecGetArray(xin,&x);CHKERRQ(<u></u>ierr);<br>
#pragma omp parallel for schedule(static) default(none) private(i) shared(n, x)<br>
  for (i=0; i<n; i++) x[i] = PetscConj(x[i]);<br>
  ierr = VecRestoreArray(xin,&x);<u></u>CHKERRQ(ierr);<br>
  PetscFunctionReturn(0);<br>
}<br>
<br>
To something like this:<br>
<br>
PetscErrorCode VecConjugate_Seq(Vec xin)<br>
{<br>
    PetscScalar *x;<br>
    PetscInt i;<br>
    PetscErrorCode ierr;<br>
    PetscFunctionBegin;<br>
    ierr = VecGetArray(xin,&x);CHKERRQ(<u></u>ierr);<br>
    PetscPragmaOMPObject(xin, parallel private(i) shared(xin, x))<br>
    {<br>
        for (i=PetscVecLowerBound(xin); i<PetscVecUpperBound(xin); i++) x[i] = PetscConj(x[i]);<br>
    }<br>
    ierr = VecRestoreArray(xin,&x);<u></u>CHKERRQ(ierr);<br>
    PetscFunctionReturn(0);<br>
}<br>
<br>
I prefer inlining the loop bounds into macro/function calls at point of use, rather than relying on magic names from a macro definition, but I'm not completely wedded to the idea.</div></blockquote></div><br><div>Not a big deal, we could optimize either variant.</div>