November 2011 - petsc-dev - mailman3.cels.anl.gov

Re: [petsc-dev] programming model for PETSc
by Jed Brown 26 Nov '11

26 Nov '11

On Fri, Nov 25, 2011 at 16:48, Matthew Knepley <knepley(a)gmail.com> wrote: > 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 > I think these are also probably too simple. Discontinuous Galerkin with overlapped flux computations and interior integration would be a somewhat better model problem. Nonlinear Gauss-Seidel in a multigrid context would be another. > > 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. > As is, this seems to assume a flat memory model and the memory access only appears in how the kernel uses threadIdx to determine what memory to operate on. If we could say something about this up-front, then the library could schedule tasks relative to memory and perhaps handle some updates for distributed memory. Can we have a way to specify the required memory access before launching the kernels?

1 0

Re: [petsc-dev] Generality of VecScatter
by Jed Brown 26 Nov '11

26 Nov '11

On Sat, Nov 26, 2011 at 10:35, Mark F. Adams <mark.adams(a)columbia.edu>wrote: > Thats what I'm asking. I assumed that there is some buffering going on > here but that is a later question. > I assumed that you would either (a) go in discrete rounds where every process does as much as it can or (b) do some predetermined amount of work, initiate the update, do another work unit, complete the last update and initiate a new one, do another work unit, etc. > Remember that my first primitive was "broadcast", which is essentially > global-to-local. The map is as a local-to-global map and the "get" is > initiated at the local space. > > > I was confused, and perhaps still am, the broadcast is the one sided pull > of ghost data to the local space. So a broadcast is passive for the > broadcaster whereas it implied an active process to me. The reduce is a > push of ghost data with the same (only) map to the global vertex. > > The way I'm writing the algorithm there is no need for a reduce, I think, > but this is a detail, this algorithm broadcasts the selected vertex status > (and deleted which I forgot in my psuedo-code) and then it infers the local > deleted w/o having to have it sent (reduced). > Sure, either way works. If the data volume for updates was big enough to be a performance limitation (instead of the bottleneck being latency), then we could convert to a message-based model initiated by the owner. It's easy to extract this when we want it, but it's simpler to skip it. > > So I'm thinking a good way to implement this MIS with your model is after > each vertex loop, you do what you call a "broadcast" of all new status. > > Now I want to see exactly how that is implemented, here is a guess: > As part of the distributed graph, you have (owner rank, index) for every ghosted point. As local setup, the communication API makes a local (where in my array the ghost points reside) and remote (where in the owner's array the points reside) MPI_Datatype. Store these in an array of length num_neighbor_ranks. > > MPI_Win_fence > loop over all of your boundary vertices that are not "done" yet, and do an > MPI_Get (I assume that these requests are buffered in someway) > Not over boundary vertices, it loops over neighbor each neighbor rank and calls one MPI_Get() to fetch the status on all vertices owned by that rank that are ghosted on this rank. > MPI_Win_fence. > The communication need only be finished once you need the result, so you can do more work before calling the second MPI_Win_fence().

1 0

Re: [petsc-dev] Generality of VecScatter
by Mark F. Adams 26 Nov '11

26 Nov '11

