Actual source code: pbvec.c
1: /*
2: This file contains routines for Parallel vector operations.
3: */
4: #include <petscsys.h>
5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
7: PETSC_INTERN PetscErrorCode VecView_MPI_Draw(Vec, PetscViewer);
9: PetscErrorCode VecPlaceArray_MPI(Vec vin, const PetscScalar *a)
10: {
11: Vec_MPI *v = (Vec_MPI *)vin->data;
13: PetscFunctionBegin;
14: PetscCheck(!v->unplacedarray, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "VecPlaceArray() was already called on this vector, without a call to VecResetArray()");
15: v->unplacedarray = v->array; /* save previous array so reset can bring it back */
16: v->array = (PetscScalar *)a;
17: if (v->localrep) PetscCall(VecPlaceArray(v->localrep, a));
18: PetscFunctionReturn(PETSC_SUCCESS);
19: }
21: PetscErrorCode VecDuplicate_MPI(Vec win, Vec *v)
22: {
23: Vec_MPI *vw, *w = (Vec_MPI *)win->data;
24: PetscScalar *array;
26: PetscFunctionBegin;
27: PetscCall(VecCreateWithLayout_Private(win->map, v));
29: PetscCall(VecCreate_MPI_Private(*v, PETSC_TRUE, w->nghost, NULL));
30: vw = (Vec_MPI *)(*v)->data;
31: (*v)->ops[0] = win->ops[0];
33: /* save local representation of the parallel vector (and scatter) if it exists */
34: if (w->localrep) {
35: PetscCall(VecGetArray(*v, &array));
36: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, PetscAbs(win->map->bs), win->map->n + w->nghost, array, &vw->localrep));
37: vw->localrep->ops[0] = w->localrep->ops[0];
38: PetscCall(VecRestoreArray(*v, &array));
40: vw->localupdate = w->localupdate;
41: if (vw->localupdate) PetscCall(PetscObjectReference((PetscObject)vw->localupdate));
43: vw->ghost = w->ghost;
44: if (vw->ghost) PetscCall(PetscObjectReference((PetscObject)vw->ghost));
45: }
47: /* New vector should inherit stashing property of parent */
48: (*v)->stash.donotstash = win->stash.donotstash;
49: (*v)->stash.ignorenegidx = win->stash.ignorenegidx;
51: PetscCall(PetscObjectListDuplicate(((PetscObject)win)->olist, &((PetscObject)*v)->olist));
52: PetscCall(PetscFunctionListDuplicate(((PetscObject)win)->qlist, &((PetscObject)*v)->qlist));
54: (*v)->bstash.bs = win->bstash.bs;
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
58: static PetscErrorCode VecDuplicateVecs_MPI_GEMV(Vec w, PetscInt m, Vec *V[])
59: {
60: Vec_MPI *wmpi = (Vec_MPI *)w->data;
62: PetscFunctionBegin;
63: // Currently only do GEMV for vectors without ghosts. Note w might be a VECMPI subclass object.
64: // This routine relies on the duplicate operation being VecDuplicate_MPI. If not, bail out to the default.
65: if (wmpi->nghost || w->ops->duplicate != VecDuplicate_MPI) {
66: w->ops->duplicatevecs = VecDuplicateVecs_Default;
67: PetscCall(VecDuplicateVecs(w, m, V));
68: } else {
69: PetscInt nlocal;
70: PetscScalar *array;
71: PetscInt64 lda; // use 64-bit as we will do "m * lda"
73: PetscCall(PetscMalloc1(m, V));
74: PetscCall(VecGetLocalSize(w, &nlocal));
75: lda = nlocal;
76: lda = ((lda + 31) / 32) * 32; // make every vector 32-elements aligned
78: PetscCall(PetscCalloc1(m * lda, &array));
79: for (PetscInt i = 0; i < m; i++) {
80: Vec v;
81: PetscCall(VecCreateMPIWithLayoutAndArray_Private(w->map, PetscSafePointerPlusOffset(array, i * lda), &v));
82: PetscCall(PetscObjectListDuplicate(((PetscObject)w)->olist, &((PetscObject)v)->olist));
83: PetscCall(PetscFunctionListDuplicate(((PetscObject)w)->qlist, &((PetscObject)v)->qlist));
84: v->ops->view = w->ops->view;
85: v->stash.donotstash = w->stash.donotstash;
86: v->stash.ignorenegidx = w->stash.ignorenegidx;
87: v->stash.bs = w->stash.bs;
88: (*V)[i] = v;
89: }
90: // So when the first vector is destroyed it will destroy the array
91: if (m) ((Vec_MPI *)(*V)[0]->data)->array_allocated = array;
92: // disable replacearray of the first vector, as freeing its memory also frees others in the group.
93: // But replacearray of others is ok, as they don't own their array.
94: if (m > 1) (*V)[0]->ops->replacearray = VecReplaceArray_Default_GEMV_Error;
95: }
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: static PetscErrorCode VecSetOption_MPI(Vec V, VecOption op, PetscBool flag)
100: {
101: Vec_MPI *v = (Vec_MPI *)V->data;
103: PetscFunctionBegin;
104: switch (op) {
105: case VEC_IGNORE_OFF_PROC_ENTRIES:
106: V->stash.donotstash = flag;
107: break;
108: case VEC_IGNORE_NEGATIVE_INDICES:
109: V->stash.ignorenegidx = flag;
110: break;
111: case VEC_SUBSET_OFF_PROC_ENTRIES:
112: v->assembly_subset = flag; /* See the same logic in MatAssembly wrt MAT_SUBSET_OFF_PROC_ENTRIES */
113: if (!v->assembly_subset) { /* User indicates "do not reuse the communication pattern" */
114: PetscCall(VecAssemblyReset_MPI(V)); /* Reset existing pattern to free memory */
115: v->first_assembly_done = PETSC_FALSE; /* Mark the first assembly is not done */
116: }
117: break;
118: }
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: PetscErrorCode VecResetArray_MPI(Vec vin)
123: {
124: Vec_MPI *v = (Vec_MPI *)vin->data;
126: PetscFunctionBegin;
127: v->array = v->unplacedarray;
128: v->unplacedarray = NULL;
129: if (v->localrep) PetscCall(VecResetArray(v->localrep));
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: static PetscErrorCode VecAssemblySend_MPI_Private(MPI_Comm comm, const PetscMPIInt tag[], PetscMPIInt rankid, PetscMPIInt rank, void *sdata, MPI_Request req[], void *ctx)
134: {
135: Vec X = (Vec)ctx;
136: Vec_MPI *x = (Vec_MPI *)X->data;
137: VecAssemblyHeader *hdr = (VecAssemblyHeader *)sdata;
138: PetscInt bs = X->map->bs;
140: PetscFunctionBegin;
141: /* x->first_assembly_done indicates we are reusing a communication network. In that case, some
142: messages can be empty, but we have to send them this time if we sent them before because the
143: receiver is expecting them.
144: */
145: if (hdr->count || (x->first_assembly_done && x->sendptrs[rankid].ints)) {
146: PetscCallMPI(MPI_Isend(x->sendptrs[rankid].ints, hdr->count, MPIU_INT, rank, tag[0], comm, &req[0]));
147: PetscCallMPI(MPI_Isend(x->sendptrs[rankid].scalars, hdr->count, MPIU_SCALAR, rank, tag[1], comm, &req[1]));
148: }
149: if (hdr->bcount || (x->first_assembly_done && x->sendptrs[rankid].intb)) {
150: PetscCallMPI(MPI_Isend(x->sendptrs[rankid].intb, hdr->bcount, MPIU_INT, rank, tag[2], comm, &req[2]));
151: PetscCallMPI(MPI_Isend(x->sendptrs[rankid].scalarb, hdr->bcount * bs, MPIU_SCALAR, rank, tag[3], comm, &req[3]));
152: }
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: static PetscErrorCode VecAssemblyRecv_MPI_Private(MPI_Comm comm, const PetscMPIInt tag[], PetscMPIInt rank, void *rdata, MPI_Request req[], void *ctx)
157: {
158: Vec X = (Vec)ctx;
159: Vec_MPI *x = (Vec_MPI *)X->data;
160: VecAssemblyHeader *hdr = (VecAssemblyHeader *)rdata;
161: PetscInt bs = X->map->bs;
162: VecAssemblyFrame *frame;
164: PetscFunctionBegin;
165: PetscCall(PetscSegBufferGet(x->segrecvframe, 1, &frame));
167: if (hdr->count) {
168: PetscCall(PetscSegBufferGet(x->segrecvint, hdr->count, &frame->ints));
169: PetscCallMPI(MPI_Irecv(frame->ints, hdr->count, MPIU_INT, rank, tag[0], comm, &req[0]));
170: PetscCall(PetscSegBufferGet(x->segrecvscalar, hdr->count, &frame->scalars));
171: PetscCallMPI(MPI_Irecv(frame->scalars, hdr->count, MPIU_SCALAR, rank, tag[1], comm, &req[1]));
172: frame->pendings = 2;
173: } else {
174: frame->ints = NULL;
175: frame->scalars = NULL;
176: frame->pendings = 0;
177: }
179: if (hdr->bcount) {
180: PetscCall(PetscSegBufferGet(x->segrecvint, hdr->bcount, &frame->intb));
181: PetscCallMPI(MPI_Irecv(frame->intb, hdr->bcount, MPIU_INT, rank, tag[2], comm, &req[2]));
182: PetscCall(PetscSegBufferGet(x->segrecvscalar, hdr->bcount * bs, &frame->scalarb));
183: PetscCallMPI(MPI_Irecv(frame->scalarb, hdr->bcount * bs, MPIU_SCALAR, rank, tag[3], comm, &req[3]));
184: frame->pendingb = 2;
185: } else {
186: frame->intb = NULL;
187: frame->scalarb = NULL;
188: frame->pendingb = 0;
189: }
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: static PetscErrorCode VecAssemblyBegin_MPI_BTS(Vec X)
194: {
195: Vec_MPI *x = (Vec_MPI *)X->data;
196: MPI_Comm comm;
197: PetscInt i, j, jb, bs;
199: PetscFunctionBegin;
200: if (X->stash.donotstash) PetscFunctionReturn(PETSC_SUCCESS);
202: PetscCall(PetscObjectGetComm((PetscObject)X, &comm));
203: PetscCall(VecGetBlockSize(X, &bs));
204: if (PetscDefined(USE_DEBUG)) {
205: InsertMode addv;
206: PetscCall(MPIU_Allreduce((PetscEnum *)&X->stash.insertmode, (PetscEnum *)&addv, 1, MPIU_ENUM, MPI_BOR, comm));
207: PetscCheck(addv != (ADD_VALUES | INSERT_VALUES), comm, PETSC_ERR_ARG_NOTSAMETYPE, "Some processors inserted values while others added");
208: }
209: X->bstash.insertmode = X->stash.insertmode; /* Block stash implicitly tracks InsertMode of scalar stash */
211: PetscCall(VecStashSortCompress_Private(&X->stash));
212: PetscCall(VecStashSortCompress_Private(&X->bstash));
214: if (!x->sendranks) {
215: PetscMPIInt nowners, bnowners, *owners, *bowners;
216: PetscInt ntmp;
217: PetscCall(VecStashGetOwnerList_Private(&X->stash, X->map, &nowners, &owners));
218: PetscCall(VecStashGetOwnerList_Private(&X->bstash, X->map, &bnowners, &bowners));
219: PetscCall(PetscMergeMPIIntArray(nowners, owners, bnowners, bowners, &ntmp, &x->sendranks));
220: x->nsendranks = ntmp;
221: PetscCall(PetscFree(owners));
222: PetscCall(PetscFree(bowners));
223: PetscCall(PetscMalloc1(x->nsendranks, &x->sendhdr));
224: PetscCall(PetscCalloc1(x->nsendranks, &x->sendptrs));
225: }
226: for (i = 0, j = 0, jb = 0; i < x->nsendranks; i++) {
227: PetscMPIInt rank = x->sendranks[i];
228: x->sendhdr[i].insertmode = X->stash.insertmode;
229: /* Initialize pointers for non-empty stashes the first time around. Subsequent assemblies with
230: * VEC_SUBSET_OFF_PROC_ENTRIES will leave the old pointers (dangling because the stash has been collected) when
231: * there is nothing new to send, so that size-zero messages get sent instead. */
232: x->sendhdr[i].count = 0;
233: if (X->stash.n) {
234: x->sendptrs[i].ints = &X->stash.idx[j];
235: x->sendptrs[i].scalars = &X->stash.array[j];
236: for (; j < X->stash.n && X->stash.idx[j] < X->map->range[rank + 1]; j++) x->sendhdr[i].count++;
237: }
238: x->sendhdr[i].bcount = 0;
239: if (X->bstash.n) {
240: x->sendptrs[i].intb = &X->bstash.idx[jb];
241: x->sendptrs[i].scalarb = &X->bstash.array[jb * bs];
242: for (; jb < X->bstash.n && X->bstash.idx[jb] * bs < X->map->range[rank + 1]; jb++) x->sendhdr[i].bcount++;
243: }
244: }
246: if (!x->segrecvint) PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 1000, &x->segrecvint));
247: if (!x->segrecvscalar) PetscCall(PetscSegBufferCreate(sizeof(PetscScalar), 1000, &x->segrecvscalar));
248: if (!x->segrecvframe) PetscCall(PetscSegBufferCreate(sizeof(VecAssemblyFrame), 50, &x->segrecvframe));
249: if (x->first_assembly_done) { /* this is not the first assembly */
250: PetscMPIInt tag[4];
251: for (i = 0; i < 4; i++) PetscCall(PetscCommGetNewTag(comm, &tag[i]));
252: for (i = 0; i < x->nsendranks; i++) PetscCall(VecAssemblySend_MPI_Private(comm, tag, i, x->sendranks[i], x->sendhdr + i, x->sendreqs + 4 * i, X));
253: for (i = 0; i < x->nrecvranks; i++) PetscCall(VecAssemblyRecv_MPI_Private(comm, tag, x->recvranks[i], x->recvhdr + i, x->recvreqs + 4 * i, X));
254: x->use_status = PETSC_TRUE;
255: } else { /* First time assembly */
256: PetscCall(PetscCommBuildTwoSidedFReq(comm, 3, MPIU_INT, x->nsendranks, x->sendranks, (PetscInt *)x->sendhdr, &x->nrecvranks, &x->recvranks, &x->recvhdr, 4, &x->sendreqs, &x->recvreqs, VecAssemblySend_MPI_Private, VecAssemblyRecv_MPI_Private, X));
257: x->use_status = PETSC_FALSE;
258: }
260: /* The first_assembly_done flag is only meaningful when x->assembly_subset is set.
261: This line says when assembly_subset is set, then we mark that the first assembly is done.
262: */
263: x->first_assembly_done = x->assembly_subset;
265: {
266: PetscInt nstash, reallocs;
267: PetscCall(VecStashGetInfo_Private(&X->stash, &nstash, &reallocs));
268: PetscCall(PetscInfo(X, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
269: PetscCall(VecStashGetInfo_Private(&X->bstash, &nstash, &reallocs));
270: PetscCall(PetscInfo(X, "Block-Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
271: }
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: static PetscErrorCode VecAssemblyEnd_MPI_BTS(Vec X)
276: {
277: Vec_MPI *x = (Vec_MPI *)X->data;
278: PetscInt bs = X->map->bs;
279: PetscMPIInt npending, *some_indices, r;
280: MPI_Status *some_statuses;
281: PetscScalar *xarray;
282: VecAssemblyFrame *frame;
284: PetscFunctionBegin;
285: if (X->stash.donotstash) {
286: X->stash.insertmode = NOT_SET_VALUES;
287: X->bstash.insertmode = NOT_SET_VALUES;
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: PetscCheck(x->segrecvframe, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing segrecvframe! Probably you forgot to call VecAssemblyBegin first");
292: PetscCall(VecGetArray(X, &xarray));
293: PetscCall(PetscSegBufferExtractInPlace(x->segrecvframe, &frame));
294: PetscCall(PetscMalloc2(4 * x->nrecvranks, &some_indices, x->use_status ? 4 * x->nrecvranks : 0, &some_statuses));
295: for (r = 0, npending = 0; r < x->nrecvranks; r++) npending += frame[r].pendings + frame[r].pendingb;
296: while (npending > 0) {
297: PetscMPIInt ndone = 0, ii;
298: /* Filling MPI_Status fields requires some resources from the MPI library. We skip it on the first assembly, or
299: * when VEC_SUBSET_OFF_PROC_ENTRIES has not been set, because we could exchange exact sizes in the initial
300: * rendezvous. When the rendezvous is elided, however, we use MPI_Status to get actual message lengths, so that
301: * subsequent assembly can set a proper subset of the values. */
302: PetscCallMPI(MPI_Waitsome(4 * x->nrecvranks, x->recvreqs, &ndone, some_indices, x->use_status ? some_statuses : MPI_STATUSES_IGNORE));
303: for (ii = 0; ii < ndone; ii++) {
304: PetscInt i = some_indices[ii] / 4, j, k;
305: InsertMode imode = (InsertMode)x->recvhdr[i].insertmode;
306: PetscInt *recvint;
307: PetscScalar *recvscalar;
308: PetscBool intmsg = (PetscBool)(some_indices[ii] % 2 == 0);
309: PetscBool blockmsg = (PetscBool)((some_indices[ii] % 4) / 2 == 1);
310: npending--;
311: if (!blockmsg) { /* Scalar stash */
312: PetscMPIInt count;
313: if (--frame[i].pendings > 0) continue;
314: if (x->use_status) {
315: PetscCallMPI(MPI_Get_count(&some_statuses[ii], intmsg ? MPIU_INT : MPIU_SCALAR, &count));
316: } else count = x->recvhdr[i].count;
317: for (j = 0, recvint = frame[i].ints, recvscalar = frame[i].scalars; j < count; j++, recvint++) {
318: PetscInt loc = *recvint - X->map->rstart;
319: PetscCheck(*recvint >= X->map->rstart && X->map->rend > *recvint, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Received vector entry %" PetscInt_FMT " out of local range [%" PetscInt_FMT ",%" PetscInt_FMT ")]", *recvint, X->map->rstart, X->map->rend);
320: switch (imode) {
321: case ADD_VALUES:
322: xarray[loc] += *recvscalar++;
323: break;
324: case INSERT_VALUES:
325: xarray[loc] = *recvscalar++;
326: break;
327: default:
328: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Insert mode not supported 0x%x", imode);
329: }
330: }
331: } else { /* Block stash */
332: PetscMPIInt count;
333: if (--frame[i].pendingb > 0) continue;
334: if (x->use_status) {
335: PetscCallMPI(MPI_Get_count(&some_statuses[ii], intmsg ? MPIU_INT : MPIU_SCALAR, &count));
336: if (!intmsg) count /= bs; /* Convert from number of scalars to number of blocks */
337: } else count = x->recvhdr[i].bcount;
338: for (j = 0, recvint = frame[i].intb, recvscalar = frame[i].scalarb; j < count; j++, recvint++) {
339: PetscInt loc = (*recvint) * bs - X->map->rstart;
340: switch (imode) {
341: case ADD_VALUES:
342: for (k = loc; k < loc + bs; k++) xarray[k] += *recvscalar++;
343: break;
344: case INSERT_VALUES:
345: for (k = loc; k < loc + bs; k++) xarray[k] = *recvscalar++;
346: break;
347: default:
348: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Insert mode not supported 0x%x", imode);
349: }
350: }
351: }
352: }
353: }
354: PetscCall(VecRestoreArray(X, &xarray));
355: PetscCallMPI(MPI_Waitall(4 * x->nsendranks, x->sendreqs, MPI_STATUSES_IGNORE));
356: PetscCall(PetscFree2(some_indices, some_statuses));
357: if (x->assembly_subset) {
358: PetscCall(PetscSegBufferExtractInPlace(x->segrecvint, NULL));
359: PetscCall(PetscSegBufferExtractInPlace(x->segrecvscalar, NULL));
360: } else {
361: PetscCall(VecAssemblyReset_MPI(X));
362: }
364: X->stash.insertmode = NOT_SET_VALUES;
365: X->bstash.insertmode = NOT_SET_VALUES;
366: PetscCall(VecStashScatterEnd_Private(&X->stash));
367: PetscCall(VecStashScatterEnd_Private(&X->bstash));
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: PetscErrorCode VecAssemblyReset_MPI(Vec X)
372: {
373: Vec_MPI *x = (Vec_MPI *)X->data;
375: PetscFunctionBegin;
376: PetscCall(PetscFree(x->sendreqs));
377: PetscCall(PetscFree(x->recvreqs));
378: PetscCall(PetscFree(x->sendranks));
379: PetscCall(PetscFree(x->recvranks));
380: PetscCall(PetscFree(x->sendhdr));
381: PetscCall(PetscFree(x->recvhdr));
382: PetscCall(PetscFree(x->sendptrs));
383: PetscCall(PetscSegBufferDestroy(&x->segrecvint));
384: PetscCall(PetscSegBufferDestroy(&x->segrecvscalar));
385: PetscCall(PetscSegBufferDestroy(&x->segrecvframe));
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: static PetscErrorCode VecSetFromOptions_MPI(Vec X, PetscOptionItems *PetscOptionsObject)
390: {
391: #if !defined(PETSC_HAVE_MPIUNI)
392: PetscBool flg = PETSC_FALSE, set;
394: PetscFunctionBegin;
395: PetscOptionsHeadBegin(PetscOptionsObject, "VecMPI Options");
396: PetscCall(PetscOptionsBool("-vec_assembly_legacy", "Use MPI 1 version of assembly", "", flg, &flg, &set));
397: if (set) {
398: X->ops->assemblybegin = flg ? VecAssemblyBegin_MPI : VecAssemblyBegin_MPI_BTS;
399: X->ops->assemblyend = flg ? VecAssemblyEnd_MPI : VecAssemblyEnd_MPI_BTS;
400: }
401: PetscOptionsHeadEnd();
402: #else
403: PetscFunctionBegin;
404: X->ops->assemblybegin = VecAssemblyBegin_MPI;
405: X->ops->assemblyend = VecAssemblyEnd_MPI;
406: #endif
407: PetscFunctionReturn(PETSC_SUCCESS);
408: }
410: static PetscErrorCode VecGetLocalToGlobalMapping_MPI_VecGhost(Vec X, ISLocalToGlobalMapping *ltg)
411: {
412: PetscInt *indices, n, nghost, rstart, i;
413: IS ghostis;
414: const PetscInt *ghostidx;
415: MPI_Comm comm;
417: PetscFunctionBegin;
418: if (X->map->mapping) {
419: *ltg = X->map->mapping;
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
422: PetscCall(VecGhostGetGhostIS(X, &ghostis));
423: PetscCall(ISGetLocalSize(ghostis, &nghost));
424: PetscCall(VecGetLocalSize(X, &n));
425: PetscCall(ISGetIndices(ghostis, &ghostidx));
426: /* set local to global mapping for ghosted vector */
427: PetscCall(PetscMalloc1(n + nghost, &indices));
428: PetscCall(VecGetOwnershipRange(X, &rstart, NULL));
429: for (i = 0; i < n; i++) { indices[i] = rstart + i; }
430: for (i = 0; i < nghost; i++) { indices[n + i] = ghostidx[i]; }
431: PetscCall(ISRestoreIndices(ghostis, &ghostidx));
432: PetscCall(PetscObjectGetComm((PetscObject)X, &comm));
433: PetscCall(ISLocalToGlobalMappingCreate(comm, 1, n + nghost, indices, PETSC_OWN_POINTER, &X->map->mapping));
434: *ltg = X->map->mapping;
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: static struct _VecOps DvOps = {PetscDesignatedInitializer(duplicate, VecDuplicate_MPI), /* 1 */
439: PetscDesignatedInitializer(duplicatevecs, VecDuplicateVecs_Default),
440: PetscDesignatedInitializer(destroyvecs, VecDestroyVecs_Default),
441: PetscDesignatedInitializer(dot, VecDot_MPI),
442: PetscDesignatedInitializer(mdot, VecMDot_MPI),
443: PetscDesignatedInitializer(norm, VecNorm_MPI),
444: PetscDesignatedInitializer(tdot, VecTDot_MPI),
445: PetscDesignatedInitializer(mtdot, VecMTDot_MPI),
446: PetscDesignatedInitializer(scale, VecScale_Seq),
447: PetscDesignatedInitializer(copy, VecCopy_Seq), /* 10 */
448: PetscDesignatedInitializer(set, VecSet_Seq),
449: PetscDesignatedInitializer(swap, VecSwap_Seq),
450: PetscDesignatedInitializer(axpy, VecAXPY_Seq),
451: PetscDesignatedInitializer(axpby, VecAXPBY_Seq),
452: PetscDesignatedInitializer(maxpy, VecMAXPY_Seq),
453: PetscDesignatedInitializer(aypx, VecAYPX_Seq),
454: PetscDesignatedInitializer(waxpy, VecWAXPY_Seq),
455: PetscDesignatedInitializer(axpbypcz, VecAXPBYPCZ_Seq),
456: PetscDesignatedInitializer(pointwisemult, VecPointwiseMult_Seq),
457: PetscDesignatedInitializer(pointwisedivide, VecPointwiseDivide_Seq),
458: PetscDesignatedInitializer(setvalues, VecSetValues_MPI), /* 20 */
459: PetscDesignatedInitializer(assemblybegin, VecAssemblyBegin_MPI_BTS),
460: PetscDesignatedInitializer(assemblyend, VecAssemblyEnd_MPI_BTS),
461: PetscDesignatedInitializer(getarray, NULL),
462: PetscDesignatedInitializer(getsize, VecGetSize_MPI),
463: PetscDesignatedInitializer(getlocalsize, VecGetSize_Seq),
464: PetscDesignatedInitializer(restorearray, NULL),
465: PetscDesignatedInitializer(max, VecMax_MPI),
466: PetscDesignatedInitializer(min, VecMin_MPI),
467: PetscDesignatedInitializer(setrandom, VecSetRandom_Seq),
468: PetscDesignatedInitializer(setoption, VecSetOption_MPI),
469: PetscDesignatedInitializer(setvaluesblocked, VecSetValuesBlocked_MPI),
470: PetscDesignatedInitializer(destroy, VecDestroy_MPI),
471: PetscDesignatedInitializer(view, VecView_MPI),
472: PetscDesignatedInitializer(placearray, VecPlaceArray_MPI),
473: PetscDesignatedInitializer(replacearray, VecReplaceArray_Seq),
474: PetscDesignatedInitializer(dot_local, VecDot_Seq),
475: PetscDesignatedInitializer(tdot_local, VecTDot_Seq),
476: PetscDesignatedInitializer(norm_local, VecNorm_Seq),
477: PetscDesignatedInitializer(mdot_local, VecMDot_Seq),
478: PetscDesignatedInitializer(mtdot_local, VecMTDot_Seq),
479: PetscDesignatedInitializer(load, VecLoad_Default),
480: PetscDesignatedInitializer(reciprocal, VecReciprocal_Default),
481: PetscDesignatedInitializer(conjugate, VecConjugate_Seq),
482: PetscDesignatedInitializer(setlocaltoglobalmapping, NULL),
483: PetscDesignatedInitializer(getlocaltoglobalmapping, VecGetLocalToGlobalMapping_MPI_VecGhost),
484: PetscDesignatedInitializer(setvalueslocal, NULL),
485: PetscDesignatedInitializer(resetarray, VecResetArray_MPI),
486: PetscDesignatedInitializer(setfromoptions, VecSetFromOptions_MPI), /*set from options */
487: PetscDesignatedInitializer(maxpointwisedivide, VecMaxPointwiseDivide_Seq),
488: PetscDesignatedInitializer(pointwisemax, VecPointwiseMax_Seq),
489: PetscDesignatedInitializer(pointwisemaxabs, VecPointwiseMaxAbs_Seq),
490: PetscDesignatedInitializer(pointwisemin, VecPointwiseMin_Seq),
491: PetscDesignatedInitializer(getvalues, VecGetValues_MPI),
492: PetscDesignatedInitializer(sqrt, NULL),
493: PetscDesignatedInitializer(abs, NULL),
494: PetscDesignatedInitializer(exp, NULL),
495: PetscDesignatedInitializer(log, NULL),
496: PetscDesignatedInitializer(shift, NULL),
497: PetscDesignatedInitializer(create, NULL), /* really? */
498: PetscDesignatedInitializer(stridegather, VecStrideGather_Default),
499: PetscDesignatedInitializer(stridescatter, VecStrideScatter_Default),
500: PetscDesignatedInitializer(dotnorm2, NULL),
501: PetscDesignatedInitializer(getsubvector, NULL),
502: PetscDesignatedInitializer(restoresubvector, NULL),
503: PetscDesignatedInitializer(getarrayread, NULL),
504: PetscDesignatedInitializer(restorearrayread, NULL),
505: PetscDesignatedInitializer(stridesubsetgather, VecStrideSubSetGather_Default),
506: PetscDesignatedInitializer(stridesubsetscatter, VecStrideSubSetScatter_Default),
507: PetscDesignatedInitializer(viewnative, VecView_MPI),
508: PetscDesignatedInitializer(loadnative, NULL),
509: PetscDesignatedInitializer(createlocalvector, NULL),
510: PetscDesignatedInitializer(getlocalvector, NULL),
511: PetscDesignatedInitializer(restorelocalvector, NULL),
512: PetscDesignatedInitializer(getlocalvectorread, NULL),
513: PetscDesignatedInitializer(restorelocalvectorread, NULL),
514: PetscDesignatedInitializer(bindtocpu, NULL),
515: PetscDesignatedInitializer(getarraywrite, NULL),
516: PetscDesignatedInitializer(restorearraywrite, NULL),
517: PetscDesignatedInitializer(getarrayandmemtype, NULL),
518: PetscDesignatedInitializer(restorearrayandmemtype, NULL),
519: PetscDesignatedInitializer(getarrayreadandmemtype, NULL),
520: PetscDesignatedInitializer(restorearrayreadandmemtype, NULL),
521: PetscDesignatedInitializer(getarraywriteandmemtype, NULL),
522: PetscDesignatedInitializer(restorearraywriteandmemtype, NULL),
523: PetscDesignatedInitializer(concatenate, NULL),
524: PetscDesignatedInitializer(sum, NULL),
525: PetscDesignatedInitializer(setpreallocationcoo, VecSetPreallocationCOO_MPI),
526: PetscDesignatedInitializer(setvaluescoo, VecSetValuesCOO_MPI),
527: PetscDesignatedInitializer(errorwnorm, NULL)};
529: /*
530: VecCreate_MPI_Private - Basic create routine called by VecCreate_MPI() (i.e. VecCreateMPI()),
531: VecCreateMPIWithArray(), VecCreate_Shared() (i.e. VecCreateShared()), VecCreateGhost(),
532: VecDuplicate_MPI(), VecCreateGhostWithArray(), VecDuplicate_MPI(), and VecDuplicate_Shared()
534: If alloc is true and array is NULL then this routine allocates the space, otherwise
535: no space is allocated.
536: */
537: PetscErrorCode VecCreate_MPI_Private(Vec v, PetscBool alloc, PetscInt nghost, const PetscScalar array[])
538: {
539: Vec_MPI *s;
540: PetscBool mdot_use_gemv = PETSC_TRUE;
541: PetscBool maxpy_use_gemv = PETSC_FALSE; // default is false as we saw bad performance with vendors' GEMV with tall skinny matrices.
543: PetscFunctionBegin;
544: PetscCall(PetscNew(&s));
545: v->data = (void *)s;
546: v->ops[0] = DvOps;
548: PetscCall(PetscOptionsGetBool(NULL, NULL, "-vec_mdot_use_gemv", &mdot_use_gemv, NULL));
549: PetscCall(PetscOptionsGetBool(NULL, NULL, "-vec_maxpy_use_gemv", &maxpy_use_gemv, NULL));
551: // allocate multiple vectors together
552: if (mdot_use_gemv || maxpy_use_gemv) v->ops[0].duplicatevecs = VecDuplicateVecs_MPI_GEMV;
554: if (mdot_use_gemv) {
555: v->ops[0].duplicatevecs = VecDuplicateVecs_MPI_GEMV;
556: v->ops[0].mdot = VecMDot_MPI_GEMV;
557: v->ops[0].mdot_local = VecMDot_Seq_GEMV;
558: v->ops[0].mtdot = VecMTDot_MPI_GEMV;
559: v->ops[0].mtdot_local = VecMTDot_Seq_GEMV;
560: }
561: if (maxpy_use_gemv) v->ops[0].maxpy = VecMAXPY_Seq_GEMV;
563: s->nghost = nghost;
564: v->petscnative = PETSC_TRUE;
565: if (array) v->offloadmask = PETSC_OFFLOAD_CPU;
567: PetscCall(PetscLayoutSetUp(v->map));
569: s->array = (PetscScalar *)array;
570: s->array_allocated = NULL;
571: if (alloc && !array) {
572: PetscInt n = v->map->n + nghost;
573: PetscCall(PetscCalloc1(n, &s->array));
574: s->array_allocated = s->array;
575: PetscCall(PetscObjectComposedDataSetReal((PetscObject)v, NormIds[NORM_2], 0));
576: PetscCall(PetscObjectComposedDataSetReal((PetscObject)v, NormIds[NORM_1], 0));
577: PetscCall(PetscObjectComposedDataSetReal((PetscObject)v, NormIds[NORM_INFINITY], 0));
578: }
580: /* By default parallel vectors do not have local representation */
581: s->localrep = NULL;
582: s->localupdate = NULL;
583: s->ghost = NULL;
585: v->stash.insertmode = NOT_SET_VALUES;
586: v->bstash.insertmode = NOT_SET_VALUES;
587: /* create the stashes. The block-size for bstash is set later when
588: VecSetValuesBlocked is called.
589: */
590: PetscCall(VecStashCreate_Private(PetscObjectComm((PetscObject)v), 1, &v->stash));
591: PetscCall(VecStashCreate_Private(PetscObjectComm((PetscObject)v), PetscAbs(v->map->bs), &v->bstash));
593: #if defined(PETSC_HAVE_MATLAB)
594: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", VecMatlabEnginePut_Default));
595: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", VecMatlabEngineGet_Default));
596: #endif
597: PetscCall(PetscObjectChangeTypeName((PetscObject)v, VECMPI));
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }
601: /*
602: Create a VECMPI with the given layout and array
604: Collective
606: Input Parameter:
607: + map - the layout
608: - array - the array on host
610: Output Parameter:
611: . V - The vector object
612: */
613: PetscErrorCode VecCreateMPIWithLayoutAndArray_Private(PetscLayout map, const PetscScalar array[], Vec *V)
614: {
615: PetscFunctionBegin;
616: PetscCall(VecCreateWithLayout_Private(map, V));
617: PetscCall(VecCreate_MPI_Private(*V, PETSC_FALSE, 0, array));
618: PetscFunctionReturn(PETSC_SUCCESS);
619: }
621: /*MC
622: VECMPI - VECMPI = "mpi" - The basic parallel vector
624: Options Database Key:
625: . -vec_type mpi - sets the vector type to `VECMPI` during a call to `VecSetFromOptions()`
627: Level: beginner
629: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecSetType()`, `VecSetFromOptions()`, `VecCreateMPIWithArray()`, `VECMPI`, `VecType`, `VecCreateMPI()`, `VecCreateMPI()`
630: M*/
632: PetscErrorCode VecCreate_MPI(Vec vv)
633: {
634: PetscFunctionBegin;
635: PetscCall(VecCreate_MPI_Private(vv, PETSC_TRUE, 0, NULL));
636: PetscFunctionReturn(PETSC_SUCCESS);
637: }
639: /*MC
640: VECSTANDARD = "standard" - A `VECSEQ` on one process and `VECMPI` on more than one process
642: Options Database Key:
643: . -vec_type standard - sets a vector type to standard on calls to `VecSetFromOptions()`
645: Level: beginner
647: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeq()`, `VecCreateMPI()`
648: M*/
650: PETSC_EXTERN PetscErrorCode VecCreate_Standard(Vec v)
651: {
652: PetscMPIInt size;
654: PetscFunctionBegin;
655: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
656: if (size == 1) {
657: PetscCall(VecSetType(v, VECSEQ));
658: } else {
659: PetscCall(VecSetType(v, VECMPI));
660: }
661: PetscFunctionReturn(PETSC_SUCCESS);
662: }
664: /*@C
665: VecCreateMPIWithArray - Creates a parallel, array-style vector,
666: where the user provides the array space to store the vector values.
668: Collective
670: Input Parameters:
671: + comm - the MPI communicator to use
672: . bs - block size, same meaning as `VecSetBlockSize()`
673: . n - local vector length, cannot be `PETSC_DECIDE`
674: . N - global vector length (or `PETSC_DETERMINE` to have calculated)
675: - array - the user provided array to store the vector values
677: Output Parameter:
678: . vv - the vector
680: Level: intermediate
682: Notes:
683: Use `VecDuplicate()` or `VecDuplicateVecs()` to form additional vectors of the
684: same type as an existing vector.
686: If the user-provided array is `NULL`, then `VecPlaceArray()` can be used
687: at a later stage to SET the array for storing the vector values.
689: PETSc does NOT free `array` when the vector is destroyed via `VecDestroy()`.
691: The user should not free `array` until the vector is destroyed.
693: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeqWithArray()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateGhost()`,
694: `VecCreateMPI()`, `VecCreateGhostWithArray()`, `VecPlaceArray()`
695: @*/
696: PetscErrorCode VecCreateMPIWithArray(MPI_Comm comm, PetscInt bs, PetscInt n, PetscInt N, const PetscScalar array[], Vec *vv)
697: {
698: PetscFunctionBegin;
699: PetscCheck(n != PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must set local size of vector");
700: PetscCall(PetscSplitOwnership(comm, &n, &N));
701: PetscCall(VecCreate(comm, vv));
702: PetscCall(VecSetSizes(*vv, n, N));
703: PetscCall(VecSetBlockSize(*vv, bs));
704: PetscCall(VecCreate_MPI_Private(*vv, PETSC_FALSE, 0, array));
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: /*@C
709: VecCreateGhostWithArray - Creates a parallel vector with ghost padding on each processor;
710: the caller allocates the array space.
712: Collective
714: Input Parameters:
715: + comm - the MPI communicator to use
716: . n - local vector length
717: . N - global vector length (or `PETSC_DETERMINE` to have calculated if `n` is given)
718: . nghost - number of local ghost points
719: . ghosts - global indices of ghost points (or `NULL` if not needed), these do not need to be in increasing order (sorted)
720: - array - the space to store the vector values (as long as n + nghost)
722: Output Parameter:
723: . vv - the global vector representation (without ghost points as part of vector)
725: Level: advanced
727: Notes:
728: Use `VecGhostGetLocalForm()` to access the local, ghosted representation
729: of the vector.
731: This also automatically sets the `ISLocalToGlobalMapping()` for this vector.
733: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`,
734: `VecCreateGhost()`, `VecCreateSeqWithArray()`, `VecCreateMPIWithArray()`,
735: `VecCreateGhostBlock()`, `VecCreateGhostBlockWithArray()`, `VecMPISetGhost()`
736: @*/
737: PetscErrorCode VecCreateGhostWithArray(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[], const PetscScalar array[], Vec *vv)
738: {
739: Vec_MPI *w;
740: PetscScalar *larray;
741: IS from, to;
743: PetscFunctionBegin;
744: *vv = NULL;
746: PetscCheck(n != PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must set local size");
747: PetscCheck(nghost != PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must set local ghost size");
748: PetscCheck(nghost >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Ghost length must be >= 0");
749: PetscCall(PetscSplitOwnership(comm, &n, &N));
750: /* Create global representation */
751: PetscCall(VecCreate(comm, vv));
752: PetscCall(VecSetSizes(*vv, n, N));
753: PetscCall(VecCreate_MPI_Private(*vv, PETSC_TRUE, nghost, array));
754: w = (Vec_MPI *)(*vv)->data;
755: /* Create local representation */
756: PetscCall(VecGetArray(*vv, &larray));
757: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, n + nghost, larray, &w->localrep));
758: PetscCall(VecRestoreArray(*vv, &larray));
760: /*
761: Create scatter context for scattering (updating) ghost values
762: */
763: PetscCall(ISCreateGeneral(comm, nghost, ghosts, PETSC_COPY_VALUES, &from));
764: PetscCall(ISCreateStride(PETSC_COMM_SELF, nghost, n, 1, &to));
765: PetscCall(VecScatterCreate(*vv, from, w->localrep, to, &w->localupdate));
766: PetscCall(ISDestroy(&to));
768: w->ghost = from;
769: (*vv)->ops->getlocaltoglobalmapping = VecGetLocalToGlobalMapping_MPI_VecGhost;
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
773: /*@
774: VecGhostGetGhostIS - Return ghosting indices of a ghost vector
776: Input Parameters:
777: . X - ghost vector
779: Output Parameter:
780: . ghost - ghosting indices
782: Level: beginner
784: .seealso: `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`
785: @*/
786: PetscErrorCode VecGhostGetGhostIS(Vec X, IS *ghost)
787: {
788: Vec_MPI *w;
789: PetscBool flg;
791: PetscFunctionBegin;
793: PetscAssertPointer(ghost, 2);
794: PetscCall(PetscObjectTypeCompare((PetscObject)X, VECMPI, &flg));
795: PetscCheck(flg, PetscObjectComm((PetscObject)X), PETSC_ERR_ARG_WRONGSTATE, "VecGhostGetGhostIS was not supported for vec type %s", ((PetscObject)X)->type_name);
796: w = (Vec_MPI *)(X)->data;
797: *ghost = w->ghost;
798: PetscFunctionReturn(PETSC_SUCCESS);
799: }
801: /*@
802: VecCreateGhost - Creates a parallel vector with ghost padding on each processor.
804: Collective
806: Input Parameters:
807: + comm - the MPI communicator to use
808: . n - local vector length
809: . N - global vector length (or `PETSC_DETERMINE` to have calculated if `n` is given)
810: . nghost - number of local ghost points
811: - ghosts - global indices of ghost points, these do not need to be in increasing order (sorted)
813: Output Parameter:
814: . vv - the global vector representation (without ghost points as part of vector)
816: Level: advanced
818: Notes:
819: Use `VecGhostGetLocalForm()` to access the local, ghosted representation
820: of the vector.
822: This also automatically sets the `ISLocalToGlobalMapping()` for this vector.
824: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeq()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateMPI()`,
825: `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`, `VecGhostUpdateBegin()`,
826: `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`, `VecGhostUpdateEnd()`,
827: `VecCreateGhostBlock()`, `VecCreateGhostBlockWithArray()`, `VecMPISetGhost()`
829: @*/
830: PetscErrorCode VecCreateGhost(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[], Vec *vv)
831: {
832: PetscFunctionBegin;
833: PetscCall(VecCreateGhostWithArray(comm, n, N, nghost, ghosts, NULL, vv));
834: PetscFunctionReturn(PETSC_SUCCESS);
835: }
837: /*@
838: VecMPISetGhost - Sets the ghost points for an MPI ghost vector
840: Collective
842: Input Parameters:
843: + vv - the MPI vector
844: . nghost - number of local ghost points
845: - ghosts - global indices of ghost points, these do not need to be in increasing order (sorted)
847: Level: advanced
849: Notes:
850: Use `VecGhostGetLocalForm()` to access the local, ghosted representation
851: of the vector.
853: This also automatically sets the `ISLocalToGlobalMapping()` for this vector.
855: You must call this AFTER you have set the type of the vector (with` VecSetType()`) and the size (with `VecSetSizes()`).
857: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeq()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateMPI()`,
858: `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`, `VecGhostUpdateBegin()`,
859: `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`, `VecGhostUpdateEnd()`,
860: `VecCreateGhostBlock()`, `VecCreateGhostBlockWithArray()`
861: @*/
862: PetscErrorCode VecMPISetGhost(Vec vv, PetscInt nghost, const PetscInt ghosts[])
863: {
864: PetscBool flg;
866: PetscFunctionBegin;
867: PetscCall(PetscObjectTypeCompare((PetscObject)vv, VECMPI, &flg));
868: /* if already fully existent VECMPI then basically destroy it and rebuild with ghosting */
869: if (flg) {
870: PetscInt n, N;
871: Vec_MPI *w;
872: PetscScalar *larray;
873: IS from, to;
874: MPI_Comm comm;
876: PetscCall(PetscObjectGetComm((PetscObject)vv, &comm));
877: n = vv->map->n;
878: N = vv->map->N;
879: PetscUseTypeMethod(vv, destroy);
880: PetscCall(VecSetSizes(vv, n, N));
881: PetscCall(VecCreate_MPI_Private(vv, PETSC_TRUE, nghost, NULL));
882: w = (Vec_MPI *)(vv)->data;
883: /* Create local representation */
884: PetscCall(VecGetArray(vv, &larray));
885: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, n + nghost, larray, &w->localrep));
886: PetscCall(VecRestoreArray(vv, &larray));
888: /*
889: Create scatter context for scattering (updating) ghost values
890: */
891: PetscCall(ISCreateGeneral(comm, nghost, ghosts, PETSC_COPY_VALUES, &from));
892: PetscCall(ISCreateStride(PETSC_COMM_SELF, nghost, n, 1, &to));
893: PetscCall(VecScatterCreate(vv, from, w->localrep, to, &w->localupdate));
894: PetscCall(ISDestroy(&to));
896: w->ghost = from;
897: vv->ops->getlocaltoglobalmapping = VecGetLocalToGlobalMapping_MPI_VecGhost;
898: } else {
899: PetscCheck(vv->ops->create != VecCreate_MPI, PetscObjectComm((PetscObject)vv), PETSC_ERR_ARG_WRONGSTATE, "Must set local or global size before setting ghosting");
900: PetscCheck(((PetscObject)vv)->type_name, PetscObjectComm((PetscObject)vv), PETSC_ERR_ARG_WRONGSTATE, "Must set type to VECMPI before ghosting");
901: }
902: PetscFunctionReturn(PETSC_SUCCESS);
903: }
905: /* ------------------------------------------------------------------------------------------*/
906: /*@C
907: VecCreateGhostBlockWithArray - Creates a parallel vector with ghost padding on each processor;
908: the caller allocates the array space. Indices in the ghost region are based on blocks.
910: Collective
912: Input Parameters:
913: + comm - the MPI communicator to use
914: . bs - block size
915: . n - local vector length
916: . N - global vector length (or `PETSC_DETERMINE` to have calculated if `n` is given)
917: . nghost - number of local ghost blocks
918: . ghosts - global indices of ghost blocks (or `NULL` if not needed), counts are by block not by index, these do not need to be in increasing order (sorted)
919: - array - the space to store the vector values (as long as `n + nghost*bs`)
921: Output Parameter:
922: . vv - the global vector representation (without ghost points as part of vector)
924: Level: advanced
926: Notes:
927: Use `VecGhostGetLocalForm()` to access the local, ghosted representation
928: of the vector.
930: n is the local vector size (total local size not the number of blocks) while nghost
931: is the number of blocks in the ghost portion, i.e. the number of elements in the ghost
932: portion is bs*nghost
934: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`,
935: `VecCreateGhost()`, `VecCreateSeqWithArray()`, `VecCreateMPIWithArray()`,
936: `VecCreateGhostWithArray()`, `VecCreateGhostBlock()`
937: @*/
938: PetscErrorCode VecCreateGhostBlockWithArray(MPI_Comm comm, PetscInt bs, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[], const PetscScalar array[], Vec *vv)
939: {
940: Vec_MPI *w;
941: PetscScalar *larray;
942: IS from, to;
943: ISLocalToGlobalMapping ltog;
944: PetscInt rstart, i, nb, *indices;
946: PetscFunctionBegin;
947: *vv = NULL;
949: PetscCheck(n != PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must set local size");
950: PetscCheck(nghost != PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must set local ghost size");
951: PetscCheck(nghost >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Ghost length must be >= 0");
952: PetscCheck(n % bs == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local size must be a multiple of block size");
953: PetscCall(PetscSplitOwnership(comm, &n, &N));
954: /* Create global representation */
955: PetscCall(VecCreate(comm, vv));
956: PetscCall(VecSetSizes(*vv, n, N));
957: PetscCall(VecSetBlockSize(*vv, bs));
958: PetscCall(VecCreate_MPI_Private(*vv, PETSC_TRUE, nghost * bs, array));
959: w = (Vec_MPI *)(*vv)->data;
960: /* Create local representation */
961: PetscCall(VecGetArray(*vv, &larray));
962: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, bs, n + bs * nghost, larray, &w->localrep));
963: PetscCall(VecRestoreArray(*vv, &larray));
965: /*
966: Create scatter context for scattering (updating) ghost values
967: */
968: PetscCall(ISCreateBlock(comm, bs, nghost, ghosts, PETSC_COPY_VALUES, &from));
969: PetscCall(ISCreateStride(PETSC_COMM_SELF, bs * nghost, n, 1, &to));
970: PetscCall(VecScatterCreate(*vv, from, w->localrep, to, &w->localupdate));
971: PetscCall(ISDestroy(&to));
972: PetscCall(ISDestroy(&from));
974: /* set local to global mapping for ghosted vector */
975: nb = n / bs;
976: PetscCall(PetscMalloc1(nb + nghost, &indices));
977: PetscCall(VecGetOwnershipRange(*vv, &rstart, NULL));
978: rstart = rstart / bs;
980: for (i = 0; i < nb; i++) indices[i] = rstart + i;
981: for (i = 0; i < nghost; i++) indices[nb + i] = ghosts[i];
983: PetscCall(ISLocalToGlobalMappingCreate(comm, bs, nb + nghost, indices, PETSC_OWN_POINTER, <og));
984: PetscCall(VecSetLocalToGlobalMapping(*vv, ltog));
985: PetscCall(ISLocalToGlobalMappingDestroy(<og));
986: PetscFunctionReturn(PETSC_SUCCESS);
987: }
989: /*@
990: VecCreateGhostBlock - Creates a parallel vector with ghost padding on each processor.
991: The indicing of the ghost points is done with blocks.
993: Collective
995: Input Parameters:
996: + comm - the MPI communicator to use
997: . bs - the block size
998: . n - local vector length
999: . N - global vector length (or `PETSC_DETERMINE` to have calculated if `n` is given)
1000: . nghost - number of local ghost blocks
1001: - ghosts - global indices of ghost blocks, counts are by block, not by individual index, these do not need to be in increasing order (sorted)
1003: Output Parameter:
1004: . vv - the global vector representation (without ghost points as part of vector)
1006: Level: advanced
1008: Notes:
1009: Use `VecGhostGetLocalForm()` to access the local, ghosted representation
1010: of the vector.
1012: `n` is the local vector size (total local size not the number of blocks) while `nghost`
1013: is the number of blocks in the ghost portion, i.e. the number of elements in the ghost
1014: portion is `bs*nghost`
1016: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeq()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateMPI()`,
1017: `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`,
1018: `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`, `VecCreateGhostBlockWithArray()`
1019: @*/
1020: PetscErrorCode VecCreateGhostBlock(MPI_Comm comm, PetscInt bs, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[], Vec *vv)
1021: {
1022: PetscFunctionBegin;
1023: PetscCall(VecCreateGhostBlockWithArray(comm, bs, n, N, nghost, ghosts, NULL, vv));
1024: PetscFunctionReturn(PETSC_SUCCESS);
1025: }