Actual source code: vscat.c
1: #include <petsc/private/sfimpl.h>
2: #include <../src/vec/is/sf/impls/basic/sfbasic.h>
3: #include <../src/vec/is/sf/impls/basic/sfpack.h>
4: #include <petsc/private/vecimpl.h>
6: PETSC_SINGLE_LIBRARY_INTERN PetscErrorCode PetscGatherNumberOfMessages_Private(MPI_Comm, const PetscMPIInt[], const PetscInt[], PetscMPIInt *);
7: PETSC_SINGLE_LIBRARY_INTERN PetscErrorCode PetscGatherMessageLengths_Private(MPI_Comm, PetscMPIInt, PetscMPIInt, const PetscInt[], PetscMPIInt **, PetscInt **);
9: typedef enum {
10: IS_INVALID,
11: IS_GENERAL,
12: IS_BLOCK,
13: IS_STRIDE
14: } ISTypeID;
16: static inline PetscErrorCode ISGetTypeID_Private(IS is, ISTypeID *id)
17: {
18: PetscBool same;
20: PetscFunctionBegin;
21: *id = IS_INVALID;
22: PetscCall(PetscObjectTypeCompare((PetscObject)is, ISGENERAL, &same));
23: if (same) {
24: *id = IS_GENERAL;
25: goto functionend;
26: }
27: PetscCall(PetscObjectTypeCompare((PetscObject)is, ISBLOCK, &same));
28: if (same) {
29: *id = IS_BLOCK;
30: goto functionend;
31: }
32: PetscCall(PetscObjectTypeCompare((PetscObject)is, ISSTRIDE, &same));
33: if (same) {
34: *id = IS_STRIDE;
35: goto functionend;
36: }
37: functionend:
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: static PetscErrorCode VecScatterBegin_Internal(VecScatter sf, Vec x, Vec y, InsertMode addv, ScatterMode mode)
42: {
43: PetscSF wsf = NULL; /* either sf or its local part */
44: MPI_Op mop = MPI_OP_NULL;
45: PetscMPIInt size;
46: PetscMemType xmtype = PETSC_MEMTYPE_HOST, ymtype = PETSC_MEMTYPE_HOST;
48: PetscFunctionBegin;
49: if (x != y) PetscCall(VecLockReadPush(x));
50: PetscCall(VecGetArrayReadAndMemType(x, &sf->vscat.xdata, &xmtype));
51: PetscCall(VecGetArrayAndMemType(y, &sf->vscat.ydata, &ymtype));
52: PetscCall(VecLockWriteSet(y, PETSC_TRUE));
54: /* SCATTER_FORWARD_LOCAL indicates ignoring inter-process communication */
55: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size));
56: if ((mode & SCATTER_FORWARD_LOCAL) && size > 1) { /* Lazy creation of sf->vscat.lsf since SCATTER_FORWARD_LOCAL is uncommon */
57: if (!sf->vscat.lsf) PetscCall(PetscSFCreateLocalSF_Private(sf, &sf->vscat.lsf));
58: wsf = sf->vscat.lsf;
59: } else {
60: wsf = sf;
61: }
63: /* Note xdata/ydata is always recorded on sf (not lsf) above */
64: if (addv == INSERT_VALUES) mop = MPI_REPLACE;
65: else if (addv == ADD_VALUES) mop = MPIU_SUM; /* PETSc defines its own MPI datatype and SUM operation for __float128 etc. */
66: else if (addv == MAX_VALUES) mop = MPIU_MAX;
67: else if (addv == MIN_VALUES) mop = MPIU_MIN;
68: else SETERRQ(PetscObjectComm((PetscObject)sf), PETSC_ERR_SUP, "Unsupported InsertMode %d in VecScatterBegin/End", addv);
70: if (mode & SCATTER_REVERSE) { /* REVERSE indicates leaves to root scatter. Note that x and y are swapped in input */
71: PetscCall(PetscSFReduceWithMemTypeBegin(wsf, sf->vscat.unit, xmtype, sf->vscat.xdata, ymtype, sf->vscat.ydata, mop));
72: } else { /* FORWARD indicates x to y scatter, where x is root and y is leaf */
73: PetscCall(PetscSFBcastWithMemTypeBegin(wsf, sf->vscat.unit, xmtype, sf->vscat.xdata, ymtype, sf->vscat.ydata, mop));
74: }
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: static PetscErrorCode VecScatterEnd_Internal(VecScatter sf, Vec x, Vec y, InsertMode addv, ScatterMode mode)
79: {
80: PetscSF wsf = NULL;
81: MPI_Op mop = MPI_OP_NULL;
82: PetscMPIInt size;
84: PetscFunctionBegin;
85: /* SCATTER_FORWARD_LOCAL indicates ignoring inter-process communication */
86: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size));
87: wsf = ((mode & SCATTER_FORWARD_LOCAL) && size > 1) ? sf->vscat.lsf : sf;
89: if (addv == INSERT_VALUES) mop = MPI_REPLACE;
90: else if (addv == ADD_VALUES) mop = MPIU_SUM;
91: else if (addv == MAX_VALUES) mop = MPIU_MAX;
92: else if (addv == MIN_VALUES) mop = MPIU_MIN;
93: else SETERRQ(PetscObjectComm((PetscObject)sf), PETSC_ERR_SUP, "Unsupported InsertMode %d in VecScatterBegin/End", addv);
95: if (mode & SCATTER_REVERSE) { /* reverse scatter sends leaves to roots. Note that x and y are swapped in input */
96: PetscCall(PetscSFReduceEnd(wsf, sf->vscat.unit, sf->vscat.xdata, sf->vscat.ydata, mop));
97: } else { /* forward scatter sends roots to leaves, i.e., x to y */
98: PetscCall(PetscSFBcastEnd(wsf, sf->vscat.unit, sf->vscat.xdata, sf->vscat.ydata, mop));
99: }
101: PetscCall(VecRestoreArrayReadAndMemType(x, &sf->vscat.xdata));
102: if (x != y) PetscCall(VecLockReadPop(x));
103: PetscCall(VecRestoreArrayAndMemType(y, &sf->vscat.ydata));
104: PetscCall(VecLockWriteSet(y, PETSC_FALSE));
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: /* VecScatterRemap provides a light way to slightly modify a VecScatter. Suppose the input sf scatters
109: x[i] to y[j], tomap gives a plan to change vscat to scatter x[tomap[i]] to y[j]. Note that in SF,
110: x is roots. That means we need to change incoming stuffs such as bas->irootloc[].
111: */
112: static PetscErrorCode VecScatterRemap_Internal(VecScatter sf, const PetscInt *tomap, const PetscInt *frommap)
113: {
114: PetscInt i, bs = sf->vscat.bs;
115: PetscMPIInt size;
116: PetscBool ident = PETSC_TRUE, isbasic, isneighbor;
117: PetscSFType type;
118: PetscSF_Basic *bas = NULL;
120: PetscFunctionBegin;
121: /* check if it is an identity map. If it is, do nothing */
122: if (tomap) {
123: for (i = 0; i < sf->nroots * bs; i++) {
124: if (i != tomap[i]) {
125: ident = PETSC_FALSE;
126: break;
127: }
128: }
129: if (ident) PetscFunctionReturn(PETSC_SUCCESS);
130: }
131: PetscCheck(!frommap, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unable to remap the FROM in scatters yet");
132: if (!tomap) PetscFunctionReturn(PETSC_SUCCESS);
134: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sf), &size));
136: /* Since the indices changed, we must also update the local SF. But we do not do it since
137: lsf is rarely used. We just destroy lsf and rebuild it on demand from updated sf.
138: */
139: PetscCall(PetscSFDestroy(&sf->vscat.lsf));
141: PetscCall(PetscSFGetType(sf, &type));
142: PetscCall(PetscObjectTypeCompare((PetscObject)sf, PETSCSFBASIC, &isbasic));
143: PetscCall(PetscObjectTypeCompare((PetscObject)sf, PETSCSFNEIGHBOR, &isneighbor));
144: PetscCheck(isbasic || isneighbor, PetscObjectComm((PetscObject)sf), PETSC_ERR_SUP, "VecScatterRemap on SF type %s is not supported", type);
146: PetscCall(PetscSFSetUp(sf)); /* to build sf->irootloc if SetUp is not yet called */
148: /* Root indices are going to be remapped. This is tricky for SF. Root indices are used in sf->rremote,
149: sf->remote and bas->irootloc. The latter one is cheap to remap, but the former two are not.
150: To remap them, we have to do a bcast from roots to leaves, to let leaves know their updated roots.
151: Since VecScatterRemap is supposed to be a cheap routine to adapt a vecscatter by only changing where
152: x[] data is taken, we do not remap sf->rremote, sf->remote. The consequence is that operations
153: accessing them (such as PetscSFCompose) may get stale info. Considering VecScatter does not need
154: that complicated SF operations, we do not remap sf->rremote, sf->remote, instead we destroy them
155: so that code accessing them (if any) will crash (instead of get silent errors). Note that BcastAndOp/Reduce,
156: which are used by VecScatter and only rely on bas->irootloc, are updated and correct.
157: */
158: sf->remote = NULL;
159: PetscCall(PetscFree(sf->remote_alloc));
160: /* Not easy to free sf->rremote since it was allocated with PetscMalloc4(), so just give it crazy values */
161: for (i = 0; i < sf->roffset[sf->nranks]; i++) sf->rremote[i] = PETSC_INT_MIN;
163: /* Indices in tomap[] are for each individual vector entry. But indices in sf are for each
164: block in the vector. So before the remapping, we have to expand indices in sf by bs, and
165: after the remapping, we have to shrink them back.
166: */
167: bas = (PetscSF_Basic *)sf->data;
168: for (i = 0; i < bas->ioffset[bas->niranks]; i++) bas->irootloc[i] = tomap[bas->irootloc[i] * bs] / bs;
169: #if PetscDefined(HAVE_DEVICE)
170: /* Free the irootloc copy on device. We allocate a new copy and get the updated value on demand. See PetscSFLinkGetRootPackOptAndIndices() */
171: for (i = 0; i < 2; i++) PetscCall(PetscSFFree(sf, PETSC_MEMTYPE_DEVICE, bas->irootloc_d[i]));
172: #endif
173: /* Destroy and then rebuild root packing optimizations since indices are changed */
174: PetscCall(PetscSFResetPackFields(sf));
175: PetscCall(PetscSFSetUpPackFields(sf));
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: /*
180: Given a parallel VecScatter context, return number of procs and vector entries involved in remote (i.e., off-process) communication
182: Input Parameters:
183: + sf - the context (must be a parallel vecscatter)
184: - send - true to select the send info (i.e., todata), otherwise to select the recv info (i.e., fromdata)
186: Output parameters:
187: + num_procs - number of remote processors
188: - num_entries - number of vector entries to send or recv
190: Notes:
191: Sometimes PETSc internally needs to use the matrix-vector-multiply vecscatter context for other purposes. The client code
192: usually only uses MPI_Send/Recv. This group of subroutines provides info needed for such uses.
194: .seealso: [](sec_scatter), `VecScatterGetRemote_Private()`, `VecScatterGetRemoteOrdered_Private()`
195: */
196: PetscErrorCode VecScatterGetRemoteCount_Private(VecScatter sf, PetscBool send, PetscInt *num_procs, PetscInt *num_entries)
197: {
198: PetscMPIInt nranks, remote_start;
199: PetscMPIInt rank;
200: const PetscInt *offset;
201: const PetscMPIInt *ranks;
203: PetscFunctionBegin;
204: PetscCall(PetscSFSetUp(sf));
205: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sf), &rank));
207: /* This routine is mainly used for MatMult's Mvctx. In Mvctx, we scatter an MPI vector x to a sequential vector lvec.
208: Remember x is roots and lvec is leaves. 'send' means roots to leaves communication. If 'send' is true, we need to
209: get info about which ranks this processor needs to send to. In other words, we need to call PetscSFGetLeafRanks().
210: If send is false, we do the opposite, calling PetscSFGetRootRanks().
211: */
212: if (send) PetscCall(PetscSFGetLeafRanks(sf, &nranks, &ranks, &offset, NULL));
213: else PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &offset, NULL, NULL));
214: if (nranks) {
215: remote_start = (rank == ranks[0]) ? 1 : 0;
216: if (num_procs) *num_procs = nranks - remote_start;
217: if (num_entries) *num_entries = offset[nranks] - offset[remote_start];
218: } else {
219: if (num_procs) *num_procs = 0;
220: if (num_entries) *num_entries = 0;
221: }
222: PetscFunctionReturn(PETSC_SUCCESS);
223: }
225: /* Given a parallel VecScatter context, return a plan that represents the remote communication.
226: Any output parameter can be NULL.
228: Input Parameters:
229: + sf - the context
230: - send - true to select the send info (i.e., todata), otherwise to select the recv info (i.e., fromdata)
232: Output parameters:
233: + n - number of remote processors
234: . starts - starting point in indices for each proc. ATTENTION: starts[0] is not necessarily zero.
235: Therefore, expressions like starts[i+1]-starts[i] and indices[starts[i]+j] work as
236: expected for a CSR structure but buf[starts[i]+j] may be out of range if buf was allocated
237: with length starts[n]-starts[0]. One should use buf[starts[i]-starts[0]+j] instead.
238: . indices - indices of entries to send/recv
239: . procs - ranks of remote processors
240: - bs - block size
242: .seealso: `VecScatterRestoreRemote_Private()`, `VecScatterGetRemoteOrdered_Private()`
243: */
244: PetscErrorCode VecScatterGetRemote_Private(VecScatter sf, PetscBool send, PetscMPIInt *n, const PetscInt **starts, const PetscInt **indices, const PetscMPIInt **procs, PetscInt *bs)
245: {
246: PetscMPIInt nranks, remote_start;
247: PetscMPIInt rank;
248: const PetscInt *offset, *location;
249: const PetscMPIInt *ranks;
251: PetscFunctionBegin;
252: PetscCall(PetscSFSetUp(sf));
253: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sf), &rank));
255: if (send) PetscCall(PetscSFGetLeafRanks(sf, &nranks, &ranks, &offset, &location));
256: else PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &offset, &location, NULL));
258: if (nranks) {
259: remote_start = (rank == ranks[0]) ? 1 : 0;
260: if (n) *n = nranks - remote_start;
261: if (starts) *starts = &offset[remote_start];
262: if (indices) *indices = location; /* not &location[offset[remote_start]]. Starts[0] may point to the middle of indices[] */
263: if (procs) *procs = &ranks[remote_start];
264: } else {
265: if (n) *n = 0;
266: if (starts) *starts = NULL;
267: if (indices) *indices = NULL;
268: if (procs) *procs = NULL;
269: }
271: if (bs) *bs = 1;
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: /* Given a parallel VecScatter context, return a plan that represents the remote communication. Ranks of remote
276: processors returned in procs must be sorted in ascending order. Any output parameter can be NULL.
278: Input Parameters:
279: + sf - the context
280: - send - true to select the send info (i.e., todata), otherwise to select the recv info (i.e., fromdata)
282: Output parameters:
283: + n - number of remote processors
284: . starts - starting point in indices for each proc. ATTENTION: starts[0] is not necessarily zero.
285: Therefore, expressions like starts[i+1]-starts[i] and indices[starts[i]+j] work as
286: expected for a CSR structure but buf[starts[i]+j] may be out of range if buf was allocated
287: with length starts[n]-starts[0]. One should use buf[starts[i]-starts[0]+j] instead.
288: . indices - indices of entries to send/recv
289: . procs - ranks of remote processors
290: - bs - block size
292: Notes:
293: Output parameters like starts, indices must also be adapted according to the sorted ranks.
295: .seealso: `VecScatterRestoreRemoteOrdered_Private()`, `VecScatterGetRemote_Private()`
296: */
297: PetscErrorCode VecScatterGetRemoteOrdered_Private(VecScatter sf, PetscBool send, PetscMPIInt *n, const PetscInt **starts, const PetscInt **indices, const PetscMPIInt **procs, PetscInt *bs)
298: {
299: PetscFunctionBegin;
300: PetscCall(VecScatterGetRemote_Private(sf, send, n, starts, indices, procs, bs));
301: if (PetscUnlikelyDebug(n && procs)) {
302: PetscMPIInt i;
303: /* from back to front to also handle cases *n=0 */
304: for (i = *n - 1; i > 0; i--) PetscCheck((*procs)[i - 1] <= (*procs)[i], PETSC_COMM_SELF, PETSC_ERR_PLIB, "procs[] are not ordered");
305: }
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: /* Given a parallel VecScatter context, restore the plan returned by VecScatterGetRemote_Private. This gives a chance for
310: an implementation to free memory allocated in the VecScatterGetRemote_Private call.
312: Input Parameters:
313: + sf - the context
314: - send - true to select the send info (i.e., todata), otherwise to select the recv info (i.e., fromdata)
316: Output parameters:
317: + n - number of remote processors
318: . starts - starting point in indices for each proc
319: . indices - indices of entries to send/recv
320: . procs - ranks of remote processors
321: - bs - block size
323: .seealso: `VecScatterGetRemote_Private()`
324: */
325: PetscErrorCode VecScatterRestoreRemote_Private(VecScatter sf, PetscBool send, PetscMPIInt *n, const PetscInt **starts, const PetscInt **indices, const PetscMPIInt **procs, PetscInt *bs)
326: {
327: PetscFunctionBegin;
328: if (starts) *starts = NULL;
329: if (indices) *indices = NULL;
330: if (procs) *procs = NULL;
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
334: /* Given a parallel VecScatter context, restore the plan returned by VecScatterGetRemoteOrdered_Private. This gives a chance for
335: an implementation to free memory allocated in the VecScatterGetRemoteOrdered_Private call.
337: Input Parameters:
338: + sf - the context
339: - send - true to select the send info (i.e., todata), otherwise to select the recv info (i.e., fromdata)
341: Output parameters:
342: + n - number of remote processors
343: . starts - starting point in indices for each proc
344: . indices - indices of entries to send/recv
345: . procs - ranks of remote processors
346: - bs - block size
348: .seealso: `VecScatterGetRemoteOrdered_Private()`
349: */
350: PetscErrorCode VecScatterRestoreRemoteOrdered_Private(VecScatter sf, PetscBool send, PetscMPIInt *n, const PetscInt **starts, const PetscInt **indices, const PetscMPIInt **procs, PetscInt *bs)
351: {
352: PetscFunctionBegin;
353: PetscCall(VecScatterRestoreRemote_Private(sf, send, n, starts, indices, procs, bs));
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: /*@
358: VecScatterSetUp - Sets up the `VecScatter` to be able to actually scatter information between vectors
360: Collective
362: Input Parameter:
363: . sf - the scatter context
365: Level: intermediate
367: .seealso: [](sec_scatter), `VecScatter`, `VecScatterCreate()`, `VecScatterCopy()`
368: @*/
369: PetscErrorCode VecScatterSetUp(VecScatter sf)
370: {
371: PetscFunctionBegin;
372: PetscCall(PetscSFSetUp(sf));
373: PetscFunctionReturn(PETSC_SUCCESS);
374: }
376: /*@
377: VecScatterSetType - Builds a vector scatter, for a particular vector scatter implementation.
379: Collective
381: Input Parameters:
382: + sf - The `VecScatter` object
383: - type - The name of the vector scatter type
385: Options Database Key:
386: . -sf_type type - Sets the `VecScatterType`
388: Level: intermediate
390: Note:
391: Use `VecScatterDuplicate()` to form additional vectors scatter of the same type as an existing vector scatter.
393: .seealso: [](sec_scatter), `VecScatter`, `VecScatterType`, `VecScatterGetType()`, `VecScatterCreate()`
394: @*/
395: PetscErrorCode VecScatterSetType(VecScatter sf, VecScatterType type)
396: {
397: PetscFunctionBegin;
398: PetscCall(PetscSFSetType(sf, type));
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: /*@
403: VecScatterGetType - Gets the vector scatter type name (as a string) from the `VecScatter`.
405: Not Collective
407: Input Parameter:
408: . sf - The vector scatter
410: Output Parameter:
411: . type - The vector scatter type name
413: Level: intermediate
415: .seealso: [](sec_scatter), `VecScatter`, `VecScatterType`, `VecScatterSetType()`, `VecScatterCreate()`
416: @*/
417: PetscErrorCode VecScatterGetType(VecScatter sf, VecScatterType *type)
418: {
419: PetscFunctionBegin;
420: PetscCall(PetscSFGetType(sf, type));
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: /*@
425: VecScatterRegister - Adds a new vector scatter component implementation
427: Not Collective
429: Input Parameters:
430: + sname - The name of a new user-defined creation routine
431: - function - The creation routine
433: Level: advanced
435: .seealso: [](sec_scatter), `VecScatter`, `VecScatterType`, `VecRegister()`
436: @*/
437: PetscErrorCode VecScatterRegister(const char sname[], PetscErrorCode (*function)(VecScatter))
438: {
439: PetscFunctionBegin;
440: PetscCall(PetscSFRegister(sname, function));
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
444: /*@
445: VecScatterGetMerged - Returns true if the scatter is completed in the `VecScatterBegin()`
446: and the `VecScatterEnd()` does nothing
448: Not Collective
450: Input Parameter:
451: . sf - scatter context created with `VecScatterCreate()`
453: Output Parameter:
454: . flg - `PETSC_TRUE` if the `VecScatterBegin()`/`VecScatterEnd()` are all done during the `VecScatterBegin()`
456: Level: developer
458: .seealso: [](sec_scatter), `VecScatter`, `VecScatterCreate()`, `VecScatterEnd()`, `VecScatterBegin()`
459: @*/
460: PetscErrorCode VecScatterGetMerged(VecScatter sf, PetscBool *flg)
461: {
462: PetscFunctionBegin;
464: if (flg) *flg = sf->vscat.beginandendtogether;
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
467: /*@
468: VecScatterDestroy - Destroys a scatter context created by `VecScatterCreate()`
470: Collective
472: Input Parameter:
473: . sf - the scatter context
475: Level: intermediate
477: .seealso: [](sec_scatter), `VecScatter`, `VecScatterCreate()`, `VecScatterCopy()`
478: @*/
479: PetscErrorCode VecScatterDestroy(VecScatter *sf)
480: {
481: PetscFunctionBegin;
482: PetscCall(PetscSFDestroy(sf));
483: PetscFunctionReturn(PETSC_SUCCESS);
484: }
486: /*@
487: VecScatterCopy - Makes a copy of a scatter context.
489: Collective
491: Input Parameter:
492: . sf - the scatter context
494: Output Parameter:
495: . newsf - the context copy
497: Level: advanced
499: .seealso: [](sec_scatter), `VecScatter`, `VecScatterType`, `VecScatterCreate()`, `VecScatterDestroy()`
500: @*/
501: PetscErrorCode VecScatterCopy(VecScatter sf, VecScatter *newsf)
502: {
503: PetscFunctionBegin;
504: PetscAssertPointer(newsf, 2);
505: PetscCall(PetscSFDuplicate(sf, PETSCSF_DUPLICATE_GRAPH, newsf));
506: PetscCall(PetscSFSetUp(*newsf));
507: PetscFunctionReturn(PETSC_SUCCESS);
508: }
510: /*@
511: VecScatterViewFromOptions - View a `VecScatter` object based on values in the options database
513: Collective
515: Input Parameters:
516: + sf - the scatter context
517: . obj - optional object that provides the prefix for the option names, pass `NULL` to use the options prefix of `sf`
518: - name - command line option
520: Options Database Key:
521: . -name viewer_specification - See `PetscOptionsCreateViewer()` for the values of `viewer_specification`
523: Level: intermediate
525: Note:
526: This checks the options database, creates the viewer on-the-fly, uses it and then destroys it. Hence it should not be called in heavily used routines,
527: rather `PetscOptionsCreateViewer()` should be used to construct the viewer once which can then be utilized in the heavily used routine.
529: .seealso: [](sec_scatter), `VecScatter`, `VecScatterView()`, `PetscObjectViewFromOptions()`, `VecScatterCreate()`, `PetscOptionsCreateViewer()`
530: @*/
531: PetscErrorCode VecScatterViewFromOptions(VecScatter sf, PetscObject obj, const char name[])
532: {
533: PetscFunctionBegin;
535: PetscCall(PetscObjectViewFromOptions((PetscObject)sf, obj, name));
536: PetscFunctionReturn(PETSC_SUCCESS);
537: }
539: /*@
540: VecScatterView - Views a vector scatter context.
542: Collective
544: Input Parameters:
545: + sf - the scatter context
546: - viewer - the viewer for displaying the context
548: Level: intermediate
550: .seealso: [](sec_scatter), `VecScatter`, `PetscViewer`, `VecScatterViewFromOptions()`, `PetscObjectViewFromOptions()`, `VecScatterCreate()`
551: @*/
552: PetscErrorCode VecScatterView(VecScatter sf, PetscViewer viewer)
553: {
554: PetscFunctionBegin;
555: PetscCall(PetscSFView(sf, viewer));
556: PetscFunctionReturn(PETSC_SUCCESS);
557: }
559: /*@
560: VecScatterRemap - Remaps the "from" and "to" indices in a
561: vector scatter context.
563: Collective
565: Input Parameters:
566: + sf - vector scatter context
567: . tomap - remapping plan for "to" indices (may be `NULL`).
568: - frommap - remapping plan for "from" indices (may be `NULL`)
570: Level: developer
572: Notes:
573: In the parallel case the todata contains indices from where the data is taken
574: (and then sent to others)! The fromdata contains indices from where the received
575: data is finally put locally.
577: In the sequential case the todata contains indices from where the data is put
578: and the fromdata contains indices from where the data is taken from.
579: This is backwards from the parallel case!
581: .seealso: [](sec_scatter), `VecScatter`, `VecScatterCreate()`
582: @*/
583: PetscErrorCode VecScatterRemap(VecScatter sf, PetscInt tomap[], PetscInt frommap[])
584: {
585: PetscFunctionBegin;
586: if (tomap) PetscAssertPointer(tomap, 2);
587: if (frommap) PetscAssertPointer(frommap, 3);
588: PetscCall(VecScatterRemap_Internal(sf, tomap, frommap));
589: PetscCheck(!frommap, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unable to remap the FROM in scatters yet");
590: /* Mark then vector lengths as unknown because we do not know the lengths of the remapped vectors */
591: sf->vscat.from_n = -1;
592: sf->vscat.to_n = -1;
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
596: /*@
597: VecScatterSetFromOptions - Configures the vector scatter from values in the options database.
599: Collective
601: Input Parameter:
602: . sf - The vector scatter
604: Notes:
605: To see all options, run your program with the -help option, or consult the users manual.
607: Must be called before `VecScatterSetUp()` and before the vector scatter is used.
609: Level: beginner
611: .seealso: [](sec_scatter), `VecScatter`, `VecScatterCreate()`, `VecScatterDestroy()`, `VecScatterSetUp()`
612: @*/
613: PetscErrorCode VecScatterSetFromOptions(VecScatter sf)
614: {
615: PetscFunctionBegin;
617: PetscObjectOptionsBegin((PetscObject)sf);
619: sf->vscat.beginandendtogether = PETSC_FALSE;
620: PetscCall(PetscOptionsBool("-vecscatter_merge", "Use combined (merged) vector scatter begin and end", "VecScatterCreate", sf->vscat.beginandendtogether, &sf->vscat.beginandendtogether, NULL));
621: if (sf->vscat.beginandendtogether) PetscCall(PetscInfo(sf, "Using combined (merged) vector scatter begin and end\n"));
622: PetscOptionsEnd();
623: PetscFunctionReturn(PETSC_SUCCESS);
624: }
626: /*@
627: VecScatterCreate - Creates a vector scatter `VecScatter` context that is used to communicate entries between two vectors `Vec`
629: Collective
631: Input Parameters:
632: + x - a vector that defines the shape (parallel data layout of the vector) of vectors from which we scatter
633: . y - a vector that defines the shape (parallel data layout of the vector) of vectors to which we scatter
634: . ix - the indices of `x` to scatter (if `NULL` scatters all values)
635: - iy - the indices of `y` to hold results (if `NULL` fills entire vector `yin` in order)
637: Output Parameter:
638: . newsf - location to store the new scatter context
640: Options Database Key:
641: . -vecscatter_merge - `VecScatterBegin()` handles all of the communication, `VecScatterEnd()` is a nop
642: eliminates the chance for overlap of computation and communication
644: Level: intermediate
646: Notes:
647: If both `x` and `y` are parallel, their communicator must be on the same
648: set of processes, but their process order can be different.
649: In calls to the scatter options you can use different vectors than the `x` and
650: `y` you used above; BUT they must have the same parallel data layout, for example,
651: they could be obtained from `VecDuplicate()`.
652: A `VecScatter` context CANNOT be used in two or more simultaneous scatters;
653: that is you cannot call a second `VecScatterBegin()` with the same scatter
654: context until the `VecScatterEnd()` has been called on the first `VecScatterBegin()`.
655: In this case a separate `VecScatter` is needed for each concurrent scatter.
657: Both `ix` and `iy` cannot be `NULL` at the same time.
659: Use `VecScatterCreateToAll()` to create a `VecScatter` that copies an MPI vector to sequential vectors on all MPI processes.
660: Use `VecScatterCreateToZero()` to create a `VecScatter` that copies an MPI vector to a sequential vector on MPI rank 0.
661: These special `VecScatter` have better performance than general ones.
663: Developer Note:
664: The implementations of most the `VecScatter` are done using `PetscSF`.
666: .seealso: [](sec_scatter), `VecScatter`, `VecScatterDestroy()`, `VecScatterCreateToAll()`, `VecScatterCreateToZero()`, `PetscSFCreate()`,
667: `VecScatterType`, `InsertMode`, `ScatterMode`, `VecScatterBegin()`, `VecScatterEnd()`, `VecScatterViewFromOptions()`, `VecScatterView()`
668: @*/
669: PetscErrorCode VecScatterCreate(Vec x, IS ix, Vec y, IS iy, VecScatter *newsf)
670: {
671: MPI_Comm xcomm, ycomm, bigcomm;
672: Vec xx, yy;
673: IS ix_old = ix, iy_old = iy, ixx, iyy;
674: PetscMPIInt xcommsize, ycommsize, rank, result;
675: PetscInt i, n, N, nroots, nleaves, *ilocal, xstart, ystart, ixsize, iysize, xlen, ylen;
676: const PetscInt *xindices, *yindices;
677: PetscSFNode *iremote;
678: PetscLayout xlayout, ylayout;
679: ISTypeID ixid, iyid;
680: PetscInt bs, bsx, bsy, min, max, m[2], mg[2], ixfirst, ixstep, iyfirst, iystep;
681: PetscBool can_do_block_opt = PETSC_FALSE;
682: PetscSF sf;
684: PetscFunctionBegin;
685: PetscAssertPointer(newsf, 5);
686: PetscCheck(ix || iy, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot pass default in for both input and output indices");
688: /* Get comm from x and y */
689: PetscCall(PetscObjectGetComm((PetscObject)x, &xcomm));
690: PetscCallMPI(MPI_Comm_size(xcomm, &xcommsize));
691: PetscCall(PetscObjectGetComm((PetscObject)y, &ycomm));
692: PetscCallMPI(MPI_Comm_size(ycomm, &ycommsize));
693: if (xcommsize > 1 && ycommsize > 1) {
694: PetscCallMPI(MPI_Comm_compare(xcomm, ycomm, &result));
695: PetscCheck(result != MPI_UNEQUAL, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMECOMM, "VecScatterCreate: parallel vectors x and y must have identical/congruent/similar communicators");
696: }
697: bs = 1; /* default, no blocking */
699: /*
700: Let P and S stand for parallel and sequential vectors respectively. There are four combinations of vecscatters: PtoP, PtoS,
701: StoP and StoS. The assumption of VecScatterCreate(Vec x,IS ix,Vec y,IS iy,VecScatter *newctx) is: if x is parallel, then ix
702: contains global indices of x. If x is sequential, ix contains local indices of x. Similarly for y and iy.
704: SF builds around concepts of local leaves and remote roots. We treat source vector x as roots and destination vector y as
705: leaves. A PtoS scatter can be naturally mapped to SF. We transform PtoP and StoP to PtoS, and treat StoS as trivial PtoS.
706: */
708: /* NULL ix or iy in VecScatterCreate(x,ix,y,iy,newctx) has special meaning. Recover them for these cases */
709: if (!ix) {
710: if (xcommsize > 1 && ycommsize == 1) { /* PtoS: null ix means the whole x will be scattered to each seq y */
711: PetscCall(VecGetSize(x, &N));
712: PetscCall(ISCreateStride(PETSC_COMM_SELF, N, 0, 1, &ix));
713: } else { /* PtoP, StoP or StoS: null ix means the whole local part of x will be scattered */
714: PetscCall(VecGetLocalSize(x, &n));
715: PetscCall(VecGetOwnershipRange(x, &xstart, NULL));
716: PetscCall(ISCreateStride(PETSC_COMM_SELF, n, xstart, 1, &ix));
717: }
718: }
720: if (!iy) {
721: if (xcommsize == 1 && ycommsize > 1) { /* StoP: null iy means the whole y will be scattered to from each seq x */
722: PetscCall(VecGetSize(y, &N));
723: PetscCall(ISCreateStride(PETSC_COMM_SELF, N, 0, 1, &iy));
724: } else { /* PtoP, StoP or StoS: null iy means the whole local part of y will be scattered to */
725: PetscCall(VecGetLocalSize(y, &n));
726: PetscCall(VecGetOwnershipRange(y, &ystart, NULL));
727: PetscCall(ISCreateStride(PETSC_COMM_SELF, n, ystart, 1, &iy));
728: }
729: }
731: /* Do error checking immediately after we have non-empty ix, iy */
732: PetscCall(ISGetLocalSize(ix, &ixsize));
733: PetscCall(ISGetLocalSize(iy, &iysize));
734: PetscCall(VecGetSize(x, &xlen));
735: PetscCall(VecGetSize(y, &ylen));
736: PetscCheck(ixsize == iysize, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Scatter sizes of ix and iy don't match locally ix=%" PetscInt_FMT " iy=%" PetscInt_FMT, ixsize, iysize);
737: PetscCall(ISGetMinMax(ix, &min, &max));
738: PetscCheck(min >= 0 && max < xlen, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Scatter indices in ix are out of range: found [%" PetscInt_FMT ",%" PetscInt_FMT "), expected in [0,%" PetscInt_FMT ")", min, max, xlen);
739: PetscCall(ISGetMinMax(iy, &min, &max));
740: PetscCheck(min >= 0 && max < ylen, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Scatter indices in iy are out of range: found [%" PetscInt_FMT ",%" PetscInt_FMT "), expected in [0,%" PetscInt_FMT ")", min, max, ylen);
742: /* Extract info about ix, iy for further test */
743: PetscCall(ISGetTypeID_Private(ix, &ixid));
744: PetscCall(ISGetTypeID_Private(iy, &iyid));
745: if (ixid == IS_BLOCK) PetscCall(ISGetBlockSize(ix, &bsx));
746: else if (ixid == IS_STRIDE) PetscCall(ISStrideGetInfo(ix, &ixfirst, &ixstep));
748: if (iyid == IS_BLOCK) PetscCall(ISGetBlockSize(iy, &bsy));
749: else if (iyid == IS_STRIDE) PetscCall(ISStrideGetInfo(iy, &iyfirst, &iystep));
751: /* Check if a PtoS is special ToAll/ToZero scatters, which can be results of VecScatterCreateToAll/Zero.
752: ToAll means a whole MPI vector is copied to a seq vector on every process. ToZero means a whole MPI
753: vector is copied to a seq vector on rank 0 and other processes do nothing(i.e.,they input empty ix,iy).
755: We can optimize these scatters with MPI collectives. We can also avoid costly analysis used for general scatters.
756: */
757: if (xcommsize > 1 && ycommsize == 1) { /* Ranks do not diverge at this if-test */
758: PetscInt pattern[2] = {0, 0}; /* A boolean array with pattern[0] for allgather-like (ToAll) and pattern[1] for gather-like (ToZero) */
759: PetscLayout map;
761: PetscCallMPI(MPI_Comm_rank(xcomm, &rank));
762: PetscCall(VecGetLayout(x, &map));
763: if (rank == 0) {
764: if (ixid == IS_STRIDE && iyid == IS_STRIDE && ixsize == xlen && ixfirst == 0 && ixstep == 1 && iyfirst == 0 && iystep == 1) {
765: /* Rank 0 scatters the whole mpi x to seq y, so it is either a ToAll or a ToZero candidate in its view */
766: pattern[0] = pattern[1] = 1;
767: }
768: } else {
769: if (ixid == IS_STRIDE && iyid == IS_STRIDE && ixsize == xlen && ixfirst == 0 && ixstep == 1 && iyfirst == 0 && iystep == 1) {
770: /* Other ranks also scatter the whole mpi x to seq y, so it is a ToAll candidate in their view */
771: pattern[0] = 1;
772: } else if (ixsize == 0) {
773: /* Other ranks do nothing, so it is a ToZero candidate */
774: pattern[1] = 1;
775: }
776: }
778: /* One stone (the expensive allreduce) two birds: pattern[] tells if it is ToAll or ToZero */
779: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, pattern, 2, MPIU_INT, MPI_LAND, xcomm));
781: if (pattern[0] || pattern[1]) {
782: PetscCall(PetscSFCreate(xcomm, &sf));
783: PetscCall(PetscSFSetFromOptions(sf));
784: PetscCall(PetscSFSetGraphWithPattern(sf, map, pattern[0] ? PETSCSF_PATTERN_ALLGATHER : PETSCSF_PATTERN_GATHER));
785: goto functionend; /* No further analysis needed. What a big win! */
786: }
787: }
789: /* Continue ...
790: Do block optimization by taking advantage of high level info available in ix, iy.
791: The block optimization is valid when all of the following conditions are met:
792: 1) ix, iy are blocked or can be blocked (i.e., strided with step=1);
793: 2) ix, iy have the same block size;
794: 3) all processors agree on one block size;
795: 4) no blocks span more than one process;
796: */
797: bigcomm = (xcommsize == 1) ? ycomm : xcomm;
799: /* Processors could go through different path in this if-else test */
800: m[0] = PETSC_INT_MAX;
801: m[1] = PETSC_INT_MIN;
802: if (ixid == IS_BLOCK && iyid == IS_BLOCK) {
803: m[0] = PetscMin(bsx, bsy);
804: m[1] = PetscMax(bsx, bsy);
805: } else if (ixid == IS_BLOCK && iyid == IS_STRIDE && iystep == 1 && iyfirst % bsx == 0) {
806: m[0] = bsx;
807: m[1] = bsx;
808: } else if (ixid == IS_STRIDE && iyid == IS_BLOCK && ixstep == 1 && ixfirst % bsy == 0) {
809: m[0] = bsy;
810: m[1] = bsy;
811: }
812: /* Get max and min of bsx,bsy over all processes in one allreduce */
813: PetscCall(PetscGlobalMinMaxInt(bigcomm, m, mg));
815: /* Since we used allreduce above, all ranks will have the same min and max. min==max
816: implies all ranks have the same bs. Do further test to see if local vectors are dividable
817: by bs on ALL ranks. If they are, we are ensured that no blocks span more than one processor.
818: */
819: if (mg[0] == mg[1] && mg[0] > 1) {
820: PetscCall(VecGetLocalSize(x, &xlen));
821: PetscCall(VecGetLocalSize(y, &ylen));
822: m[0] = xlen % mg[0];
823: m[1] = ylen % mg[0];
824: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, m, 2, MPIU_INT, MPI_LOR, bigcomm));
825: if (!m[0] && !m[1]) can_do_block_opt = PETSC_TRUE;
826: }
828: /* If can_do_block_opt, then shrink x, y, ix and iy by bs to get xx, yy, ixx and iyy, whose indices
829: and layout are actually used in building SF. Suppose blocked ix representing {0,1,2,6,7,8} has
830: indices {0,2} and bs=3, then ixx = {0,2}; suppose strided iy={3,4,5,6,7,8}, then iyy={1,2}.
832: xx is a little special. If x is seq, then xx is the concatenation of seq x's on ycomm. In this way,
833: we can treat PtoP and StoP uniformly as PtoS.
834: */
835: if (can_do_block_opt) {
836: const PetscInt *indices;
838: /* Shrink x and ix */
839: bs = mg[0];
840: PetscCall(VecCreateMPIWithArray(bigcomm, 1, xlen / bs, PETSC_DECIDE, NULL, &xx)); /* We only care xx's layout */
841: if (ixid == IS_BLOCK) {
842: PetscCall(ISBlockGetIndices(ix, &indices));
843: PetscCall(ISBlockGetLocalSize(ix, &ixsize));
844: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ixsize, indices, PETSC_COPY_VALUES, &ixx));
845: PetscCall(ISBlockRestoreIndices(ix, &indices));
846: } else { /* ixid == IS_STRIDE */
847: PetscCall(ISGetLocalSize(ix, &ixsize));
848: PetscCall(ISCreateStride(PETSC_COMM_SELF, ixsize / bs, ixfirst / bs, 1, &ixx));
849: }
851: /* Shrink y and iy */
852: PetscCall(VecCreateMPIWithArray(ycomm, 1, ylen / bs, PETSC_DECIDE, NULL, &yy));
853: if (iyid == IS_BLOCK) {
854: PetscCall(ISBlockGetIndices(iy, &indices));
855: PetscCall(ISBlockGetLocalSize(iy, &iysize));
856: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, iysize, indices, PETSC_COPY_VALUES, &iyy));
857: PetscCall(ISBlockRestoreIndices(iy, &indices));
858: } else { /* iyid == IS_STRIDE */
859: PetscCall(ISGetLocalSize(iy, &iysize));
860: PetscCall(ISCreateStride(PETSC_COMM_SELF, iysize / bs, iyfirst / bs, 1, &iyy));
861: }
862: } else {
863: ixx = ix;
864: iyy = iy;
865: yy = y;
866: if (xcommsize == 1) PetscCall(VecCreateMPIWithArray(bigcomm, 1, xlen, PETSC_DECIDE, NULL, &xx));
867: else xx = x;
868: }
870: /* Now it is ready to build SF with preprocessed (xx, yy) and (ixx, iyy) */
871: PetscCall(ISGetIndices(ixx, &xindices));
872: PetscCall(ISGetIndices(iyy, &yindices));
873: PetscCall(VecGetLayout(xx, &xlayout));
875: if (ycommsize > 1) {
876: /* PtoP or StoP */
878: /* Below is a piece of complex code with a very simple goal: move global index pairs (xindices[i], yindices[i]),
879: to owner process of yindices[i] according to ylayout, i = 0..n.
881: I did it through a temp sf, but later I thought the old design was inefficient and also distorted log view.
882: We want to map one VecScatterCreate() call to one PetscSFCreate() call. The old design mapped to three
883: PetscSFCreate() calls. This code is on critical path of VecScatterSetUp and is used by every VecScatterCreate.
884: So I commented it out and did another optimized implementation. The commented code is left here for reference.
885: */
886: #if 0
887: const PetscInt *degree;
888: PetscSF tmpsf;
889: PetscInt inedges=0,*leafdata,*rootdata;
891: PetscCall(VecGetOwnershipRange(xx,&xstart,NULL));
892: PetscCall(VecGetLayout(yy,&ylayout));
893: PetscCall(VecGetOwnershipRange(yy,&ystart,NULL));
895: PetscCall(VecGetLocalSize(yy,&nroots));
896: PetscCall(ISGetLocalSize(iyy,&nleaves));
897: PetscCall(PetscMalloc2(nleaves,&iremote,nleaves*2,&leafdata));
899: for (i=0; i<nleaves; i++) {
900: PetscCall(PetscLayoutFindOwnerIndex(ylayout,yindices[i],&iremote[i].rank,&iremote[i].index));
901: leafdata[2*i] = yindices[i];
902: leafdata[2*i+1] = (xcommsize > 1)? xindices[i] : xindices[i] + xstart;
903: }
905: PetscCall(PetscSFCreate(ycomm,&tmpsf));
906: PetscCall(PetscSFSetGraph(tmpsf,nroots,nleaves,NULL,PETSC_USE_POINTER,iremote,PETSC_USE_POINTER));
908: PetscCall(PetscSFComputeDegreeBegin(tmpsf,°ree));
909: PetscCall(PetscSFComputeDegreeEnd(tmpsf,°ree));
911: for (i=0; i<nroots; i++) inedges += degree[i];
912: PetscCall(PetscMalloc1(inedges*2,&rootdata));
913: PetscCall(PetscSFGatherBegin(tmpsf,MPIU_2INT,leafdata,rootdata));
914: PetscCall(PetscSFGatherEnd(tmpsf,MPIU_2INT,leafdata,rootdata));
916: PetscCall(PetscFree2(iremote,leafdata));
917: PetscCall(PetscSFDestroy(&tmpsf));
919: /* rootdata contains global index pairs (i, j). j's are owned by the current process, but i's can point to anywhere.
920: We convert j to local, and convert i to (rank, index). In the end, we get an PtoS suitable for building SF.
921: */
922: nleaves = inedges;
923: PetscCall(VecGetLocalSize(xx,&nroots));
924: PetscCall(PetscMalloc1(nleaves,&ilocal));
925: PetscCall(PetscMalloc1(nleaves,&iremote));
927: for (i=0; i<inedges; i++) {
928: ilocal[i] = rootdata[2*i] - ystart; /* convert y's global index to local index */
929: PetscCall(PetscLayoutFindOwnerIndex(xlayout,rootdata[2*i+1],&iremote[i].rank,&iremote[i].index)); /* convert x's global index to (rank, index) */
930: }
931: PetscCall(PetscFree(rootdata));
932: #else
933: PetscInt j, k, n, disp, rlentotal, *sstart, *xindices_sorted, *yindices_sorted;
934: const PetscInt *yrange;
935: PetscMPIInt nsend, nrecv, nreq, yrank, *sendto, *recvfrom, tag1, tag2;
936: PetscInt *slens, *rlens, count;
937: PetscInt *rxindices, *ryindices;
938: MPI_Request *reqs, *sreqs, *rreqs;
940: /* Sorting makes code simpler, faster and also helps getting rid of many O(P) arrays, which hurt scalability at large scale
941: yindices_sorted - sorted yindices
942: xindices_sorted - xindices sorted along with yindces
943: */
944: PetscCall(ISGetLocalSize(ixx, &n)); /*ixx, iyy have the same local size */
945: PetscCall(PetscMalloc2(n, &xindices_sorted, n, &yindices_sorted));
946: PetscCall(PetscArraycpy(xindices_sorted, xindices, n));
947: PetscCall(PetscArraycpy(yindices_sorted, yindices, n));
948: PetscCall(PetscSortIntWithArray(n, yindices_sorted, xindices_sorted));
949: PetscCall(VecGetOwnershipRange(xx, &xstart, NULL));
950: if (xcommsize == 1) {
951: for (i = 0; i < n; i++) xindices_sorted[i] += xstart;
952: } /* Convert to global indices */
954: /*
955: Calculate info about messages I need to send
956: nsend - number of non-empty messages to send
957: sendto - [nsend] ranks I will send messages to
958: sstart - [nsend+1] sstart[i] is the start index in xsindices_sorted[] I send to rank sendto[i]
959: slens - [ycommsize] I want to send slens[i] entries to rank i.
960: */
961: PetscCall(VecGetLayout(yy, &ylayout));
962: PetscCall(PetscLayoutGetRanges(ylayout, &yrange));
963: PetscCall(PetscCalloc1(ycommsize, &slens)); /* The only O(P) array in this algorithm */
965: i = j = nsend = 0;
966: while (i < n) {
967: if (yindices_sorted[i] >= yrange[j + 1]) { /* If i-th index is out of rank j's bound */
968: do {
969: j++;
970: } while (yindices_sorted[i] >= yrange[j + 1] && j < ycommsize); /* Increase j until i-th index falls in rank j's bound */
971: PetscCheck(j != ycommsize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " not owned by any process, upper bound %" PetscInt_FMT, yindices_sorted[i], yrange[ycommsize]);
972: }
973: i++;
974: if (!slens[j]++) nsend++;
975: }
977: PetscCall(PetscMalloc2(nsend + 1, &sstart, nsend, &sendto));
979: sstart[0] = 0;
980: for (i = j = 0; i < ycommsize; i++) {
981: if (slens[i]) {
982: PetscCall(PetscMPIIntCast(i, &sendto[j]));
983: sstart[j + 1] = sstart[j] + slens[i];
984: j++;
985: }
986: }
988: /*
989: Calculate the reverse info about messages I will recv
990: nrecv - number of messages I will recv
991: recvfrom - [nrecv] ranks I recv from
992: rlens - [nrecv] I will recv rlens[i] entries from rank recvfrom[i]
993: rlentotal - sum of rlens[]
994: rxindices - [rlentotal] recv buffer for xindices_sorted
995: ryindices - [rlentotal] recv buffer for yindices_sorted
996: */
997: PetscCall(PetscGatherNumberOfMessages_Private(ycomm, NULL, slens, &nrecv));
998: PetscCall(PetscGatherMessageLengths_Private(ycomm, nsend, nrecv, slens, &recvfrom, &rlens));
999: PetscCall(PetscFree(slens)); /* Free the O(P) array ASAP */
1000: rlentotal = 0;
1001: for (i = 0; i < nrecv; i++) rlentotal += rlens[i];
1003: /*
1004: Communicate with processors in recvfrom[] to populate rxindices and ryindices
1005: */
1006: PetscCall(PetscCommGetNewTag(ycomm, &tag1));
1007: PetscCall(PetscCommGetNewTag(ycomm, &tag2));
1008: PetscCall(PetscMalloc2(rlentotal, &rxindices, rlentotal, &ryindices));
1009: PetscCall(PetscMPIIntCast((nsend + nrecv) * 2, &nreq));
1010: PetscCall(PetscMalloc1(nreq, &reqs));
1011: sreqs = reqs;
1012: rreqs = PetscSafePointerPlusOffset(reqs, nsend * 2);
1014: for (i = disp = 0; i < nrecv; i++) {
1015: count = rlens[i];
1016: PetscCallMPI(MPIU_Irecv(rxindices + disp, count, MPIU_INT, recvfrom[i], tag1, ycomm, rreqs + i));
1017: PetscCallMPI(MPIU_Irecv(ryindices + disp, count, MPIU_INT, recvfrom[i], tag2, ycomm, rreqs + nrecv + i));
1018: disp += rlens[i];
1019: }
1021: for (i = 0; i < nsend; i++) {
1022: count = sstart[i + 1] - sstart[i];
1023: PetscCallMPI(MPIU_Isend(xindices_sorted + sstart[i], count, MPIU_INT, sendto[i], tag1, ycomm, sreqs + i));
1024: PetscCallMPI(MPIU_Isend(yindices_sorted + sstart[i], count, MPIU_INT, sendto[i], tag2, ycomm, sreqs + nsend + i));
1025: }
1026: PetscCallMPI(MPI_Waitall(nreq, reqs, MPI_STATUSES_IGNORE));
1028: /* Transform VecScatter into SF */
1029: nleaves = rlentotal;
1030: PetscCall(PetscMalloc1(nleaves, &ilocal));
1031: PetscCall(PetscMalloc1(nleaves, &iremote));
1032: PetscCallMPI(MPI_Comm_rank(ycomm, &yrank));
1033: for (i = disp = 0; i < nrecv; i++) {
1034: for (j = 0; j < rlens[i]; j++) {
1035: k = disp + j; /* k-th index pair */
1036: ilocal[k] = ryindices[k] - yrange[yrank]; /* Convert y's global index to local index */
1037: PetscCall(PetscLayoutFindOwnerIndex(xlayout, rxindices[k], &rank, &iremote[k].index)); /* Convert x's global index to (rank, index) */
1038: iremote[k].rank = rank;
1039: }
1040: disp += rlens[i];
1041: }
1043: PetscCall(PetscFree2(sstart, sendto));
1044: PetscCall(PetscFree(rlens));
1045: PetscCall(PetscFree(recvfrom));
1046: PetscCall(PetscFree(reqs));
1047: PetscCall(PetscFree2(rxindices, ryindices));
1048: PetscCall(PetscFree2(xindices_sorted, yindices_sorted));
1049: #endif
1050: } else {
1051: /* PtoS or StoS */
1052: PetscCall(ISGetLocalSize(iyy, &nleaves));
1053: PetscCall(PetscMalloc1(nleaves, &ilocal));
1054: PetscCall(PetscMalloc1(nleaves, &iremote));
1055: PetscCall(PetscArraycpy(ilocal, yindices, nleaves));
1056: for (i = 0; i < nleaves; i++) {
1057: PetscCall(PetscLayoutFindOwnerIndex(xlayout, xindices[i], &rank, &iremote[i].index));
1058: iremote[i].rank = rank;
1059: }
1060: }
1062: /* MUST build SF on xx's comm, which is not necessarily identical to yy's comm.
1063: In SF's view, xx contains the roots (i.e., the remote) and iremote[].rank are ranks in xx's comm.
1064: yy contains leaves, which are local and can be thought as part of PETSC_COMM_SELF. */
1065: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)xx), &sf));
1066: sf->allow_multi_leaves = PETSC_TRUE;
1067: PetscCall(PetscSFSetFromOptions(sf));
1068: PetscCall(VecGetLocalSize(xx, &nroots));
1069: PetscCall(PetscSFSetGraph(sf, nroots, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER)); /* Give ilocal/iremote to PETSc and no need to free them here */
1071: /* Free memory no longer needed */
1072: PetscCall(ISRestoreIndices(ixx, &xindices));
1073: PetscCall(ISRestoreIndices(iyy, &yindices));
1074: if (can_do_block_opt) {
1075: PetscCall(VecDestroy(&xx));
1076: PetscCall(VecDestroy(&yy));
1077: PetscCall(ISDestroy(&ixx));
1078: PetscCall(ISDestroy(&iyy));
1079: } else if (xcommsize == 1) {
1080: PetscCall(VecDestroy(&xx));
1081: }
1083: functionend:
1084: sf->vscat.bs = bs;
1085: if (sf->vscat.bs > 1) {
1086: PetscMPIInt ibs;
1088: PetscCall(PetscMPIIntCast(sf->vscat.bs, &ibs));
1089: PetscCallMPI(MPI_Type_contiguous(ibs, MPIU_SCALAR, &sf->vscat.unit));
1090: PetscCallMPI(MPI_Type_commit(&sf->vscat.unit));
1091: } else {
1092: sf->vscat.unit = MPIU_SCALAR;
1093: }
1094: PetscCall(VecGetLocalSize(x, &sf->vscat.from_n));
1095: PetscCall(VecGetLocalSize(y, &sf->vscat.to_n));
1096: if (!ix_old) PetscCall(ISDestroy(&ix)); /* We created helper ix, iy. Free them */
1097: if (!iy_old) PetscCall(ISDestroy(&iy));
1099: /* Set default */
1100: PetscCall(VecScatterSetFromOptions(sf));
1101: PetscCall(PetscSFSetUp(sf));
1103: *newsf = sf;
1104: PetscFunctionReturn(PETSC_SUCCESS);
1105: }
1107: /*@
1108: VecScatterCreateToAll - Creates a vector and a scatter context that copies all
1109: vector values to each processor
1111: Collective
1113: Input Parameter:
1114: . vin - an `MPIVEC`
1116: Output Parameters:
1117: + ctx - scatter context
1118: - vout - output `SEQVEC` that is large enough to scatter into
1120: Level: intermediate
1122: Example Usage:
1123: .vb
1124: VecScatterCreateToAll(vin, &ctx, &vout);
1126: // scatter as many times as you need
1127: VecScatterBegin(ctx, vin, vout, INSERT_VALUES, SCATTER_FORWARD);
1128: VecScatterEnd(ctx, vin, vout, INSERT_VALUES, SCATTER_FORWARD);
1130: // destroy scatter context and local vector when no longer needed
1131: VecScatterDestroy(&ctx);
1132: VecDestroy(&vout);
1133: .ve
1135: Notes:
1136: `vout` may be `NULL` [`PETSC_NULL_VEC` from Fortran] if you do not
1137: need to have it created
1139: Do NOT create a vector and then pass it in as the final argument `vout`! `vout` is created by this routine
1140: automatically (unless you pass `NULL` in for that argument if you do not need it).
1142: .seealso: [](sec_scatter), `VecScatter`, `VecScatterCreate()`, `VecScatterCreateToZero()`, `VecScatterBegin()`, `VecScatterEnd()`
1143: @*/
1144: PetscErrorCode VecScatterCreateToAll(Vec vin, VecScatter *ctx, Vec *vout)
1145: {
1146: PetscInt N;
1147: IS is;
1148: Vec tmp;
1149: Vec *tmpv;
1150: PetscBool tmpvout = PETSC_FALSE;
1151: VecType roottype;
1153: PetscFunctionBegin;
1156: PetscAssertPointer(ctx, 2);
1157: if (vout) {
1158: PetscAssertPointer(vout, 3);
1159: tmpv = vout;
1160: } else {
1161: tmpvout = PETSC_TRUE;
1162: tmpv = &tmp;
1163: }
1165: /* Create seq vec on each proc, with the same size of the original vec */
1166: PetscCall(VecGetSize(vin, &N));
1167: PetscCall(VecGetRootType_Private(vin, &roottype));
1168: PetscCall(VecCreate(PETSC_COMM_SELF, tmpv));
1169: PetscCall(VecSetSizes(*tmpv, N, PETSC_DECIDE));
1170: PetscCall(VecSetType(*tmpv, roottype));
1171: /* Create the VecScatter ctx with the communication info */
1172: PetscCall(ISCreateStride(PETSC_COMM_SELF, N, 0, 1, &is));
1173: PetscCall(VecScatterCreate(vin, is, *tmpv, is, ctx));
1174: PetscCall(ISDestroy(&is));
1175: if (tmpvout) PetscCall(VecDestroy(tmpv));
1176: PetscFunctionReturn(PETSC_SUCCESS);
1177: }
1179: /*@
1180: VecScatterCreateToZero - Creates an output vector and a scatter context used to
1181: copy all vector values into the output vector on the zeroth processor
1183: Collective
1185: Input Parameter:
1186: . vin - `Vec` of type `MPIVEC`
1188: Output Parameters:
1189: + ctx - scatter context
1190: - vout - output `SEQVEC` that is large enough to scatter into on processor 0 and
1191: of length zero on all other processors
1193: Level: intermediate
1195: Example Usage:
1196: .vb
1197: VecScatterCreateToZero(vin, &ctx, &vout);
1199: // scatter as many times as you need
1200: VecScatterBegin(ctx, vin, vout, INSERT_VALUES, SCATTER_FORWARD);
1201: VecScatterEnd(ctx, vin, vout, INSERT_VALUES, SCATTER_FORWARD);
1203: // destroy scatter context and local vector when no longer needed
1204: VecScatterDestroy(&ctx);
1205: VecDestroy(&vout);
1206: .ve
1208: Notes:
1209: vout may be `NULL` [`PETSC_NULL_VEC` from Fortran] if you do not
1210: need to have it created
1212: Do NOT create a vector and then pass it in as the final argument `vout`! `vout` is created by this routine
1213: automatically (unless you pass `NULL` in for that argument if you do not need it).
1215: .seealso: [](sec_scatter), `VecScatter`, `VecScatterCreate()`, `VecScatterCreateToAll()`, `VecScatterBegin()`, `VecScatterEnd()`
1216: @*/
1217: PetscErrorCode VecScatterCreateToZero(Vec vin, VecScatter *ctx, Vec *vout)
1218: {
1219: PetscInt N;
1220: PetscMPIInt rank;
1221: IS is;
1222: Vec tmp;
1223: Vec *tmpv;
1224: PetscBool tmpvout = PETSC_FALSE;
1225: VecType roottype;
1227: PetscFunctionBegin;
1230: PetscAssertPointer(ctx, 2);
1231: if (vout) {
1232: PetscAssertPointer(vout, 3);
1233: tmpv = vout;
1234: } else {
1235: tmpvout = PETSC_TRUE;
1236: tmpv = &tmp;
1237: }
1239: /* Create vec on each proc, with the same size of the original vec all on process 0 */
1240: PetscCall(VecGetSize(vin, &N));
1241: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)vin), &rank));
1242: if (rank) N = 0;
1243: PetscCall(VecGetRootType_Private(vin, &roottype));
1244: PetscCall(VecCreate(PETSC_COMM_SELF, tmpv));
1245: PetscCall(VecSetSizes(*tmpv, N, PETSC_DECIDE));
1246: PetscCall(VecSetType(*tmpv, roottype));
1247: /* Create the VecScatter ctx with the communication info */
1248: PetscCall(ISCreateStride(PETSC_COMM_SELF, N, 0, 1, &is));
1249: PetscCall(VecScatterCreate(vin, is, *tmpv, is, ctx));
1250: PetscCall(ISDestroy(&is));
1251: if (tmpvout) PetscCall(VecDestroy(tmpv));
1252: PetscFunctionReturn(PETSC_SUCCESS);
1253: }
1255: /*@
1256: VecScatterBegin - Begins a generalized scatter from one vector to
1257: another. Complete the scattering phase with `VecScatterEnd()`.
1259: Neighbor-wise Collective
1261: Input Parameters:
1262: + sf - scatter context generated by `VecScatterCreate()`
1263: . x - the vector from which we scatter
1264: . y - the vector to which we scatter
1265: . addv - either `ADD_VALUES`, `MAX_VALUES`, `MIN_VALUES` or `INSERT_VALUES`, with `INSERT_VALUES` mode any location
1266: not scattered to retains its old value; i.e. the vector is NOT first zeroed.
1267: - mode - the scattering mode, usually `SCATTER_FORWARD`. The available modes are: `SCATTER_FORWARD` or `SCATTER_REVERSE`
1269: Level: intermediate
1271: Notes:
1272: The vectors `x` and `y` need not be the same vectors used in the call
1273: to `VecScatterCreate()`, but `x` must have the same parallel data layout
1274: as that passed in as the `x` to `VecScatterCreate()`, similarly for the `y`.
1275: Most likely they have been obtained from `VecDuplicate()`.
1277: You cannot change the values in the input vector between the calls to `VecScatterBegin()`
1278: and `VecScatterEnd()`.
1280: If you use `SCATTER_REVERSE` the two arguments `x` and `y` should be reversed, from
1281: the `SCATTER_FORWARD`.
1283: .vb
1284: y[iy[i]] = x[ix[i]], for i=0,...,ni-1
1285: .ve
1287: This scatter is far more general than the conventional
1288: scatter, since it can be a gather or a scatter or a combination,
1289: depending on the indices ix and iy. If x is a parallel vector and y
1290: is sequential, `VecScatterBegin()` can serve to gather values to a
1291: single processor. Similarly, if `y` is parallel and `x` sequential, the
1292: routine can scatter from one processor to many processors.
1294: .seealso: [](sec_scatter), `VecScatter`, `VecScatterCreate()`, `VecScatterEnd()`, `InsertMode`, `ScatterMode`
1295: @*/
1296: PetscErrorCode VecScatterBegin(VecScatter sf, Vec x, Vec y, InsertMode addv, ScatterMode mode)
1297: {
1298: PetscInt to_n, from_n;
1300: PetscFunctionBegin;
1304: if (PetscDefined(USE_DEBUG)) {
1305: /*
1306: Error checking to make sure these vectors match the vectors used
1307: to create the vector scatter context. -1 in the from_n and to_n indicate the
1308: vector lengths are unknown (for example with mapped scatters) and thus
1309: no error checking is performed.
1310: */
1311: if (sf->vscat.from_n >= 0 && sf->vscat.to_n >= 0) {
1312: PetscCall(VecGetLocalSize(x, &from_n));
1313: PetscCall(VecGetLocalSize(y, &to_n));
1314: if (mode & SCATTER_REVERSE) {
1315: PetscCheck(to_n == sf->vscat.from_n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Vector wrong size %" PetscInt_FMT " for scatter %" PetscInt_FMT " (scatter reverse and vector to != sf from size)", to_n, sf->vscat.from_n);
1316: PetscCheck(from_n == sf->vscat.to_n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Vector wrong size %" PetscInt_FMT " for scatter %" PetscInt_FMT " (scatter reverse and vector from != sf to size)", from_n, sf->vscat.to_n);
1317: } else {
1318: PetscCheck(to_n == sf->vscat.to_n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Vector wrong size %" PetscInt_FMT " for scatter %" PetscInt_FMT " (scatter forward and vector to != sf to size)", to_n, sf->vscat.to_n);
1319: PetscCheck(from_n == sf->vscat.from_n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Vector wrong size %" PetscInt_FMT " for scatter %" PetscInt_FMT " (scatter forward and vector from != sf from size)", from_n, sf->vscat.from_n);
1320: }
1321: }
1322: }
1324: sf->vscat.logging = PETSC_TRUE;
1325: PetscCall(PetscLogEventBegin(VEC_ScatterBegin, sf, x, y, 0));
1326: PetscCall(VecScatterBegin_Internal(sf, x, y, addv, mode));
1327: if (sf->vscat.beginandendtogether) PetscCall(VecScatterEnd_Internal(sf, x, y, addv, mode));
1328: PetscCall(PetscLogEventEnd(VEC_ScatterBegin, sf, x, y, 0));
1329: sf->vscat.logging = PETSC_FALSE;
1330: PetscFunctionReturn(PETSC_SUCCESS);
1331: }
1333: /*@
1334: VecScatterEnd - Ends a generalized scatter from one vector to another. Call
1335: after first calling `VecScatterBegin()`.
1337: Neighbor-wise Collective
1339: Input Parameters:
1340: + sf - scatter context generated by `VecScatterCreate()`
1341: . x - the vector from which we scatter
1342: . y - the vector to which we scatter
1343: . addv - one of `ADD_VALUES`, `MAX_VALUES`, `MIN_VALUES` or `INSERT_VALUES`
1344: - mode - the scattering mode, usually `SCATTER_FORWARD`. The available modes are: `SCATTER_FORWARD`, `SCATTER_REVERSE`
1346: Level: intermediate
1348: Notes:
1349: If you use `SCATTER_REVERSE` the arguments `x` and `y` should be reversed, from the `SCATTER_FORWARD`.
1351: y[iy[i]] = x[ix[i]], for i=0,...,ni-1
1353: .seealso: [](sec_scatter), `VecScatter`, `VecScatterBegin()`, `VecScatterCreate()`
1354: @*/
1355: PetscErrorCode VecScatterEnd(VecScatter sf, Vec x, Vec y, InsertMode addv, ScatterMode mode)
1356: {
1357: PetscFunctionBegin;
1361: if (!sf->vscat.beginandendtogether) {
1362: sf->vscat.logging = PETSC_TRUE;
1363: PetscCall(PetscLogEventBegin(VEC_ScatterEnd, sf, x, y, 0));
1364: PetscCall(VecScatterEnd_Internal(sf, x, y, addv, mode));
1365: PetscCall(PetscLogEventEnd(VEC_ScatterEnd, sf, x, y, 0));
1366: sf->vscat.logging = PETSC_FALSE;
1367: }
1368: PetscFunctionReturn(PETSC_SUCCESS);
1369: }