[petsc-dev] PetscSection

[Barry Smith]

On Sep 11, 2013, at 3:47 AM, Matthew Knepley <knepley at gmail.com> wrote:

> On Tue, Sep 10, 2013 at 11:25 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
> On Sep 10, 2013, at 8:55 PM, Jed Brown <jedbrown at mcs.anl.gov> wrote:
> 
> > Matthew Knepley <knepley at gmail.com> writes:
> >>> PetscSF cannot currently scatter individual bytes (4 bytes minimum), and
> >>> even if it could, it's a horribly inefficient representation (4 or 8
> >>> bytes of metadata for each byte of payload).  The quick fix at that
> >>> moment was to send in units of size PetscInt (the struct was always
> >>> going to be divisible by that size).
> 
>    Since when are PETSc people interested in quick fixes, we are only concerned with doing things right; quick fixes lead to bad code and bad ideas.
> 
> >>>
> >>
> >> The way I understand it, he never had a single MPI_BYTE, but always a bunch.
> >> Shouldn't that packing handle that case?
> >
> > In DMPlexDistributeData, your fieldSF is blown up so that every unit has
> > its own metadata (rank, offset).  That's a terrible mechanism for moving
> > structs with hundreds of bytes.
> >
> > I would rather not make PetscSF deal with fully-heterogeneous data
> > (where each node can have a different size) because it makes indexing a
> > nightmare (you need a PetscSection or similar just to get in the door; I
> > don't want PetscSF to depend on PetscSection).
> 
>    Having PetscSF depend on PetscSection would be insane (after all PetscSection has information about constrained variables and other silly nonsense) but having PetscSF depend on XX makes perfect sense.
> 
> A review:
> 
> (current situation; with unimportant stuff removed)
> 
> +  sf - star forest
> .  nroots - number of root vertices on the current process (these are possible targets for other process to attach leaves)
> .  nleaves - number of leaf vertices on the current process, each of these references a root on any process
> .  ilocal - locations of leaves in leafdata buffers, pass NULL for contiguous storage
> .  iremote - remote locations of root vertices for each leaf on the current process
> 
> PetscSFSetGraph(PetscSF sf,PetscInt nroots,PetscInt nleaves,PetscInt *ilocal,,PetscSFNode *iremote,)
> 
> typedef struct {
>   PetscInt rank;                /* Rank of owner */  (why is this not PetscMPIInt?)
>   PetscInt index;               /* Index of node on rank */
> } PetscSFNode;
> 
> abstractly nleaves, ilocal  are the indices (into an array) you are copying too; iremote are the indices (into an array) you are copying from
> 
> Now if PETSc had an abstract concept of indexing (funny no one ever put that in PETSc decades ago) it would look like
> 
> PetscSFSetGraph(PetscSF sf,PetscInt nroots,  toIndices  , fromIndices)
> 
> but wait; PETSc does have have an abstraction for indexing (into regular arrays) called IS, come to think of it, PETSc has another abstraction for indexing (into arrays with different sized items) called XX.  Actually thinking a tiny bit more one realizes that there is a third: PetscSFNode is a way of indexing (into regular arrays) on a bunch of different processes. We come up with a couple more related, but inconsistent syntaxes for indexing, we'll have code has good as Trilinos.
> 
>    So let's fix this up! One abstraction for indexing that handles both regular arrays, arrays with different size items and both kinds of arrays on a bunch of different processes; with simple enough syntax so it can be used whenever indexing is needed including passing to PetscSF! (We can ignore the need for MPI datatypes in communication for now).
> 
>    Concrete examples to demonstrate this need not be so difficult. A completely general set of indices for this situation could be handled by
> 
> typedef struct {
>    PetscInt         nindices;
>    PetscMPIInt   *ranks;
>    PetscInt          *offsets;
>    PetscInt          *sizes;
> }  XXData;  Hardly more than PetscSFNode.  Or could be handled as an array of structs.
> 
> special cases would include:
>     all on one process so ranks is not needed
>     all items the same size so sizes not needed
>     contiquous memory (or strided) so offsets  not needed
> 
> If we decided not to go with an abstract object but instead just a concrete data structure that could handle all cases it could look something like
> 
> typedef struct {
>    PetscInt         nindices;
>    PetscMPIInt   *ranks;
>    PetscInt          *offsets;
>    PetscInt          *sizes;
>    PetscInt          chunksize,strideoffset;
> }  XXData;
> 
> For completeness I note with IS (for example in VecScatter creating with parallel vectors) we don't use rank,offset notation we use global offset == rstart[rank] + offset but that is something I'm sure we can deal with.
> 
> 
>     QED
> 
>     Comments?
> 
> This sounds like a good idea. So
> 
>   - we always store indices in blocks (get rid of regular IS)
>   - we can optimize for strided (get rid of strided IS)
>   - they can be local or nonlocal (get rid of PetscSFNode)
> 
> Right now, I have a GlobalSection, which uses XX but puts in different values. It has negative offsets
> and sizes for nonlocal indices. Now I could just allocate the ranks array and check for different rank
> rather than negative offset.
> 
> This would increase the size of GlobalSection, but then I could shared indices with PetscSF, so overall
> we could have a decrease. So it seems like PetscSection and PetscSF could share a substantial
> amount of code if we go to a common data structure.

   Sometimes it is more convenient to use "global indices" than rank/local offset. Thus we need a "scalable" utility to convert back and forth between them. Possibly
having interfaces to supply the information in either format.

For example in VecScatterCreate_PtoS()

 j      = 0;
  nsends = 0;
  for (i=0; i<nx; i++) {
    idx = bs*inidx[i];
    if (idx < owners[j]) j = 0;
    for (; j<size; j++) {
      if (idx < owners[j+1]) {
        if (!nprocs[j]++) nsends++;
        owner[i] = j;
        break;
      }
    }
    if (j == size) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ith %D block entry %D not owned by any process, upper bound %D",i,idx,owners[size]);
  }

Call the new unified beast PetscIS?

Barry


> 
> I believe all this is orthogonal to the packing issue, which I think can be done before we unify indices.
> In PetscSectionCreateFieldSF() I choose a granularity for communication. It would be easy to check
> for a common multiple of sizes here, and only check if MPI_Type is small. This is likely to be what we
> want to do in general, even creating our own packed types for communication. I will try and code it up.
> 
>    Matt
> 
> -- 
> 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