On Thu, Nov 24, 2011 at 7:28 PM, Matthew Knepley <[email protected]> wrote:

On Thu, Nov 24, 2011 at 6:37 PM, Barry Smith <[email protected]> wrote:

...

On Nov 24, 2011, at 4:41 PM, Matthew Knepley wrote:

...

On Thu, Nov 24, 2011 at 4:09 PM, Barry Smith <[email protected]> wrote:

Jed,

Let's stop arguing about whether MPI is or is not a good base for the next generation of HPC software but instead start a new conversation on what API (implemented on top of or not on top of MPI/pthreads etc etc) we want to build PETSc on to scale PETSc up to millions of cores with large NUMA nodes and GPU like accelerators.

What do you want in the API?

Let's start with the "lowest" level, or at least the smallest. I think the only sane way to program for portable performance here is using CUDA-type vectorization. This SIMT style is explained well here http://www.yosefk.com/blog/simd-simt-smt-parallelism-in-nvidia-gpus.html I think this is much easier and more portable than the intrinsics for Intel, and more performant and less error prone than threads. I think you can show that it will accomplish anything we want to do. OpenCL seems to have capitulated on this point. Do we agree here?

What syntax do you suggest for writing the code that is "vectorized"? What tools exist, could exist, for mapping from that syntax to what is needed by the various compilers/hardware?

Some history. Brook is the basis for CUDA, but like any good foundation, almost everything the creator thought was important was thrown away to make something usable. Brook is a streaming language, much like Thrust. When problems fit this paradigm, it is fantastic. However, CUDA is not a streaming language. The programmer decides exactly what memory to put where, so what is there?

Both CUDA and OpenCL inherit the process grid. Of course, MPI already has this (rank), so what is different? All threads in a vector have access to shared memory. I guess you could get the same thing in MPI if you had an idea of a "neighborhood" which had shared memory. This is exactly how OpenCL handles it. You specify a "vector length" (our neighborhood size) and the compiler tries it ass off to vectorize the code. In fact, you could probably manage everything I want with subcommunicators and a runtime code generator.

Synopsis of what I said before to elicit comment: 1) I think the only thing we can learn from Brook, CUDA, OpenCL is that you identify threads by a grid ID. 2) Things like BLAS are so easy that you can move up to the streaming model, but this does not work for - FD and FEM residual evaluation (Jed has an FD example with Aron, SNES ex52 is my FEM example) - FD and FEM Jacobian evaluation 3) If you look at ex52 I do a "thread transposition" meaning threads start working on different areas of memory which looks like a transpose on a 2D grid. I can do this using shared memory for the vector group. The API is very simple. Give grid indices to the thread, and its done in CUDA and OpenCL essentially the same way. Matt

For daxpy() the syntax doesn't really matter, anything will do. For

...

other kernels: maxpy, sparse matrix vector product, triangular solves, P*A*Pt, mesh operations, sorts, indirect access .... the choice of syntax likely matters a great deal. We should test the syntax out on a wide range of kernels. For example look at VecMAXPY_kernel vs VecMAXPY_Seq vs VecMAXPY_VecCUSPMAXPY4 and the three delegators VecMAXPY_SeqPThread, VecMAXPY_MPI, and VecMAXPY_SeqCUSP;

For things as easy as BLAS, you can go all the way to the streaming-type kernel (Thrust, PyCUDA, etc). I think this is basically solved. I am more interested in kernels like FD residual (ask Jed), FEM residual, and FEM Jacobian. I promise to finish this paper by mid-December.

...

How does data layout relate to the vectorization you are going to do on that data and vis-versa?

That is the crux. Vectorization is about execution layout (the CUDA thread grid). Somehow we must also layout memory and match them up. This is all of programming (which is why I initially like the stuff from Victor).

Matt

...

Barry

...

Matt

Barry

-- What most experimenters take for granted before they begin their

experiments is infinitely more interesting than any results to which their experiments lead.

...

-- Norbert Wiener

-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener

-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener