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