On Nov 25, 2011, at 9:13 PM, Jed Brown wrote: > On Fri, Nov 25, 2011 at 20:06, Mark F. Adams <mark.adams(a)columbia.edu> wrote: >> MIS is straightforward. To define the parallel graph, let's say that ghosting processes have (owner rank, index) of neighbors. An owned point is only a candidate for selection if all neighbors have the same rank or smaller (according to Rule 2 of your paper). >> >> while (candidates remain): >> forall v in owned_candidates with neighbor ranks smaller than p: > > This is not quite right, and I want to see how complex or cumbersome this one sided model will get when you get all the details in. > > You could say loop: > > for all owned_candidates v > for all neighbors q of v > if q is selected > delete v > skip v > else if q.proc > myproc and q is not deleted > skip v > endif > end for > if not skip v then > select v > delete all q > broadcast v status ???? > > Why broadcast the single value immediately? Thats what I'm asking. I assumed that there is some buffering going on here but that is a later question. > I figured that we would just store the status until we are done with what we can do locally, then we'll do reduce-and-broadcast and move on to the next round. > > endif > end for > reduce all status to owners ??? > > Yeah, we reduce the status from the local spaces (which may have marked a point as deleted), then broadcast. The communication graph is stored on the local side. > > > I thought in your one-sided model you could not broadcast, or push data. You can only pull data and only want to pull data. But you have a broadcast??? > > Remember that my first primitive was "broadcast", which is essentially global-to-local. The map is as a local-to-global map and the "get" is initiated at the local space. I was confused, and perhaps still am, the broadcast is the one sided pull of ghost data to the local space. So a broadcast is passive for the broadcaster whereas it implied an active process to me. The reduce is a push of ghost data with the same (only) map to the global vertex. The way I'm writing the algorithm there is no need for a reduce, I think, but this is a detail, this algorithm broadcasts the selected vertex status (and deleted which I forgot in my psuedo-code) and then it infers the local deleted w/o having to have it sent (reduced). So I'm thinking a good way to implement this MIS with your model is after each vertex loop, you do what you call a "broadcast" of all new status. Now I want to see exactly how that is implemented, here is a guess: MPI_Win_fence loop over all of your boundary vertices that are not "done" yet, and do an MPI_Get (I assume that these requests are buffered in someway) MPI_Win_fence. This is not complete but am I on the right track? Mark

1 0

Re: [petsc-dev] Generality of VecScatter
by Jed Brown 25 Nov '11

25 Nov '11

On Fri, Nov 25, 2011 at 20:06, Mark F. Adams <mark.adams(a)columbia.edu>wrote: > MIS is straightforward. To define the parallel graph, let's say that > ghosting processes have (owner rank, index) of neighbors. An owned point is > only a candidate for selection if all neighbors have the same rank or > smaller (according to Rule 2 of your paper). > > while (candidates remain): > forall v in owned_candidates with neighbor ranks smaller than p: > > > This is not quite right, and I want to see how complex or cumbersome this > one sided model will get when you get all the details in. > > You could say loop: > > for all owned_candidates v > for all neighbors q of v > if q is selected > delete v > skip v > else if q.proc > myproc and q is not deleted > skip v > endif > end for > if not skip v then > select v > delete all q > broadcast v status ???? > Why broadcast the single value immediately? I figured that we would just store the status until we are done with what we can do locally, then we'll do reduce-and-broadcast and move on to the next round. > endif > end for > reduce all status to owners ??? > Yeah, we reduce the status from the local spaces (which may have marked a point as deleted), then broadcast. The communication graph is stored on the local side. > > I thought in your one-sided model you could not broadcast, or push data. > You can only pull data and only want to pull data. But you have a > broadcast??? > Remember that my first primitive was "broadcast", which is essentially global-to-local. The map is as a local-to-global map and the "get" is initiated at the local space.

1 0

Re: [petsc-dev] Generality of VecScatter
by Mark F. Adams 25 Nov '11

25 Nov '11

