xref: /petsc/src/vec/is/sf/impls/basic/sfbasic.c (revision d083f849a86f1f43e18d534ee43954e2786cb29a)
1 
2 #include <petsc/private/sfimpl.h> /*I "petscsf.h" I*/
3 
4 typedef struct _n_PetscSFBasicPack *PetscSFBasicPack;
5 struct _n_PetscSFBasicPack {
6   void (*Pack)(PetscInt,PetscInt,const PetscInt*,const void*,void*);
7   void (*UnpackInsert)(PetscInt,PetscInt,const PetscInt*,void*,const void*);
8   void (*UnpackAdd)(PetscInt,PetscInt,const PetscInt*,void*,const void*);
9   void (*UnpackMin)(PetscInt,PetscInt,const PetscInt*,void*,const void*);
10   void (*UnpackMax)(PetscInt,PetscInt,const PetscInt*,void*,const void*);
11   void (*UnpackMinloc)(PetscInt,PetscInt,const PetscInt*,void*,const void*);
12   void (*UnpackMaxloc)(PetscInt,PetscInt,const PetscInt*,void*,const void*);
13   void (*UnpackMult)(PetscInt,PetscInt,const PetscInt*,void*,const void *);
14   void (*UnpackLAND)(PetscInt,PetscInt,const PetscInt*,void*,const void *);
15   void (*UnpackBAND)(PetscInt,PetscInt,const PetscInt*,void*,const void *);
16   void (*UnpackLOR)(PetscInt,PetscInt,const PetscInt*,void*,const void *);
17   void (*UnpackBOR)(PetscInt,PetscInt,const PetscInt*,void*,const void *);
18   void (*UnpackLXOR)(PetscInt,PetscInt,const PetscInt*,void*,const void *);
19   void (*UnpackBXOR)(PetscInt,PetscInt,const PetscInt*,void*,const void *);
20   void (*FetchAndInsert)(PetscInt,PetscInt,const PetscInt*,void*,void*);
21   void (*FetchAndAdd)(PetscInt,PetscInt,const PetscInt*,void*,void*);
22   void (*FetchAndMin)(PetscInt,PetscInt,const PetscInt*,void*,void*);
23   void (*FetchAndMax)(PetscInt,PetscInt,const PetscInt*,void*,void*);
24   void (*FetchAndMinloc)(PetscInt,PetscInt,const PetscInt*,void*,void*);
25   void (*FetchAndMaxloc)(PetscInt,PetscInt,const PetscInt*,void*,void*);
26   void (*FetchAndMult)(PetscInt,PetscInt,const PetscInt*,void*,void*);
27   void (*FetchAndLAND)(PetscInt,PetscInt,const PetscInt*,void*,void*);
28   void (*FetchAndBAND)(PetscInt,PetscInt,const PetscInt*,void*,void*);
29   void (*FetchAndLOR)(PetscInt,PetscInt,const PetscInt*,void*,void*);
30   void (*FetchAndBOR)(PetscInt,PetscInt,const PetscInt*,void*,void*);
31   void (*FetchAndLXOR)(PetscInt,PetscInt,const PetscInt*,void*,void*);
32   void (*FetchAndBXOR)(PetscInt,PetscInt,const PetscInt*,void*,void*);
33 
34   MPI_Datatype     unit;
35   PetscBool        isbuiltin;   /* Is unit an MPI builtin datatype? */
36   size_t           unitbytes;   /* Number of bytes in a unit */
37   PetscInt         bs;          /* Number of basic units in a unit */
38   const void       *key;        /* Array used as key for operation */
39   char             **root;      /* Packed root data, indexed by leaf rank */
40   char             **leaf;      /* Packed leaf data, indexed by root rank */
41   MPI_Request      *requests;   /* Array of root requests followed by leaf requests */
42   PetscSFBasicPack next;
43 };
44 
45 typedef struct {
46   PetscMPIInt      tag;
47   PetscMPIInt      niranks;     /* Number of incoming ranks (ranks accessing my roots) */
48   PetscMPIInt      ndiranks;    /* Number of incoming ranks (ranks accessing my roots) in distinguished set */
49   PetscMPIInt      *iranks;     /* Array of ranks that reference my roots */
50   PetscInt         itotal;      /* Total number of graph edges referencing my roots */
51   PetscInt         *ioffset;    /* Array of length niranks+1 holding offset in irootloc[] for each rank */
52   PetscInt         *irootloc;   /* Incoming roots referenced by ranks starting at ioffset[rank] */
53   PetscSFBasicPack avail;       /* One or more entries per MPI Datatype, lazily constructed */
54   PetscSFBasicPack inuse;       /* Buffers being used for transactions that have not yet completed */
55 } PetscSF_Basic;
56 
57 #if !defined(PETSC_HAVE_MPI_TYPE_DUP)
58 PETSC_STATIC_INLINE int MPI_Type_dup(MPI_Datatype datatype,MPI_Datatype *newtype)
59 {
60   int ierr;
61   ierr = MPI_Type_contiguous(1,datatype,newtype); if (ierr) return ierr;
62   ierr = MPI_Type_commit(newtype); if (ierr) return ierr;
63   return MPI_SUCCESS;
64 }
65 #endif
66 
67 /*
68  * MPI_Reduce_local is not really useful because it can't handle sparse data and it vectorizes "in the wrong direction",
69  * therefore we pack data types manually. This section defines packing routines for the standard data types.
70  */
71 
72 #define CPPJoin2_exp(a,b) a ## b
73 #define CPPJoin2(a,b) CPPJoin2_exp(a,b)
74 #define CPPJoin3_exp_(a,b,c) a ## b ## _ ## c
75 #define CPPJoin3_(a,b,c) CPPJoin3_exp_(a,b,c)
76 
77 /* Basic types without addition */
78 #define DEF_PackNoInit(type,BS)                                         \
79   static void CPPJoin3_(Pack_,type,BS)(PetscInt n,PetscInt bs,const PetscInt *idx,const void *unpacked,void *packed) { \
80     const type *u = (const type*)unpacked;                              \
81     type *p = (type*)packed;                                            \
82     PetscInt i,j,k;                                                     \
83     for (i=0; i<n; i++)                                                 \
84       for (j=0; j<bs; j+=BS)                                            \
85         for (k=j; k<j+BS; k++)                                          \
86           p[i*bs+k] = u[idx[i]*bs+k];                                   \
87   }                                                                     \
88   static void CPPJoin3_(UnpackInsert_,type,BS)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,const void *packed) { \
89     type *u = (type*)unpacked;                                          \
90     const type *p = (const type*)packed;                                \
91     PetscInt i,j,k;                                                     \
92     for (i=0; i<n; i++)                                                 \
93       for (j=0; j<bs; j+=BS)                                            \
94         for (k=j; k<j+BS; k++)                                          \
95           u[idx[i]*bs+k] = p[i*bs+k];                                   \
96   }                                                                     \
97   static void CPPJoin3_(FetchAndInsert_,type,BS)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,void *packed) { \
98     type *u = (type*)unpacked;                                          \
99     type *p = (type*)packed;                                            \
100     PetscInt i,j,k;                                                     \
101     for (i=0; i<n; i++) {                                               \
102       PetscInt ii = idx[i];                                             \
103       for (j=0; j<bs; j+=BS)                                            \
104         for (k=j; k<j+BS; k++) {                                        \
105           type t = u[ii*bs+k];                                          \
106           u[ii*bs+k] = p[i*bs+k];                                       \
107           p[i*bs+k] = t;                                                \
108         }                                                               \
109     }                                                                   \
110   }
111 
112 /* Basic types defining addition */
113 #define DEF_PackAddNoInit(type,BS)                                      \
114   DEF_PackNoInit(type,BS)                                               \
115   static void CPPJoin3_(UnpackAdd_,type,BS)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,const void *packed) { \
116     type *u = (type*)unpacked;                                          \
117     const type *p = (const type*)packed;                                \
118     PetscInt i,j,k;                                                     \
119     for (i=0; i<n; i++)                                                 \
120       for (j=0; j<bs; j+=BS)                                            \
121         for (k=j; k<j+BS; k++)                                          \
122           u[idx[i]*bs+k] += p[i*bs+k];                                  \
123   }                                                                     \
124   static void CPPJoin3_(FetchAndAdd_,type,BS)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,void *packed) { \
125     type *u = (type*)unpacked;                                          \
126     type *p = (type*)packed;                                            \
127     PetscInt i,j,k;                                                     \
128     for (i=0; i<n; i++) {                                               \
129       PetscInt ii = idx[i];                                             \
130       for (j=0; j<bs; j+=BS)                                            \
131         for (k=j; k<j+BS; k++) {                                        \
132           type t = u[ii*bs+k];                                          \
133           u[ii*bs+k] = t + p[i*bs+k];                                   \
134           p[i*bs+k] = t;                                                \
135         }                                                               \
136     }                                                                   \
137   }                                                                     \
138   static void CPPJoin3_(UnpackMult_,type,BS)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,const void *packed) { \
139     type *u = (type*)unpacked;                                          \
140     const type *p = (const type*)packed;                                \
141     PetscInt i,j,k;                                                     \
142     for (i=0; i<n; i++)                                                 \
143       for (j=0; j<bs; j+=BS)                                            \
144         for (k=j; k<j+BS; k++)                                          \
145           u[idx[i]*bs+k] *= p[i*bs+k];                                  \
146   }                                                                     \
147   static void CPPJoin3_(FetchAndMult_,type,BS)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,void *packed) { \
148     type *u = (type*)unpacked;                                          \
149     type *p = (type*)packed;                                            \
150     PetscInt i,j,k;                                                     \
151     for (i=0; i<n; i++) {                                               \
152       PetscInt ii = idx[i];                                             \
153       for (j=0; j<bs; j+=BS)                                            \
154         for (k=j; k<j+BS; k++) {                                        \
155           type t = u[ii*bs+k];                                          \
156           u[ii*bs+k] = t * p[i*bs+k];                                   \
157           p[i*bs+k] = t;                                                \
158         }                                                               \
159     }                                                                   \
160   }
161 #define DEF_Pack(type,BS)                                               \
162   DEF_PackAddNoInit(type,BS)                                            \
163   static void CPPJoin3_(PackInit_,type,BS)(PetscSFBasicPack link) {     \
164     link->Pack = CPPJoin3_(Pack_,type,BS);                              \
165     link->UnpackInsert = CPPJoin3_(UnpackInsert_,type,BS);              \
166     link->UnpackAdd = CPPJoin3_(UnpackAdd_,type,BS);                    \
167     link->UnpackMult = CPPJoin3_(UnpackMult_,type,BS);                  \
168     link->FetchAndInsert = CPPJoin3_(FetchAndInsert_,type,BS);          \
169     link->FetchAndAdd = CPPJoin3_(FetchAndAdd_,type,BS);                \
170     link->FetchAndMult = CPPJoin3_(FetchAndMult_,type,BS);              \
171     link->unitbytes = sizeof(type);                                     \
172   }
173 /* Comparable types */
174 #define DEF_PackCmp(type)                                               \
175   DEF_PackAddNoInit(type,1)                                             \
176   static void CPPJoin2(UnpackMax_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,const void *packed) { \
177     type *u = (type*)unpacked;                                          \
178     const type *p = (const type*)packed;                                \
179     PetscInt i;                                                         \
180     for (i=0; i<n; i++) {                                               \
181       type v = u[idx[i]];                                               \
182       u[idx[i]] = PetscMax(v,p[i]);                                     \
183     }                                                                   \
184   }                                                                     \
185   static void CPPJoin2(UnpackMin_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,const void *packed) { \
186     type *u = (type*)unpacked;                                          \
187     const type *p = (const type*)packed;                                \
188     PetscInt i;                                                         \
189     for (i=0; i<n; i++) {                                               \
190       type v = u[idx[i]];                                               \
191       u[idx[i]] = PetscMin(v,p[i]);                                     \
192     }                                                                   \
193   }                                                                     \
194   static void CPPJoin2(FetchAndMax_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,void *packed) { \
195     type *u = (type*)unpacked;                                          \
196     type *p = (type*)packed;                                            \
197     PetscInt i;                                                         \
198     for (i=0; i<n; i++) {                                               \
199       PetscInt j = idx[i];                                              \
200       type v = u[j];                                                    \
201       u[j] = PetscMax(v,p[i]);                                          \
202       p[i] = v;                                                         \
203     }                                                                   \
204   }                                                                     \
205   static void CPPJoin2(FetchAndMin_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,void *packed) { \
206     type *u = (type*)unpacked;                                          \
207     type *p = (type*)packed;                                            \
208     PetscInt i;                                                         \
209     for (i=0; i<n; i++) {                                               \
210       PetscInt j = idx[i];                                              \
211       type v = u[j];                                                    \
212       u[j] = PetscMin(v,p[i]);                                          \
213       p[i] = v;                                                         \
214     }                                                                   \
215   }                                                                     \
216   static void CPPJoin2(PackInit_,type)(PetscSFBasicPack link) {         \
217     link->Pack = CPPJoin3_(Pack_,type,1);                               \
218     link->UnpackInsert = CPPJoin3_(UnpackInsert_,type,1);               \
219     link->UnpackAdd  = CPPJoin3_(UnpackAdd_,type,1);                    \
220     link->UnpackMax  = CPPJoin2(UnpackMax_,type);                       \
221     link->UnpackMin  = CPPJoin2(UnpackMin_,type);                       \
222     link->UnpackMult = CPPJoin3_(UnpackMult_,type,1);                   \
223     link->FetchAndInsert = CPPJoin3_(FetchAndInsert_,type,1);           \
224     link->FetchAndAdd = CPPJoin3_(FetchAndAdd_ ,type,1);                \
225     link->FetchAndMax = CPPJoin2(FetchAndMax_ ,type);                   \
226     link->FetchAndMin = CPPJoin2(FetchAndMin_ ,type);                   \
227     link->FetchAndMult = CPPJoin3_(FetchAndMult_,type,1);               \
228     link->unitbytes = sizeof(type);                                     \
229   }
230 
231 /* Logical Types */
232 #define DEF_PackLog(type)                                               \
233   static void CPPJoin2(UnpackLAND_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,const void *packed) { \
234     type *u = (type*)unpacked;                                          \
235     const type *p = (const type*)packed;                                \
236     PetscInt i;                                                         \
237     for (i=0; i<n; i++) {                                               \
238       type v = u[idx[i]];                                               \
239       u[idx[i]] = v && p[i];                                            \
240     }                                                                   \
241   }                                                                     \
242   static void CPPJoin2(UnpackLOR_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,const void *packed) { \
243     type *u = (type*)unpacked;                                          \
244     const type *p = (const type*)packed;                                \
245     PetscInt i;                                                         \
246     for (i=0; i<n; i++) {                                               \
247       type v = u[idx[i]];                                               \
248       u[idx[i]] = v || p[i];                                            \
249     }                                                                   \
250   }                                                                     \
251   static void CPPJoin2(UnpackLXOR_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,const void *packed) { \
252     type *u = (type*)unpacked;                                          \
253     const type *p = (const type*)packed;                                \
254     PetscInt i;                                                         \
255     for (i=0; i<n; i++) {                                               \
256       type v = u[idx[i]];                                               \
257       u[idx[i]] = (!v)!=(!p[i]);                                        \
258     }                                                                   \
259   }                                                                     \
260   static void CPPJoin2(FetchAndLAND_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,void *packed) { \
261     type *u = (type*)unpacked;                                          \
262     type *p = (type*)packed;                                            \
263     PetscInt i;                                                         \
264     for (i=0; i<n; i++) {                                               \
265       PetscInt j = idx[i];                                              \
266       type v = u[j];                                                    \
267       u[j] = v && p[i];                                                 \
268       p[i] = v;                                                         \
269     }                                                                   \
270   }                                                                     \
271   static void CPPJoin2(FetchAndLOR_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,void *packed) { \
272     type *u = (type*)unpacked;                                          \
273     type *p = (type*)packed;                                            \
274     PetscInt i;                                                         \
275     for (i=0; i<n; i++) {                                               \
276       PetscInt j = idx[i];                                              \
277       type v = u[j];                                                    \
278       u[j] = v || p[i];                                                 \
279       p[i] = v;                                                         \
280     }                                                                   \
281   }                                                                     \
282   static void CPPJoin2(FetchAndLXOR_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,void *packed) { \
283     type *u = (type*)unpacked;                                          \
284     type *p = (type*)packed;                                            \
285     PetscInt i;                                                         \
286     for (i=0; i<n; i++) {                                               \
287       PetscInt j = idx[i];                                              \
288       type v = u[j];                                                    \
289       u[j] = (!v)!=(!p[i]);                                             \
290       p[i] = v;                                                         \
291     }                                                                   \
292   }                                                                     \
293   static void CPPJoin2(PackInit_Logical_,type)(PetscSFBasicPack link) { \
294     link->UnpackLAND = CPPJoin2(UnpackLAND_,type);                      \
295     link->UnpackLOR  = CPPJoin2(UnpackLOR_,type);                       \
296     link->UnpackLXOR = CPPJoin2(UnpackLXOR_,type);                      \
297     link->FetchAndLAND = CPPJoin2(FetchAndLAND_,type);                  \
298     link->FetchAndLOR  = CPPJoin2(FetchAndLOR_,type);                   \
299     link->FetchAndLXOR = CPPJoin2(FetchAndLXOR_,type);                  \
300   }
301 
302 
303 /* Bitwise Types */
304 #define DEF_PackBit(type)                                               \
305   static void CPPJoin2(UnpackBAND_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,const void *packed) { \
306     type *u = (type*)unpacked;                                          \
307     const type *p = (const type*)packed;                                \
308     PetscInt i;                                                         \
309     for (i=0; i<n; i++) {                                               \
310       type v = u[idx[i]];                                               \
311       u[idx[i]] = v & p[i];                                             \
312     }                                                                   \
313   }                                                                     \
314   static void CPPJoin2(UnpackBOR_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,const void *packed) { \
315     type *u = (type*)unpacked;                                          \
316     const type *p = (const type*)packed;                                \
317     PetscInt i;                                                         \
318     for (i=0; i<n; i++) {                                               \
319       type v = u[idx[i]];                                               \
320       u[idx[i]] = v | p[i];                                             \
321     }                                                                   \
322   }                                                                     \
323   static void CPPJoin2(UnpackBXOR_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,const void *packed) { \
324     type *u = (type*)unpacked;                                          \
325     const type *p = (const type*)packed;                                \
326     PetscInt i;                                                         \
327     for (i=0; i<n; i++) {                                               \
328       type v = u[idx[i]];                                               \
329       u[idx[i]] = v^p[i];                                               \
330     }                                                                   \
331   }                                                                     \
332   static void CPPJoin2(FetchAndBAND_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,void *packed) { \
333     type *u = (type*)unpacked;                                          \
334     type *p = (type*)packed;                                            \
335     PetscInt i;                                                         \
336     for (i=0; i<n; i++) {                                               \
337       PetscInt j = idx[i];                                              \
338       type v = u[j];                                                    \
339       u[j] = v & p[i];                                                  \
340       p[i] = v;                                                         \
341     }                                                                   \
342   }                                                                     \
343   static void CPPJoin2(FetchAndBOR_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,void *packed) { \
344     type *u = (type*)unpacked;                                          \
345     type *p = (type*)packed;                                            \
346     PetscInt i;                                                         \
347     for (i=0; i<n; i++) {                                               \
348       PetscInt j = idx[i];                                              \
349       type v = u[j];                                                    \
350       u[j] = v | p[i];                                                  \
351       p[i] = v;                                                         \
352     }                                                                   \
353   }                                                                     \
354   static void CPPJoin2(FetchAndBXOR_,type)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,void *packed) { \
355     type *u = (type*)unpacked;                                          \
356     type *p = (type*)packed;                                            \
357     PetscInt i;                                                         \
358     for (i=0; i<n; i++) {                                               \
359       PetscInt j = idx[i];                                              \
360       type v = u[j];                                                    \
361       u[j] = v^p[i];                                                    \
362       p[i] = v;                                                         \
363     }                                                                   \
364   }                                                                     \
365   static void CPPJoin2(PackInit_Bitwise_,type)(PetscSFBasicPack link) { \
366     link->UnpackBAND = CPPJoin2(UnpackBAND_,type);                      \
367     link->UnpackBOR  = CPPJoin2(UnpackBOR_,type);                       \
368     link->UnpackBXOR = CPPJoin2(UnpackBXOR_,type);                      \
369     link->FetchAndBAND = CPPJoin2(FetchAndBAND_,type);                  \
370     link->FetchAndBOR  = CPPJoin2(FetchAndBOR_,type);                   \
371     link->FetchAndBXOR = CPPJoin2(FetchAndBXOR_,type);                  \
372   }
373 
374 /* Pair types */
375 #define CPPJoinloc_exp(base,op,t1,t2) base ## op ## loc_ ## t1 ## _ ## t2
376 #define CPPJoinloc(base,op,t1,t2) CPPJoinloc_exp(base,op,t1,t2)
377 #define PairType(type1,type2) CPPJoin3_(_pairtype_,type1,type2)
378 #define DEF_UnpackXloc(type1,type2,locname,op)                              \
379   static void CPPJoinloc(Unpack,locname,type1,type2)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,const void *packed) { \
380     PairType(type1,type2) *u = (PairType(type1,type2)*)unpacked;        \
381     const PairType(type1,type2) *p = (const PairType(type1,type2)*)packed; \
382     PetscInt i;                                                         \
383     for (i=0; i<n; i++) {                                               \
384       PetscInt j = idx[i];                                              \
385       if (p[i].a op u[j].a) {                                           \
386         u[j].a = p[i].a;                                                \
387         u[j].b = p[i].b;                                                \
388       } else if (u[j].a == p[i].a) {                                    \
389         u[j].b = PetscMin(u[j].b,p[i].b);                               \
390       }                                                                 \
391     }                                                                   \
392   }                                                                     \
393   static void CPPJoinloc(FetchAnd,locname,type1,type2)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,void *packed) { \
394     PairType(type1,type2) *u = (PairType(type1,type2)*)unpacked;        \
395     PairType(type1,type2) *p = (PairType(type1,type2)*)packed;          \
396     PetscInt i;                                                         \
397     for (i=0; i<n; i++) {                                               \
398       PetscInt j = idx[i];                                              \
399       PairType(type1,type2) v;                                          \
400       v.a = u[j].a;                                                     \
401       v.b = u[j].b;                                                     \
402       if (p[i].a op u[j].a) {                                           \
403         u[j].a = p[i].a;                                                \
404         u[j].b = p[i].b;                                                \
405       } else if (u[j].a == p[i].a) {                                    \
406         u[j].b = PetscMin(u[j].b,p[i].b);                               \
407       }                                                                 \
408       p[i].a = v.a;                                                     \
409       p[i].b = v.b;                                                     \
410     }                                                                   \
411   }
412 #define DEF_PackPair(type1,type2)                                       \
413   typedef struct {type1 a; type2 b;} PairType(type1,type2);             \
414   static void CPPJoin3_(Pack_,type1,type2)(PetscInt n,PetscInt bs,const PetscInt *idx,const void *unpacked,void *packed) { \
415     const PairType(type1,type2) *u = (const PairType(type1,type2)*)unpacked; \
416     PairType(type1,type2) *p = (PairType(type1,type2)*)packed;          \
417     PetscInt i;                                                         \
418     for (i=0; i<n; i++) {                                               \
419       p[i].a = u[idx[i]].a;                                             \
420       p[i].b = u[idx[i]].b;                                             \
421     }                                                                   \
422   }                                                                     \
423   static void CPPJoin3_(UnpackInsert_,type1,type2)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,const void *packed) { \
424     PairType(type1,type2) *u = (PairType(type1,type2)*)unpacked;       \
425     const PairType(type1,type2) *p = (const PairType(type1,type2)*)packed; \
426     PetscInt i;                                                         \
427     for (i=0; i<n; i++) {                                               \
428       u[idx[i]].a = p[i].a;                                             \
429       u[idx[i]].b = p[i].b;                                             \
430     }                                                                   \
431   }                                                                     \
432   static void CPPJoin3_(UnpackAdd_,type1,type2)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,const void *packed) { \
433     PairType(type1,type2) *u = (PairType(type1,type2)*)unpacked;       \
434     const PairType(type1,type2) *p = (const PairType(type1,type2)*)packed; \
435     PetscInt i;                                                         \
436     for (i=0; i<n; i++) {                                               \
437       u[idx[i]].a += p[i].a;                                            \
438       u[idx[i]].b += p[i].b;                                            \
439     }                                                                   \
440   }                                                                     \
441   static void CPPJoin3_(FetchAndInsert_,type1,type2)(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,void *packed) { \
442     PairType(type1,type2) *u = (PairType(type1,type2)*)unpacked;        \
443     PairType(type1,type2) *p = (PairType(type1,type2)*)packed;          \
444     PetscInt i;                                                         \
445     for (i=0; i<n; i++) {                                               \
446       PetscInt j = idx[i];                                              \
447       PairType(type1,type2) v;                                          \
448       v.a = u[j].a;                                                     \
449       v.b = u[j].b;                                                     \
450       u[j].a = p[i].a;                                                  \
451       u[j].b = p[i].b;                                                  \
452       p[i].a = v.a;                                                     \
453       p[i].b = v.b;                                                     \
454     }                                                                   \
455   }                                                                     \
456   static void FetchAndAdd_ ## type1 ## _ ## type2(PetscInt n,PetscInt bs,const PetscInt *idx,void *unpacked,void *packed) { \
457     PairType(type1,type2) *u = (PairType(type1,type2)*)unpacked;       \
458     PairType(type1,type2) *p = (PairType(type1,type2)*)packed;         \
459     PetscInt i;                                                         \
460     for (i=0; i<n; i++) {                                               \
461       PetscInt j = idx[i];                                              \
462       PairType(type1,type2) v;                                          \
463       v.a = u[j].a;                                                     \
464       v.b = u[j].b;                                                     \
465       u[j].a = v.a + p[i].a;                                            \
466       u[j].b = v.b + p[i].b;                                            \
467       p[i].a = v.a;                                                     \
468       p[i].b = v.b;                                                     \
469     }                                                                   \
470   }                                                                     \
471   DEF_UnpackXloc(type1,type2,Max,>)                                     \
472   DEF_UnpackXloc(type1,type2,Min,<)                                     \
473   static void CPPJoin3_(PackInit_,type1,type2)(PetscSFBasicPack link) { \
474     link->Pack = CPPJoin3_(Pack_,type1,type2);                          \
475     link->UnpackInsert = CPPJoin3_(UnpackInsert_,type1,type2);          \
476     link->UnpackAdd = CPPJoin3_(UnpackAdd_,type1,type2);                \
477     link->UnpackMaxloc = CPPJoin3_(UnpackMaxloc_,type1,type2);          \
478     link->UnpackMinloc = CPPJoin3_(UnpackMinloc_,type1,type2);          \
479     link->FetchAndInsert = CPPJoin3_(FetchAndInsert_,type1,type2);      \
480     link->FetchAndAdd = CPPJoin3_(FetchAndAdd_,type1,type2);            \
481     link->FetchAndMaxloc = CPPJoin3_(FetchAndMaxloc_,type1,type2);      \
482     link->FetchAndMinloc = CPPJoin3_(FetchAndMinloc_,type1,type2);      \
483     link->unitbytes = sizeof(PairType(type1,type2));                    \
484   }
485 
486 /* Currently only dumb blocks of data */
487 #define BlockType(unit,count) CPPJoin3_(_blocktype_,unit,count)
488 #define DEF_Block(unit,count)                                           \
489   typedef struct {unit v[count];} BlockType(unit,count);                \
490   DEF_PackNoInit(BlockType(unit,count),1)                               \
491   static void CPPJoin3_(PackInit_block_,unit,count)(PetscSFBasicPack link) { \
492     link->Pack = CPPJoin3_(Pack_,BlockType(unit,count),1);               \
493     link->UnpackInsert = CPPJoin3_(UnpackInsert_,BlockType(unit,count),1); \
494     link->FetchAndInsert = CPPJoin3_(FetchAndInsert_,BlockType(unit,count),1); \
495     link->unitbytes = sizeof(BlockType(unit,count));                    \
496   }
497 
498 /* The typedef is used to get a typename without space that CPPJoin can handle */
499 typedef signed char SignedChar;
500 typedef unsigned char UnsignedChar;
501 
502 DEF_PackCmp(SignedChar)
503 DEF_PackBit(SignedChar)
504 DEF_PackLog(SignedChar)
505 DEF_PackCmp(UnsignedChar)
506 DEF_PackBit(UnsignedChar)
507 DEF_PackLog(UnsignedChar)
508 DEF_PackCmp(int)
509 DEF_PackBit(int)
510 DEF_PackLog(int)
511 DEF_PackCmp(PetscInt)
512 DEF_PackBit(PetscInt)
513 DEF_PackLog(PetscInt)
514 DEF_Pack(PetscInt,2)
515 DEF_Pack(PetscInt,3)
516 DEF_Pack(PetscInt,4)
517 DEF_Pack(PetscInt,5)
518 DEF_Pack(PetscInt,7)
519 DEF_PackCmp(PetscReal)
520 DEF_PackLog(PetscReal)
521 DEF_Pack(PetscReal,2)
522 DEF_Pack(PetscReal,3)
523 DEF_Pack(PetscReal,4)
524 DEF_Pack(PetscReal,5)
525 DEF_Pack(PetscReal,7)
526 #if defined(PETSC_HAVE_COMPLEX)
527 DEF_Pack(PetscComplex,1)
528 DEF_Pack(PetscComplex,2)
529 DEF_Pack(PetscComplex,3)
530 DEF_Pack(PetscComplex,4)
531 DEF_Pack(PetscComplex,5)
532 DEF_Pack(PetscComplex,7)
533 #endif
534 DEF_PackPair(int,int)
535 DEF_PackPair(PetscInt,PetscInt)
536 DEF_Block(int,1)
537 DEF_Block(int,2)
538 DEF_Block(int,3)
539 DEF_Block(int,4)
540 DEF_Block(int,5)
541 DEF_Block(int,6)
542 DEF_Block(int,7)
543 DEF_Block(int,8)
544 DEF_Block(char,1)
545 DEF_Block(char,2)
546 DEF_Block(char,3)
547 #if PETSC_SIZEOF_INT == 8
548 DEF_Block(char,4)
549 DEF_Block(char,5)
550 DEF_Block(char,6)
551 DEF_Block(char,7)
552 #endif
553 
554 static PetscErrorCode PetscSFSetUp_Basic(PetscSF sf)
555 {
556   PetscSF_Basic  *bas = (PetscSF_Basic*)sf->data;
557   PetscErrorCode ierr;
558   PetscInt       *rlengths,*ilengths,i;
559   PetscMPIInt    rank,niranks,*iranks;
560   MPI_Comm       comm;
561   MPI_Group      group;
562   PetscMPIInt    nreqs = 0;
563   MPI_Request    *reqs;
564 
565   PetscFunctionBegin;
566   ierr = MPI_Comm_group(PETSC_COMM_SELF,&group);CHKERRQ(ierr);
567   ierr = PetscSFSetUpRanks(sf,group);CHKERRQ(ierr);
568   ierr = MPI_Group_free(&group);CHKERRQ(ierr);
569   ierr = PetscObjectGetComm((PetscObject)sf,&comm);CHKERRQ(ierr);
570   ierr = PetscObjectGetNewTag((PetscObject)sf,&bas->tag);CHKERRQ(ierr);
571   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
572   /*
573    * Inform roots about how many leaves and from which ranks
574    */
575   ierr = PetscMalloc1(sf->nranks,&rlengths);CHKERRQ(ierr);
576   /* Determine number, sending ranks, and length of incoming */
577   for (i=0; i<sf->nranks; i++) {
578     rlengths[i] = sf->roffset[i+1] - sf->roffset[i]; /* Number of roots referenced by my leaves; for rank sf->ranks[i] */
579   }
580   ierr = PetscCommBuildTwoSided(comm,1,MPIU_INT,sf->nranks-sf->ndranks,sf->ranks+sf->ndranks,rlengths+sf->ndranks,&niranks,&iranks,&ilengths);CHKERRQ(ierr);
581 
582   /* Sort iranks. See use of VecScatterGetRemoteOrdered_Private() in MatGetBrowsOfAoCols_MPIAIJ() on why.
583      We could sort ranks there at the price of allocating extra working arrays. Presumably, niranks is
584      small and the sorting is cheap.
585    */
586   ierr = PetscSortMPIIntWithIntArray(niranks,iranks,ilengths);CHKERRQ(ierr);
587 
588   /* Partition into distinguished and non-distinguished incoming ranks */
589   bas->ndiranks = sf->ndranks;
590   bas->niranks = bas->ndiranks + niranks;
591   ierr = PetscMalloc2(bas->niranks,&bas->iranks,bas->niranks+1,&bas->ioffset);CHKERRQ(ierr);
592   bas->ioffset[0] = 0;
593   for (i=0; i<bas->ndiranks; i++) {
594     bas->iranks[i] = sf->ranks[i];
595     bas->ioffset[i+1] = bas->ioffset[i] + rlengths[i];
596   }
597   for (i=bas->ndiranks; i<bas->niranks; i++) {
598     bas->iranks[i] = iranks[i-bas->ndiranks];
599     bas->ioffset[i+1] = bas->ioffset[i] + ilengths[i-bas->ndiranks];
600   }
601   bas->itotal = bas->ioffset[i];
602   ierr = PetscFree(rlengths);CHKERRQ(ierr);
603   ierr = PetscFree(iranks);CHKERRQ(ierr);
604   ierr = PetscFree(ilengths);CHKERRQ(ierr);
605 
606   /* Sanity checks for distinguished ranks */
607   if (sf->ndranks != bas->ndiranks) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Broken setup for shared ranks");
608   if (sf->ndranks > 1 || (sf->ndranks == 1 && sf->ranks[0] != rank)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Broken setup for shared ranks");
609   if (bas->ndiranks > 1 || (bas->ndiranks == 1 && bas->iranks[0] != rank)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Broken setup for shared ranks");
610 
611   /* Send leaf identities to roots */
612   ierr = PetscMalloc1(bas->itotal,&bas->irootloc);CHKERRQ(ierr);
613   ierr = PetscMalloc1(bas->niranks-bas->ndiranks+sf->nranks-sf->ndranks,&reqs);CHKERRQ(ierr);
614   for (i=sf->ndranks; i<sf->nranks; i++, nreqs++) {
615     PetscMPIInt npoints;
616     ierr = PetscMPIIntCast(sf->roffset[i+1]-sf->roffset[i],&npoints);CHKERRQ(ierr);
617     ierr = MPI_Isend(sf->rremote+sf->roffset[i],npoints,MPIU_INT,sf->ranks[i],bas->tag,comm,&reqs[nreqs]);CHKERRQ(ierr);
618   }
619   for (i=bas->ndiranks; i<bas->niranks; i++, nreqs++) {
620     PetscMPIInt npoints;
621     ierr = PetscMPIIntCast(bas->ioffset[i+1]-bas->ioffset[i],&npoints);CHKERRQ(ierr);
622     ierr = MPI_Irecv(bas->irootloc+bas->ioffset[i],npoints,MPIU_INT,bas->iranks[i],bas->tag,comm,&reqs[nreqs]);CHKERRQ(ierr);
623   }
624   for (i=0; i<sf->ndranks; i++) {
625     PetscInt npoints = sf->roffset[i+1]-sf->roffset[i];
626     if (npoints != bas->ioffset[i+1]-bas->ioffset[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Distinguished rank exchange has mismatched lengths");
627     ierr = PetscMemcpy(bas->irootloc+bas->ioffset[i],sf->rremote+sf->roffset[i],npoints*sizeof(bas->irootloc[0]));CHKERRQ(ierr);
628   }
629   ierr = MPI_Waitall(nreqs,reqs,MPI_STATUSES_IGNORE);CHKERRQ(ierr);
630   ierr = PetscFree(reqs);CHKERRQ(ierr);
631   PetscFunctionReturn(0);
632 }
633 
634 static PetscErrorCode PetscSFBasicPackTypeSetup(PetscSFBasicPack link,MPI_Datatype unit)
635 {
636   PetscErrorCode ierr;
637   PetscBool      isInt,isPetscInt,isPetscReal,is2Int,is2PetscInt,isSignedChar,isUnsignedChar;
638   PetscInt       nPetscIntContig,nPetscRealContig;
639   PetscMPIInt    ni,na,nd,combiner;
640 #if defined(PETSC_HAVE_COMPLEX)
641   PetscBool isPetscComplex;
642   PetscInt nPetscComplexContig;
643 #endif
644 
645   PetscFunctionBegin;
646   ierr = MPIPetsc_Type_compare(unit,MPI_SIGNED_CHAR,&isSignedChar);CHKERRQ(ierr);
647   ierr = MPIPetsc_Type_compare(unit,MPI_UNSIGNED_CHAR,&isUnsignedChar);CHKERRQ(ierr);
648   /* MPI_CHAR is treated below as a dumb block type that does not support reduction according to MPI standard */
649   ierr = MPIPetsc_Type_compare(unit,MPI_INT,&isInt);CHKERRQ(ierr);
650   ierr = MPIPetsc_Type_compare(unit,MPIU_INT,&isPetscInt);CHKERRQ(ierr);
651   ierr = MPIPetsc_Type_compare_contig(unit,MPIU_INT,&nPetscIntContig);CHKERRQ(ierr);
652   ierr = MPIPetsc_Type_compare(unit,MPIU_REAL,&isPetscReal);CHKERRQ(ierr);
653   ierr = MPIPetsc_Type_compare_contig(unit,MPIU_REAL,&nPetscRealContig);CHKERRQ(ierr);
654 #if defined(PETSC_HAVE_COMPLEX)
655   ierr = MPIPetsc_Type_compare(unit,MPIU_COMPLEX,&isPetscComplex);CHKERRQ(ierr);
656   ierr = MPIPetsc_Type_compare_contig(unit,MPIU_COMPLEX,&nPetscComplexContig);CHKERRQ(ierr);
657 #endif
658   ierr = MPIPetsc_Type_compare(unit,MPI_2INT,&is2Int);CHKERRQ(ierr);
659   ierr = MPIPetsc_Type_compare(unit,MPIU_2INT,&is2PetscInt);CHKERRQ(ierr);
660   ierr = MPI_Type_get_envelope(unit,&ni,&na,&nd,&combiner);CHKERRQ(ierr);
661   link->isbuiltin = (combiner == MPI_COMBINER_NAMED) ? PETSC_TRUE : PETSC_FALSE;
662   link->bs = 1;
663 
664   if (isSignedChar) {PackInit_SignedChar(link); PackInit_Logical_SignedChar(link); PackInit_Bitwise_SignedChar(link);}
665   else if (isUnsignedChar) {PackInit_UnsignedChar(link); PackInit_Logical_UnsignedChar(link); PackInit_Bitwise_UnsignedChar(link);}
666   else if (isInt) {PackInit_int(link); PackInit_Logical_int(link); PackInit_Bitwise_int(link);}
667   else if (isPetscInt) {PackInit_PetscInt(link); PackInit_Logical_PetscInt(link); PackInit_Bitwise_PetscInt(link);}
668   else if (isPetscReal) {PackInit_PetscReal(link); PackInit_Logical_PetscReal(link);}
669 #if defined(PETSC_HAVE_COMPLEX)
670   else if (isPetscComplex) PackInit_PetscComplex_1(link);
671 #endif
672   else if (is2Int) PackInit_int_int(link);
673   else if (is2PetscInt) PackInit_PetscInt_PetscInt(link);
674   else if (nPetscIntContig) {
675     if (nPetscIntContig%7 == 0) PackInit_PetscInt_7(link);
676     else if (nPetscIntContig%5 == 0) PackInit_PetscInt_5(link);
677     else if (nPetscIntContig%4 == 0) PackInit_PetscInt_4(link);
678     else if (nPetscIntContig%3 == 0) PackInit_PetscInt_3(link);
679     else if (nPetscIntContig%2 == 0) PackInit_PetscInt_2(link);
680     else PackInit_PetscInt(link);
681     link->bs = nPetscIntContig;
682     link->unitbytes *= nPetscIntContig;
683   } else if (nPetscRealContig) {
684     if (nPetscRealContig%7 == 0) PackInit_PetscReal_7(link);
685     else if (nPetscRealContig%5 == 0) PackInit_PetscReal_5(link);
686     else if (nPetscRealContig%4 == 0) PackInit_PetscReal_4(link);
687     else if (nPetscRealContig%3 == 0) PackInit_PetscReal_3(link);
688     else if (nPetscRealContig%2 == 0) PackInit_PetscReal_2(link);
689     else PackInit_PetscReal(link);
690     link->bs = nPetscRealContig;
691     link->unitbytes *= nPetscRealContig;
692 #if defined(PETSC_HAVE_COMPLEX)
693   } else if (nPetscComplexContig) {
694     if (nPetscComplexContig%7 == 0) PackInit_PetscComplex_7(link);
695     else if (nPetscComplexContig%5 == 0) PackInit_PetscComplex_5(link);
696     else if (nPetscComplexContig%4 == 0) PackInit_PetscComplex_4(link);
697     else if (nPetscComplexContig%3 == 0) PackInit_PetscComplex_3(link);
698     else if (nPetscComplexContig%2 == 0) PackInit_PetscComplex_2(link);
699     else PackInit_PetscComplex_1(link);
700     link->bs = nPetscComplexContig;
701     link->unitbytes *= nPetscComplexContig;
702 #endif
703   } else {
704     MPI_Aint lb,bytes;
705     ierr = MPI_Type_get_extent(unit,&lb,&bytes);CHKERRQ(ierr);
706     if (lb != 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Datatype with nonzero lower bound %ld\n",(long)lb);
707     if (bytes % sizeof(int)) { /* If the type size is not multiple of int */
708 #if PETSC_SIZEOF_INT == 8
709       if      (bytes%7 == 0) {PackInit_block_char_7(link); link->bs = bytes/7;} /* Note the basic type is char[7] */
710       else if (bytes%6 == 0) {PackInit_block_char_6(link); link->bs = bytes/6;}
711       else if (bytes%5 == 0) {PackInit_block_char_5(link); link->bs = bytes/5;}
712       else if (bytes%4 == 0) {PackInit_block_char_4(link); link->bs = bytes/4;}
713       else
714 #endif
715       if      (bytes%3 == 0) {PackInit_block_char_3(link); link->bs = bytes/3;}
716       else if (bytes%2 == 0) {PackInit_block_char_2(link); link->bs = bytes/2;}
717       else                   {PackInit_block_char_1(link); link->bs = bytes/1;}
718       link->unitbytes = bytes;
719     } else {
720       PetscInt nInt = bytes / sizeof(int);
721       if      (nInt%8 == 0)  {PackInit_block_int_8(link);  link->bs = nInt/8;} /* Note the basic type is int[8] */
722       else if (nInt%7 == 0)  {PackInit_block_int_7(link);  link->bs = nInt/7;}
723       else if (nInt%6 == 0)  {PackInit_block_int_6(link);  link->bs = nInt/6;}
724       else if (nInt%5 == 0)  {PackInit_block_int_5(link);  link->bs = nInt/5;}
725       else if (nInt%4 == 0)  {PackInit_block_int_4(link);  link->bs = nInt/4;}
726       else if (nInt%3 == 0)  {PackInit_block_int_3(link);  link->bs = nInt/3;}
727       else if (nInt%2 == 0)  {PackInit_block_int_2(link);  link->bs = nInt/2;}
728       else                   {PackInit_block_int_1(link);  link->bs = nInt/1;}
729       link->unitbytes = bytes;
730     }
731   }
732   if (link->isbuiltin) link->unit = unit; /* builtin datatypes are common. Make it fast */
733   else {ierr = MPI_Type_dup(unit,&link->unit);CHKERRQ(ierr);}
734   PetscFunctionReturn(0);
735 }
736 
737 static PetscErrorCode PetscSFBasicPackGetUnpackOp(PetscSF sf,PetscSFBasicPack link,MPI_Op op,void (**UnpackOp)(PetscInt,PetscInt,const PetscInt*,void*,const void*))
738 {
739   PetscFunctionBegin;
740   *UnpackOp = NULL;
741   if (op == MPIU_REPLACE) *UnpackOp = link->UnpackInsert;
742   else if (op == MPI_SUM || op == MPIU_SUM) *UnpackOp = link->UnpackAdd;
743   else if (op == MPI_PROD) *UnpackOp = link->UnpackMult;
744   else if (op == MPI_MAX || op == MPIU_MAX) *UnpackOp = link->UnpackMax;
745   else if (op == MPI_MIN || op == MPIU_MIN) *UnpackOp = link->UnpackMin;
746   else if (op == MPI_LAND) *UnpackOp = link->UnpackLAND;
747   else if (op == MPI_BAND) *UnpackOp = link->UnpackBAND;
748   else if (op == MPI_LOR) *UnpackOp = link->UnpackLOR;
749   else if (op == MPI_BOR) *UnpackOp = link->UnpackBOR;
750   else if (op == MPI_LXOR) *UnpackOp = link->UnpackLXOR;
751   else if (op == MPI_BXOR) *UnpackOp = link->UnpackBXOR;
752   else if (op == MPI_MAXLOC) *UnpackOp = link->UnpackMaxloc;
753   else if (op == MPI_MINLOC) *UnpackOp = link->UnpackMinloc;
754   else *UnpackOp = NULL;
755   PetscFunctionReturn(0);
756 }
757 static PetscErrorCode PetscSFBasicPackGetFetchAndOp(PetscSF sf,PetscSFBasicPack link,MPI_Op op,void (**FetchAndOp)(PetscInt,PetscInt,const PetscInt*,void*,void*))
758 {
759   PetscFunctionBegin;
760   *FetchAndOp = NULL;
761   if (op == MPIU_REPLACE) *FetchAndOp = link->FetchAndInsert;
762   else if (op == MPI_SUM || op == MPIU_SUM) *FetchAndOp = link->FetchAndAdd;
763   else if (op == MPI_MAX || op == MPIU_MAX) *FetchAndOp = link->FetchAndMax;
764   else if (op == MPI_MIN || op == MPIU_MIN) *FetchAndOp = link->FetchAndMin;
765   else if (op == MPI_MAXLOC) *FetchAndOp = link->FetchAndMaxloc;
766   else if (op == MPI_MINLOC) *FetchAndOp = link->FetchAndMinloc;
767   else if (op == MPI_PROD)   *FetchAndOp = link->FetchAndMult;
768   else if (op == MPI_LAND)   *FetchAndOp = link->FetchAndLAND;
769   else if (op == MPI_BAND)   *FetchAndOp = link->FetchAndBAND;
770   else if (op == MPI_LOR)    *FetchAndOp = link->FetchAndLOR;
771   else if (op == MPI_BOR)    *FetchAndOp = link->FetchAndBOR;
772   else if (op == MPI_LXOR)   *FetchAndOp = link->FetchAndLXOR;
773   else if (op == MPI_BXOR)   *FetchAndOp = link->FetchAndBXOR;
774   else SETERRQ(PetscObjectComm((PetscObject)sf),PETSC_ERR_SUP,"No support for MPI_Op");
775   PetscFunctionReturn(0);
776 }
777 
778 typedef enum {PETSC_SF_LEAF2ROOT_REDUCE, PETSC_SF_ROOT2LEAF_BCAST} PetscSFDirection;
779 
780 static PetscErrorCode PetscSFBasicPackGetReqs(PetscSF sf,PetscSFBasicPack link,PetscSFDirection direction,MPI_Request **rootreqs,MPI_Request **leafreqs)
781 {
782   PetscSF_Basic *bas   = (PetscSF_Basic*)sf->data;
783   PetscInt       shift = (direction == PETSC_SF_LEAF2ROOT_REDUCE)? 0 : (sf->nranks + bas->niranks); /* reduce reqs are in the front, bcast reqs are at the end */
784 
785   PetscFunctionBegin;
786   if (rootreqs) *rootreqs = link->requests + shift;
787   if (leafreqs) *leafreqs = link->requests + (bas->niranks - bas->ndiranks) + shift;
788   PetscFunctionReturn(0);
789 }
790 
791 static PetscErrorCode PetscSFBasicPackWaitall(PetscSF sf,PetscSFBasicPack link,PetscSFDirection direction)
792 {
793   PetscSF_Basic  *bas  = (PetscSF_Basic*)sf->data;
794   PetscInt       shift = (direction == PETSC_SF_LEAF2ROOT_REDUCE)? 0 : (sf->nranks + bas->niranks);
795   PetscErrorCode ierr;
796 
797   PetscFunctionBegin;
798   ierr = MPI_Waitall(bas->niranks+sf->nranks-(bas->ndiranks+sf->ndranks),link->requests+shift,MPI_STATUSES_IGNORE);CHKERRQ(ierr);
799   PetscFunctionReturn(0);
800 }
801 
802 static PetscErrorCode PetscSFBasicGetRootInfo(PetscSF sf,PetscInt *nrootranks,PetscInt *ndrootranks,const PetscMPIInt **rootranks,const PetscInt **rootoffset,const PetscInt **rootloc)
803 {
804   PetscSF_Basic *bas = (PetscSF_Basic*)sf->data;
805 
806   PetscFunctionBegin;
807   if (nrootranks)  *nrootranks  = bas->niranks;
808   if (ndrootranks) *ndrootranks = bas->ndiranks;
809   if (rootranks)   *rootranks   = bas->iranks;
810   if (rootoffset)  *rootoffset  = bas->ioffset;
811   if (rootloc)     *rootloc     = bas->irootloc;
812   PetscFunctionReturn(0);
813 }
814 
815 static PetscErrorCode PetscSFBasicGetLeafInfo(PetscSF sf,PetscInt *nleafranks,PetscInt *ndleafranks,const PetscMPIInt **leafranks,const PetscInt **leafoffset,const PetscInt **leafloc)
816 {
817   PetscFunctionBegin;
818   if (nleafranks)  *nleafranks  = sf->nranks;
819   if (ndleafranks) *ndleafranks = sf->ndranks;
820   if (leafranks)   *leafranks   = sf->ranks;
821   if (leafoffset)  *leafoffset  = sf->roffset;
822   if (leafloc)     *leafloc     = sf->rmine;
823   PetscFunctionReturn(0);
824 }
825 
826 static PetscErrorCode PetscSFBasicGetPack(PetscSF sf,MPI_Datatype unit,const void *key,PetscSFBasicPack *mylink)
827 {
828   PetscSF_Basic    *bas = (PetscSF_Basic*)sf->data;
829   PetscErrorCode   ierr;
830   PetscSFBasicPack link,*p;
831   PetscInt         nrootranks,ndrootranks,nleafranks,ndleafranks,i,half;
832   const PetscInt   *rootoffset,*leafoffset;
833   MPI_Comm         comm;
834   PetscMPIInt      n;
835   MPI_Request      *rootreqs,*leafreqs;
836 
837   PetscFunctionBegin;
838   /* Look for types in cache */
839   for (p=&bas->avail; (link=*p); p=&link->next) {
840     PetscBool match;
841     ierr = MPIPetsc_Type_compare(unit,link->unit,&match);CHKERRQ(ierr);
842     if (match) {
843       *p = link->next;          /* Remove from available list */
844       goto found;
845     }
846   }
847 
848   /* Create new composite types for each send rank */
849   ierr = PetscSFBasicGetRootInfo(sf,&nrootranks,&ndrootranks,NULL,&rootoffset,NULL);CHKERRQ(ierr);
850   ierr = PetscSFBasicGetLeafInfo(sf,&nleafranks,&ndleafranks,NULL,&leafoffset,NULL);CHKERRQ(ierr);
851   ierr = PetscNew(&link);CHKERRQ(ierr);
852   ierr = PetscSFBasicPackTypeSetup(link,unit);CHKERRQ(ierr);
853   ierr = PetscMalloc2(nrootranks,&link->root,nleafranks,&link->leaf);CHKERRQ(ierr);
854   /* Double the requests. First half are used for reduce (leaf to root) communication, second half for bcast (root to leaf) communication */
855   half     = nrootranks + nleafranks;
856   ierr     = PetscCalloc1(half*2,&link->requests);CHKERRQ(ierr);
857   rootreqs = link->requests;
858   leafreqs = link->requests + bas->niranks - bas->ndiranks;
859   comm     = PetscObjectComm((PetscObject)sf);
860 
861   /* Allocate buffer and then init the persistent communcation */
862   for (i=0; i<nrootranks; i++) {
863     ierr = PetscMalloc((rootoffset[i+1]-rootoffset[i])*link->unitbytes,&link->root[i]);CHKERRQ(ierr);
864     if (i >= ndrootranks) {
865       ierr = PetscMPIIntCast(rootoffset[i+1]-rootoffset[i],&n);CHKERRQ(ierr);
866       ierr = MPI_Recv_init(link->root[i],n,unit,bas->iranks[i],bas->tag,comm,&rootreqs[i-ndrootranks]);CHKERRQ(ierr);      /* reduce */
867       ierr = MPI_Send_init(link->root[i],n,unit,bas->iranks[i],bas->tag,comm,&rootreqs[i-ndrootranks+half]);CHKERRQ(ierr); /* bcast  */
868     }
869   }
870   for (i=0; i<nleafranks; i++) {
871     if (i < ndleafranks) {      /* Leaf buffers for distinguished ranks are pointers directly into root buffers */
872       if (ndrootranks != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cannot match distinguished ranks");
873       link->leaf[i] = link->root[0];
874       continue;
875     }
876     ierr = PetscMalloc((leafoffset[i+1]-leafoffset[i])*link->unitbytes,&link->leaf[i]);CHKERRQ(ierr);
877     ierr = PetscMPIIntCast(leafoffset[i+1]-leafoffset[i],&n);CHKERRQ(ierr);
878     ierr = MPI_Send_init(link->leaf[i],n,unit,sf->ranks[i],bas->tag,comm,&leafreqs[i-ndleafranks]);CHKERRQ(ierr);      /* reduce */
879     ierr = MPI_Recv_init(link->leaf[i],n,unit,sf->ranks[i],bas->tag,comm,&leafreqs[i-ndleafranks+half]);CHKERRQ(ierr); /* bcast  */
880   }
881 
882 found:
883   link->key  = key;
884   link->next = bas->inuse;
885   bas->inuse = link;
886 
887   *mylink = link;
888   PetscFunctionReturn(0);
889 }
890 
891 static PetscErrorCode PetscSFBasicGetPackInUse(PetscSF sf,MPI_Datatype unit,const void *key,PetscCopyMode cmode,PetscSFBasicPack *mylink)
892 {
893   PetscSF_Basic    *bas = (PetscSF_Basic*)sf->data;
894   PetscErrorCode   ierr;
895   PetscSFBasicPack link,*p;
896 
897   PetscFunctionBegin;
898   /* Look for types in cache */
899   for (p=&bas->inuse; (link=*p); p=&link->next) {
900     PetscBool match;
901     ierr = MPIPetsc_Type_compare(unit,link->unit,&match);CHKERRQ(ierr);
902     if (match && (key == link->key)) {
903       switch (cmode) {
904       case PETSC_OWN_POINTER: *p = link->next; break; /* Remove from inuse list */
905       case PETSC_USE_POINTER: break;
906       default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"invalid cmode");
907       }
908       *mylink = link;
909       PetscFunctionReturn(0);
910     }
911   }
912   SETERRQ(PetscObjectComm((PetscObject)sf),PETSC_ERR_ARG_WRONGSTATE,"Could not find pack");
913   PetscFunctionReturn(0);
914 }
915 
916 static PetscErrorCode PetscSFBasicReclaimPack(PetscSF sf,PetscSFBasicPack *link)
917 {
918   PetscSF_Basic *bas = (PetscSF_Basic*)sf->data;
919 
920   PetscFunctionBegin;
921   (*link)->key  = NULL;
922   (*link)->next = bas->avail;
923   bas->avail    = *link;
924   *link         = NULL;
925   PetscFunctionReturn(0);
926 }
927 
928 static PetscErrorCode PetscSFSetFromOptions_Basic(PetscOptionItems *PetscOptionsObject,PetscSF sf)
929 {
930   PetscErrorCode ierr;
931 
932   PetscFunctionBegin;
933   ierr = PetscOptionsHead(PetscOptionsObject,"PetscSF Basic options");CHKERRQ(ierr);
934   ierr = PetscOptionsTail();CHKERRQ(ierr);
935   PetscFunctionReturn(0);
936 }
937 
938 static PetscErrorCode PetscSFReset_Basic(PetscSF sf)
939 {
940   PetscSF_Basic    *bas = (PetscSF_Basic*)sf->data;
941   PetscErrorCode   ierr;
942   PetscSFBasicPack link,next;
943 
944   PetscFunctionBegin;
945   if (bas->inuse) SETERRQ(PetscObjectComm((PetscObject)sf),PETSC_ERR_ARG_WRONGSTATE,"Outstanding operation has not been completed");
946   ierr = PetscFree2(bas->iranks,bas->ioffset);CHKERRQ(ierr);
947   ierr = PetscFree(bas->irootloc);CHKERRQ(ierr);
948   for (link=bas->avail; link; link=next) {
949     PetscInt i;
950     next = link->next;
951     if (!link->isbuiltin) {ierr = MPI_Type_free(&link->unit);CHKERRQ(ierr);}
952     for (i=0; i<bas->niranks; i++) {ierr = PetscFree(link->root[i]);CHKERRQ(ierr);}
953     for (i=sf->ndranks; i<sf->nranks; i++) {ierr = PetscFree(link->leaf[i]);CHKERRQ(ierr);} /* Free only non-distinguished leaf buffers */
954     ierr = PetscFree2(link->root,link->leaf);CHKERRQ(ierr);
955     /* Free persistent requests using MPI_Request_free */
956     for (i=0; i<sf->nranks+bas->niranks-(sf->ndranks+bas->ndiranks); i++) {
957       ierr = MPI_Request_free(&link->requests[i]);CHKERRQ(ierr); /* used in reduce */
958       ierr = MPI_Request_free(&link->requests[sf->nranks+bas->niranks+i]);CHKERRQ(ierr); /* used in bcast */
959     }
960     ierr = PetscFree(link->requests);CHKERRQ(ierr);
961     ierr = PetscFree(link);CHKERRQ(ierr);
962   }
963   bas->avail = NULL;
964   PetscFunctionReturn(0);
965 }
966 
967 static PetscErrorCode PetscSFDestroy_Basic(PetscSF sf)
968 {
969   PetscErrorCode ierr;
970 
971   PetscFunctionBegin;
972   ierr = PetscSFReset_Basic(sf);CHKERRQ(ierr);
973   ierr = PetscFree(sf->data);CHKERRQ(ierr);
974   PetscFunctionReturn(0);
975 }
976 
977 static PetscErrorCode PetscSFView_Basic(PetscSF sf,PetscViewer viewer)
978 {
979   /* PetscSF_Basic *bas = (PetscSF_Basic*)sf->data; */
980   PetscErrorCode ierr;
981   PetscBool      iascii;
982 
983   PetscFunctionBegin;
984   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
985   if (iascii) {
986     ierr = PetscViewerASCIIPrintf(viewer,"  sort=%s\n",sf->rankorder ? "rank-order" : "unordered");CHKERRQ(ierr);
987   }
988   PetscFunctionReturn(0);
989 }
990 
991 static PetscErrorCode PetscSFBcastAndOpBegin_Basic(PetscSF sf,MPI_Datatype unit,const void *rootdata,void *leafdata,MPI_Op op)
992 {
993   PetscErrorCode    ierr;
994   PetscSFBasicPack  link;
995   PetscInt          i,nrootranks,ndrootranks,nleafranks,ndleafranks;
996   const PetscInt    *rootoffset,*leafoffset,*rootloc,*leafloc;
997   const PetscMPIInt *rootranks,*leafranks;
998   MPI_Request       *rootreqs,*leafreqs;
999   PetscMPIInt       n;
1000 
1001   PetscFunctionBegin;
1002   ierr = PetscSFBasicGetRootInfo(sf,&nrootranks,&ndrootranks,&rootranks,&rootoffset,&rootloc);CHKERRQ(ierr);
1003   ierr = PetscSFBasicGetLeafInfo(sf,&nleafranks,&ndleafranks,&leafranks,&leafoffset,&leafloc);CHKERRQ(ierr);
1004   ierr = PetscSFBasicGetPack(sf,unit,rootdata,&link);CHKERRQ(ierr);
1005 
1006   ierr = PetscSFBasicPackGetReqs(sf,link,PETSC_SF_ROOT2LEAF_BCAST,&rootreqs,&leafreqs);CHKERRQ(ierr);
1007   /* Eagerly post leaf receives, but only from non-distinguished ranks -- distinguished ranks will receive via shared memory */
1008   ierr = PetscMPIIntCast(leafoffset[nleafranks]-leafoffset[ndleafranks],&n);CHKERRQ(ierr);
1009   ierr = MPI_Startall_irecv(n,unit,nleafranks-ndleafranks,leafreqs);CHKERRQ(ierr);
1010 
1011   /* Pack and send root data */
1012   for (i=0; i<nrootranks; i++) {
1013     void *packstart = link->root[i];
1014     ierr = PetscMPIIntCast(rootoffset[i+1]-rootoffset[i],&n);CHKERRQ(ierr);
1015     (*link->Pack)(n,link->bs,rootloc+rootoffset[i],rootdata,packstart);
1016     if (i < ndrootranks) continue; /* shared memory */
1017     ierr = MPI_Start_isend(n,unit,&rootreqs[i-ndrootranks]);CHKERRQ(ierr);
1018   }
1019   PetscFunctionReturn(0);
1020 }
1021 
1022 PetscErrorCode PetscSFBcastAndOpEnd_Basic(PetscSF sf,MPI_Datatype unit,const void *rootdata,void *leafdata,MPI_Op op)
1023 {
1024   PetscErrorCode   ierr;
1025   PetscSFBasicPack link;
1026   PetscInt         i,nleafranks,ndleafranks;
1027   const PetscInt   *leafoffset,*leafloc;
1028   void             (*UnpackOp)(PetscInt,PetscInt,const PetscInt*,void*,const void*);
1029   PetscMPIInt      typesize = -1;
1030 
1031   PetscFunctionBegin;
1032   ierr = PetscSFBasicGetPackInUse(sf,unit,rootdata,PETSC_OWN_POINTER,&link);CHKERRQ(ierr);
1033   ierr = PetscSFBasicPackWaitall(sf,link,PETSC_SF_ROOT2LEAF_BCAST);CHKERRQ(ierr);
1034   ierr = PetscSFBasicGetLeafInfo(sf,&nleafranks,&ndleafranks,NULL,&leafoffset,&leafloc);CHKERRQ(ierr);
1035   ierr = PetscSFBasicPackGetUnpackOp(sf,link,op,&UnpackOp);CHKERRQ(ierr);
1036 
1037   if (UnpackOp) { typesize = link->unitbytes; }
1038   else { ierr = MPI_Type_size(unit,&typesize);CHKERRQ(ierr); }
1039 
1040   for (i=0; i<nleafranks; i++) {
1041     PetscMPIInt n   = leafoffset[i+1] - leafoffset[i];
1042     char *packstart = (char *) link->leaf[i];
1043     if (UnpackOp) { (*UnpackOp)(n,link->bs,leafloc+leafoffset[i],leafdata,(const void *)packstart); }
1044 #if defined(PETSC_HAVE_MPI_REDUCE_LOCAL)
1045     else if (n) { /* the op should be defined to operate on the whole datatype, so we ignore link->bs */
1046       PetscInt j;
1047       for (j=0; j<n; j++) { ierr = MPI_Reduce_local(packstart+j*typesize,((char *) leafdata)+(leafloc[leafoffset[i]+j])*typesize,1,unit,op);CHKERRQ(ierr); }
1048     }
1049 #else
1050     else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No unpacking reduction operation for this MPI_Op");
1051 #endif
1052   }
1053 
1054   ierr = PetscSFBasicReclaimPack(sf,&link);CHKERRQ(ierr);
1055   PetscFunctionReturn(0);
1056 }
1057 
1058 /* Send from roots to leaves */
1059 static PetscErrorCode PetscSFBcastBegin_Basic(PetscSF sf,MPI_Datatype unit,const void *rootdata,void *leafdata)
1060 {
1061   PetscErrorCode   ierr;
1062 
1063   PetscFunctionBegin;
1064   ierr = PetscSFBcastAndOpBegin_Basic(sf,unit,rootdata,leafdata,MPI_REPLACE);CHKERRQ(ierr);
1065   PetscFunctionReturn(0);
1066 }
1067 
1068 PetscErrorCode PetscSFBcastEnd_Basic(PetscSF sf,MPI_Datatype unit,const void *rootdata,void *leafdata)
1069 {
1070   PetscErrorCode   ierr;
1071 
1072   PetscFunctionBegin;
1073   ierr = PetscSFBcastAndOpEnd_Basic(sf,unit,rootdata,leafdata,MPI_REPLACE);CHKERRQ(ierr);
1074   PetscFunctionReturn(0);
1075 }
1076 
1077 /* leaf -> root with reduction */
1078 PetscErrorCode PetscSFReduceBegin_Basic(PetscSF sf,MPI_Datatype unit,const void *leafdata,void *rootdata,MPI_Op op)
1079 {
1080   PetscSFBasicPack  link;
1081   PetscErrorCode    ierr;
1082   PetscInt          i,nrootranks,ndrootranks,nleafranks,ndleafranks;
1083   const PetscInt    *rootoffset,*leafoffset,*rootloc,*leafloc;
1084   const PetscMPIInt *rootranks,*leafranks;
1085   MPI_Request       *rootreqs,*leafreqs;
1086   PetscMPIInt       n;
1087 
1088   PetscFunctionBegin;
1089   ierr = PetscSFBasicGetRootInfo(sf,&nrootranks,&ndrootranks,&rootranks,&rootoffset,&rootloc);CHKERRQ(ierr);
1090   ierr = PetscSFBasicGetLeafInfo(sf,&nleafranks,&ndleafranks,&leafranks,&leafoffset,&leafloc);CHKERRQ(ierr);
1091   ierr = PetscSFBasicGetPack(sf,unit,leafdata,&link);CHKERRQ(ierr);
1092 
1093   ierr = PetscSFBasicPackGetReqs(sf,link,PETSC_SF_LEAF2ROOT_REDUCE,&rootreqs,&leafreqs);CHKERRQ(ierr);
1094   /* Eagerly post root receives for non-distinguished ranks */
1095   ierr = PetscMPIIntCast(rootoffset[nrootranks]-rootoffset[ndrootranks],&n);CHKERRQ(ierr);
1096   ierr = MPI_Startall_irecv(n,unit,nrootranks-ndrootranks,rootreqs);CHKERRQ(ierr);
1097 
1098   /* Pack and send leaf data */
1099   for (i=0; i<nleafranks; i++) {
1100     void *packstart = link->leaf[i];
1101     ierr = PetscMPIIntCast(leafoffset[i+1]-leafoffset[i],&n);CHKERRQ(ierr);
1102     (*link->Pack)(n,link->bs,leafloc+leafoffset[i],leafdata,packstart);
1103     if (i < ndleafranks) continue; /* shared memory */
1104     ierr = MPI_Start_isend(n,unit,&leafreqs[i-ndleafranks]);CHKERRQ(ierr);
1105   }
1106   PetscFunctionReturn(0);
1107 }
1108 
1109 static PetscErrorCode PetscSFReduceEnd_Basic(PetscSF sf,MPI_Datatype unit,const void *leafdata,void *rootdata,MPI_Op op)
1110 {
1111   void             (*UnpackOp)(PetscInt,PetscInt,const PetscInt*,void*,const void*);
1112   PetscErrorCode   ierr;
1113   PetscSFBasicPack link;
1114   PetscInt         i,nrootranks;
1115   PetscMPIInt      typesize = -1;
1116   const PetscInt   *rootoffset,*rootloc;
1117 
1118   PetscFunctionBegin;
1119   ierr = PetscSFBasicGetPackInUse(sf,unit,leafdata,PETSC_OWN_POINTER,&link);CHKERRQ(ierr);
1120   /* This implementation could be changed to unpack as receives arrive, at the cost of non-determinism */
1121   ierr = PetscSFBasicPackWaitall(sf,link,PETSC_SF_LEAF2ROOT_REDUCE);CHKERRQ(ierr);
1122   ierr = PetscSFBasicGetRootInfo(sf,&nrootranks,NULL,NULL,&rootoffset,&rootloc);CHKERRQ(ierr);
1123   ierr = PetscSFBasicPackGetUnpackOp(sf,link,op,&UnpackOp);CHKERRQ(ierr);
1124   if (UnpackOp) {
1125     typesize = link->unitbytes;
1126   }
1127   else {
1128     ierr = MPI_Type_size(unit,&typesize);CHKERRQ(ierr);
1129   }
1130   for (i=0; i<nrootranks; i++) {
1131     PetscMPIInt n   = rootoffset[i+1] - rootoffset[i];
1132     char *packstart = (char *) link->root[i];
1133 
1134     if (UnpackOp) {
1135       (*UnpackOp)(n,link->bs,rootloc+rootoffset[i],rootdata,(const void *)packstart);
1136     }
1137 #if defined(PETSC_HAVE_MPI_REDUCE_LOCAL)
1138     else if (n) { /* the op should be defined to operate on the whole datatype, so we ignore link->bs */
1139       PetscInt j;
1140 
1141       for (j = 0; j < n; j++) {
1142         ierr = MPI_Reduce_local(packstart+j*typesize,((char *) rootdata)+(rootloc[rootoffset[i]+j])*typesize,1,unit,op);CHKERRQ(ierr);
1143       }
1144     }
1145 #else
1146     else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No unpacking reduction operation for this MPI_Op");
1147 #endif
1148   }
1149   ierr = PetscSFBasicReclaimPack(sf,&link);CHKERRQ(ierr);
1150   PetscFunctionReturn(0);
1151 }
1152 
1153 static PetscErrorCode PetscSFFetchAndOpBegin_Basic(PetscSF sf,MPI_Datatype unit,void *rootdata,const void *leafdata,void *leafupdate,MPI_Op op)
1154 {
1155   PetscErrorCode ierr;
1156 
1157   PetscFunctionBegin;
1158   ierr = PetscSFReduceBegin_Basic(sf,unit,leafdata,rootdata,op);CHKERRQ(ierr);
1159   PetscFunctionReturn(0);
1160 }
1161 
1162 static PetscErrorCode PetscSFFetchAndOpEnd_Basic(PetscSF sf,MPI_Datatype unit,void *rootdata,const void *leafdata,void *leafupdate,MPI_Op op)
1163 {
1164   void              (*FetchAndOp)(PetscInt,PetscInt,const PetscInt*,void*,void*);
1165   PetscErrorCode    ierr;
1166   PetscSFBasicPack  link;
1167   PetscInt          i,nrootranks,ndrootranks,nleafranks,ndleafranks;
1168   const PetscInt    *rootoffset,*leafoffset,*rootloc,*leafloc;
1169   const PetscMPIInt *rootranks,*leafranks;
1170   MPI_Request       *rootreqs,*leafreqs;
1171   PetscMPIInt       n;
1172 
1173   PetscFunctionBegin;
1174   ierr = PetscSFBasicGetPackInUse(sf,unit,leafdata,PETSC_OWN_POINTER,&link);CHKERRQ(ierr);
1175   /* This implementation could be changed to unpack as receives arrive, at the cost of non-determinism */
1176   ierr = PetscSFBasicPackWaitall(sf,link,PETSC_SF_LEAF2ROOT_REDUCE);CHKERRQ(ierr);
1177   ierr = PetscSFBasicGetRootInfo(sf,&nrootranks,&ndrootranks,&rootranks,&rootoffset,&rootloc);CHKERRQ(ierr);
1178   ierr = PetscSFBasicGetLeafInfo(sf,&nleafranks,&ndleafranks,&leafranks,&leafoffset,&leafloc);CHKERRQ(ierr);
1179   ierr = PetscSFBasicPackGetReqs(sf,link,PETSC_SF_ROOT2LEAF_BCAST,&rootreqs,&leafreqs);CHKERRQ(ierr);
1180   /* Post leaf receives */
1181   ierr = PetscMPIIntCast(leafoffset[nleafranks]-leafoffset[ndleafranks],&n);CHKERRQ(ierr);
1182   ierr = MPI_Startall_irecv(n,unit,nleafranks-ndleafranks,leafreqs);CHKERRQ(ierr);
1183 
1184   /* Process local fetch-and-op, post root sends */
1185   ierr = PetscSFBasicPackGetFetchAndOp(sf,link,op,&FetchAndOp);CHKERRQ(ierr);
1186   for (i=0; i<nrootranks; i++) {
1187     void *packstart = link->root[i];
1188     ierr = PetscMPIIntCast(rootoffset[i+1]-rootoffset[i],&n);CHKERRQ(ierr);
1189     (*FetchAndOp)(n,link->bs,rootloc+rootoffset[i],rootdata,packstart);
1190     if (i < ndrootranks) continue; /* shared memory */
1191     ierr = MPI_Start_isend(n,unit,&rootreqs[i-ndrootranks]);CHKERRQ(ierr);
1192   }
1193   ierr = PetscSFBasicPackWaitall(sf,link,PETSC_SF_ROOT2LEAF_BCAST);CHKERRQ(ierr);
1194   for (i=0; i<nleafranks; i++) {
1195     const void  *packstart = link->leaf[i];
1196     ierr = PetscMPIIntCast(leafoffset[i+1]-leafoffset[i],&n);CHKERRQ(ierr);
1197     (*link->UnpackInsert)(n,link->bs,leafloc+leafoffset[i],leafupdate,packstart);
1198   }
1199   ierr = PetscSFBasicReclaimPack(sf,&link);CHKERRQ(ierr);
1200   PetscFunctionReturn(0);
1201 }
1202 
1203 static PetscErrorCode PetscSFGetLeafRanks_Basic(PetscSF sf,PetscInt *niranks,const PetscMPIInt **iranks,const PetscInt **ioffset,const PetscInt **irootloc)
1204 {
1205   PetscSF_Basic *bas = (PetscSF_Basic*)sf->data;
1206 
1207   PetscFunctionBegin;
1208   if (niranks)  *niranks  = bas->niranks;
1209   if (iranks)   *iranks   = bas->iranks;
1210   if (ioffset)  *ioffset  = bas->ioffset;
1211   if (irootloc) *irootloc = bas->irootloc;
1212   PetscFunctionReturn(0);
1213 }
1214 
1215 PETSC_EXTERN PetscErrorCode PetscSFCreate_Basic(PetscSF sf)
1216 {
1217   PetscSF_Basic  *bas = (PetscSF_Basic*)sf->data;
1218   PetscErrorCode ierr;
1219 
1220   PetscFunctionBegin;
1221   sf->ops->SetUp           = PetscSFSetUp_Basic;
1222   sf->ops->SetFromOptions  = PetscSFSetFromOptions_Basic;
1223   sf->ops->Reset           = PetscSFReset_Basic;
1224   sf->ops->Destroy         = PetscSFDestroy_Basic;
1225   sf->ops->View            = PetscSFView_Basic;
1226   sf->ops->BcastBegin      = PetscSFBcastBegin_Basic;
1227   sf->ops->BcastEnd        = PetscSFBcastEnd_Basic;
1228   sf->ops->BcastAndOpBegin = PetscSFBcastAndOpBegin_Basic;
1229   sf->ops->BcastAndOpEnd   = PetscSFBcastAndOpEnd_Basic;
1230   sf->ops->ReduceBegin     = PetscSFReduceBegin_Basic;
1231   sf->ops->ReduceEnd       = PetscSFReduceEnd_Basic;
1232   sf->ops->FetchAndOpBegin = PetscSFFetchAndOpBegin_Basic;
1233   sf->ops->FetchAndOpEnd   = PetscSFFetchAndOpEnd_Basic;
1234   sf->ops->GetLeafRanks    = PetscSFGetLeafRanks_Basic;
1235 
1236   ierr = PetscNewLog(sf,&bas);CHKERRQ(ierr);
1237   sf->data = (void*)bas;
1238   PetscFunctionReturn(0);
1239 }
1240