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