Dear Jed,<br><br>For multithread support in PETSc, my question is whether KSP and/or PC work when Vec and Mat use multithread mode.<br><br>Thanks,<br>Yujie<br><br><div class="gmail_quote">On Sun, Feb 26, 2012 at 6:39 PM, Jed Brown <span dir="ltr"><<a href="mailto:jedbrown@mcs.anl.gov">jedbrown@mcs.anl.gov</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0pt 0pt 0pt 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="gmail_quote"><div class="im">On Sun, Feb 26, 2012 at 04:07, Gerard Gorman <span dir="ltr"><<a href="mailto:g.gorman@imperial.ac.uk" target="_blank">g.gorman@imperial.ac.uk</a>></span> wrote:<br>
</div><div class="im"><blockquote class="gmail_quote" style="margin:0pt 0pt 0pt 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div>> Did you post a repository yet? I'd like to have a look at the code.<br>
<br>
</div>It's on Barry's favourite collaborative software development site of<br>
course ;-)<br>
<br>
<a href="https://bitbucket.org/wence/petsc-dev-omp/overview" target="_blank">https://bitbucket.org/wence/petsc-dev-omp/overview</a><br></blockquote></div></div><br><div>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).</div>
<div><br></div><div>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.</div>
<div><br></div><div>It's really not acceptable to insert unguarded</div><div><br></div><div>#pragma omp ...</div><div><br></div><div>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</div>
<div><br></div><div>#define PetscPragmatize(x) _Pragma(#x)</div><div>#if defined(PETSC_USE_OPENMP)</div><div># define PetscPragmaOMP(x) PetscPragmatize(omp x)</div><div>#else</div><div># define PetscPragmaOMP(x)</div><div>
#endif</div><div><br></div><div>then use</div><div><br></div><div>PetscPragmaOMP(parallel for ...)</div><div><br></div><div>We should probably use a variant for object-based threading</div><div><br></div><div>#define PetscPragmaOMPObject(obj,x) PetscPragmaOMP(x num_threads((obj)->nthreads))</div>
<div><br></div><div>In the case of multiple objects, I think you usually want the object being modified to control the number of threads.</div><div><br></div><div>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</div>
<div><br></div><div>PragmaOMP(parallel) {</div><div> int nthreads = omp_get_num_threads();</div><div> int tnum = omp_get_thread_num();</div><div> int start,end;</div><div> // compute start and end</div><div> for (int i=start; i<end; i++) {</div>
<div> // the work</div><div> }</div><div>}</div><div><br></div><div>We could perhaps capture some of this common logic in a macro:</div><div><br></div><div>#define VecOMPParallelBegin(X,args) do { \</div><div> PragmaOMPObject(X,parallel args) { \</div>
<div> PetscInt _start, _end; \</div><div> VecOMPGetThreadLocalPart(X,&_start,&_end); \</div><div> { do {} while(0)</div><div><br></div><div>#define VecOMPParallelEnd() }}} while(0)</div><div><br></div><div>VecOMPParallelBegin(X, shared/private ...);</div>
<div>{</div><div> PetscInt i;</div><div> for (i=_start; i<_end; i++) {</div><div> // the work</div><div> }</div><div>}</div><div>VecOMPParallelEnd();</div><div><br></div><div>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. 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.)</div>
</blockquote></div><br>