On Nov 25, 2011, at 4:22 PM, Jed Brown wrote: > On Fri, Nov 25, 2011 at 13:07, Mark F. Adams <mark.adams(a)columbia.edu> wrote: > OK, thats a good start. You do need two communication steps, 1) get offsets, 2) get data (right?). > > Yes, to send variable-sized data, you need two steps. For constant-sized data, it is just one. > > You avoid packing the send data but you need to unpack the receive data, converting the raw graph data into the reduced ghost data. Not simple but writing from one side only and with only one-to-one communication is attractive. > > That depends a bit on what representation you want for the ghost data. You never "unpack" in a by-rank way since it gets put into the correct part of the domain according to where an MPI_Datatype says to put it. (This was part of the local setup when defining the MPI_Get() calls, but those are internal to the library.) > > If you are increasing the level of ghosting, you probably want to number the new ghost nodes that you acquired and make the local data structure point directly at them instead of at the (rank, offset) at which they are owned. This is a global-to-local translation of the connectivity, so you can do it with a hash table or sorting. I think you can also skip sorting by pushing the new indices back through the global space. > > That is, I effectively have two local-to-global maps, both pointing at the same global space: L1toG and L2toG. I want to translate some local information defined on L1 into L2. I can do this by fetching global indices and inverting the local part of the map, or (with the assistance of some memory provided by the global owner) I can do a gather L2 -> G* (G* has a slot for each process sharing via this map) and then get those indices back in L1 using scatter G* -> L1. I think the hash is usually better, but with a fast network and/or highly vectorized hardware, there might be an advantage to pushing the indices L2 -> G* -> L1 without the need for hashing or sorting. > > > So keep going. AMR is good but I'm not sure I'm up to specifying a complete unstructured AMR algorithm at the moment. My other example of AMG coarse grid process aggregation is another.... > > Or perhaps a parallel maximal independent set algorithm, preferably my algorithm which is implemented in GAMG and documented in an old paper. > > MIS is straightforward. To define the parallel graph, let's say that ghosting processes have (owner rank, index) of neighbors. An owned point is only a candidate for selection if all neighbors have the same rank or smaller (according to Rule 2 of your paper). > > while (candidates remain): > forall v in owned_candidates with neighbor ranks smaller than p: This is not quite right, and I want to see how complex or cumbersome this one sided model will get when you get all the details in. You could say loop: for all owned_candidates v for all neighbors q of v if q is selected delete v skip v else if q.proc > myproc and q is not deleted skip v endif end for if not skip v then select v delete all q broadcast v status ???? endif end for reduce all status to owners ??? I thought in your one-sided model you could not broadcast, or push data. You can only pull data and only want to pull data. But you have a broadcast??? Mark > mark v as selected, mark neighbors(v) as deleted (in local space) > reduce status to owner > broadcast status to ghosters > > In any round, multiple procs might have marked the vertex for deletion, but only if the vertex has a neighbor with larger owner rank, in which case it could not have been selected.

1 0

Re: [petsc-dev] programming model for PETSc
by Matthew Knepley 25 Nov '11

25 Nov '11

On Thu, Nov 24, 2011 at 7:28 PM, Matthew Knepley <knepley(a)gmail.com> wrote: > On Thu, Nov 24, 2011 at 6:37 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > >> >> On Nov 24, 2011, at 4:41 PM, Matthew Knepley wrote: >> >> > On Thu, Nov 24, 2011 at 4:09 PM, Barry Smith <bsmith(a)mcs.anl.gov> >> 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

1 0

Re: [petsc-dev] Generality of VecScatter
by Jed Brown 25 Nov '11

25 Nov '11

On Fri, Nov 25, 2011 at 13:07, Mark F. Adams <mark.adams(a)columbia.edu>wrote: > OK, thats a good start. You do need two communication steps, 1) get > offsets, 2) get data (right?). > Yes, to send variable-sized data, you need two steps. For constant-sized data, it is just one. > You avoid packing the send data but you need to unpack the receive data, > converting the raw graph data into the reduced ghost data. Not simple but > writing from one side only and with only one-to-one communication is > attractive. > That depends a bit on what representation you want for the ghost data. You never "unpack" in a by-rank way since it gets put into the correct part of the domain according to where an MPI_Datatype says to put it. (This was part of the local setup when defining the MPI_Get() calls, but those are internal to the library.) If you are increasing the level of ghosting, you probably want to number the new ghost nodes that you acquired and make the local data structure point directly at them instead of at the (rank, offset) at which they are owned. This is a global-to-local translation of the connectivity, so you can do it with a hash table or sorting. I think you can also skip sorting by pushing the new indices back through the global space. That is, I effectively have two local-to-global maps, both pointing at the same global space: L1toG and L2toG. I want to translate some local information defined on L1 into L2. I can do this by fetching global indices and inverting the local part of the map, or (with the assistance of some memory provided by the global owner) I can do a gather L2 -> G* (G* has a slot for each process sharing via this map) and then get those indices back in L1 using scatter G* -> L1. I think the hash is usually better, but with a fast network and/or highly vectorized hardware, there might be an advantage to pushing the indices L2 -> G* -> L1 without the need for hashing or sorting. > So keep going. AMR is good but I'm not sure I'm up to specifying a > complete unstructured AMR algorithm at the moment. My other example of AMG > coarse grid process aggregation is another.... > > Or perhaps a parallel maximal independent set algorithm, preferably my > algorithm which is implemented in GAMG and documented in an old paper. > MIS is straightforward. To define the parallel graph, let's say that ghosting processes have (owner rank, index) of neighbors. An owned point is only a candidate for selection if all neighbors have the same rank or smaller (according to Rule 2 of your paper). while (candidates remain): forall v in owned_candidates with neighbor ranks smaller than p: mark v as selected, mark neighbors(v) as deleted (in local space) reduce status to owner broadcast status to ghosters In any round, multiple procs might have marked the vertex for deletion, but only if the vertex has a neighbor with larger owner rank, in which case it could not have been selected.

