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