1 0

Re: [petsc-dev] Generality of VecScatter
by Mark F. Adams 25 Nov '11

25 Nov '11

On Nov 25, 2011, at 1:15 PM, Jed Brown wrote: > On Fri, Nov 25, 2011 at 12:00, Mark F. Adams <mark.adams(a)columbia.edu> wrote: >> With my model, the owner never needs to be informed that some procs are ghosting another level, nor be aware what those ghosted nodes are. It may place some data in an array for the "pointwise broadcast" of connectivity, but it doesn't need semantic knowledge that that information will used to increase the ghosting. Similarly, any process can stop ghosting a point without informing the owner in any way. > > I don't see how you can do that, you must have a different data model then me. We may need a white board for this but if I want to get an extra layer of ghosts I need to have the remote process tell me what they are. I have a distributed graph so I need to be told who my new ghosts are. > > Suppose for simplicity that the remote process stores it's connectivity as a directed graph. In the local data structure, each vertex has on offset into a connectivity array that lists the other vertices that it is connected to. This is like CSR storage without the weights. We will communicate this CSR storage in-place, without packing and without knowledge of how many remote processes accessed it. > > In order to ghost the original points, the ghosters needed to know (owner rank, index). (This is my "native" representation for ghosting.) That means that I can fetch offset and row length directly from the owner's "row starts" array. With that, I can fetch the rows ("column indices") directly from the owner's storage. Underneath my (thin) API, we'll just be using MPI_Get() for these things. > OK, thats a good start. You do need two communication steps, 1) get offsets, 2) get data (right?). You avoid packing the send data but you need to unpack the receive data, converting the raw graph data into the reduced ghost data. Not simple but writing from one side only and with only one-to-one communication is attractive. So keep going. AMR is good but I'm not sure I'm up to specifying a complete unstructured AMR algorithm at the moment. My other example of AMG coarse grid process aggregation is another.... Or perhaps a parallel maximal independent set algorithm, preferably my algorithm which is implemented in GAMG and documented in an old paper. Mark > The overall semantics are collective, in that the owner needs to provide a send buffer and call MPI_Win_fence(), but it only provides one send buffer (no packing), and each process gets what it needs out of that buffer (by creating a suitable MPI_Datatype for MPI_Get()). The owner does not know how many procs accessed the data or what they accessed. > > I was trying to avoid specifying an algorithm, but complete repartitioning is inherently complex. I was thinking of a diffusive kind of thing with nearest neighbors. > > Sure, we can do that sort of thing. There are a variety of ways to "claim" an interface vertex, reducing with MAXLOC is one way.

1 0

Re: [petsc-dev] Generality of VecScatter
by Jed Brown 25 Nov '11

25 Nov '11

On Fri, Nov 25, 2011 at 12:00, Mark F. Adams <mark.adams(a)columbia.edu>wrote: > With my model, the owner never needs to be informed that some procs are > ghosting another level, nor be aware what those ghosted nodes are. It may > place some data in an array for the "pointwise broadcast" of connectivity, > but it doesn't need semantic knowledge that that information will used to > increase the ghosting. Similarly, any process can stop ghosting a point > without informing the owner in any way. > > > I don't see how you can do that, you must have a different data model then > me. We may need a white board for this but if I want to get an extra layer > of ghosts I need to have the remote process tell me what they are. I have > a distributed graph so I need to be told who my new ghosts are. > Suppose for simplicity that the remote process stores it's connectivity as a directed graph. In the local data structure, each vertex has on offset into a connectivity array that lists the other vertices that it is connected to. This is like CSR storage without the weights. We will communicate this CSR storage in-place, without packing and without knowledge of how many remote processes accessed it. In order to ghost the original points, the ghosters needed to know (owner rank, index). (This is my "native" representation for ghosting.) That means that I can fetch offset and row length directly from the owner's "row starts" array. With that, I can fetch the rows ("column indices") directly from the owner's storage. Underneath my (thin) API, we'll just be using MPI_Get() for these things. The overall semantics are collective, in that the owner needs to provide a send buffer and call MPI_Win_fence(), but it only provides one send buffer (no packing), and each process gets what it needs out of that buffer (by creating a suitable MPI_Datatype for MPI_Get()). The owner does not know how many procs accessed the data or what they accessed. > I was trying to avoid specifying an algorithm, but complete repartitioning > is inherently complex. I was thinking of a diffusive kind of thing with > nearest neighbors. > Sure, we can do that sort of thing. There are a variety of ways to "claim" an interface vertex, reducing with MAXLOC is one way.

1 0

Re: [petsc-dev] Generality of VecScatter
by Mark F. Adams 25 Nov '11

25 Nov '11

On Nov 25, 2011, at 10:51 AM, Jed Brown wrote: > On Fri, Nov 25, 2011 at 09:38, Mark F. Adams <mark.adams(a)columbia.edu> wrote: > My point is that at some point the rubber has to hit the road and you need two side information. You can then set up pointers for some fast one-sided stuff, and that is fine but you seem to agree that this is trivial. Its not clear to me that your model is not simply removing some (redundant) data from the user space, while it is probably in the MPI library. > > Whether this information is ever visible inside the MPI stack is dependent on the network. In general, it is never stored at once within the stack and the operations may be done without involving the CPU on the owning process (some network cards can do these operations directly on memory). If the MPI stack wanted to, it could log these operations, but I really meant it when I said it's not needed. > > Even with some funky implementation that built full two-sided information internally, I think the simpler user-level abstraction would be worthwhile. > > >> This information is communicated to the ghosters of that point (again without knowing how many there are), > > But you knew how many there were, and chose to forget it, and the MPI layer must know how many there are so it knows when its done communicating. I'm not sure how you insure correctness with your model. > > I never knew. We can ask: what is required to add another level of ghosting? OK, that sounds like a place to start. > With my model, the owner never needs to be informed that some procs are ghosting another level, nor be aware what those ghosted nodes are. It may place some data in an array for the "pointwise broadcast" of connectivity, but it doesn't need semantic knowledge that that information will used to increase the ghosting. Similarly, any process can stop ghosting a point without informing the owner in any way. I don't see how you can do that, you must have a different data model then me. We may need a white board for this but if I want to get an extra layer of ghosts I need to have the remote process tell me what they are. I have a distributed graph so I need to be told who my new ghosts are. > > Completion comes through MPI_Win_fence() which is collective, but not globally synchronizing (so neighbor communication must have completed, but it's not like an MPI_Barrier or MPI_Allreduce. >> I can assume that my fine grid was well partitioned, and to keep the process local assume the fine grid was partitioned with a space filling curve, or something, such that when I aggregate processors (eg, (0-3), (4-7) ... if I reduced the number of active processor by 4) I still have decent data locality. Eventually I would like a local algorithm to improve load balance at this point, but that is more than you might want to start with. >> >> If you have a local graph partitioner that can tell you what to move for better load balance, then I think I can move it efficiently and without redundant storage or specification. > > How do you load balance with a local partitioner? > > You said you wanted a local algorithm to improve load balance. I took that to mean that you would redistribute within a subcommunicator Ok, I see. > and, with simple communication, inform the rest of the processes of the new partition. I think this will work fine. > > What did you have in mind? I was trying to avoid specifying an algorithm, but complete repartitioning is inherently complex. I was thinking of a diffusive kind of thing with nearest neighbors.

1 0