1 /*
2 This file contains routines for Parallel vector operations.
3 */
4 #include <petscsys.h>
5 #include <../src/vec/vec/impls/mpi/pvecimpl.h> /*I "petscvec.h" I*/
6
7 PETSC_INTERN PetscErrorCode VecView_MPI_Draw(Vec, PetscViewer);
8
VecPlaceArray_MPI(Vec vin,const PetscScalar * a)9 PetscErrorCode VecPlaceArray_MPI(Vec vin, const PetscScalar *a)
10 {
11 Vec_MPI *v = (Vec_MPI *)vin->data;
12
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 }
20
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.
VecDuplicateWithArray_MPI(Vec win,const PetscScalar * array,Vec * v)24 PetscErrorCode VecDuplicateWithArray_MPI(Vec win, const PetscScalar *array, Vec *v)
25 {
26 Vec_MPI *vw, *w = (Vec_MPI *)win->data;
27
28 PetscFunctionBegin;
29 PetscCall(VecCreateWithLayout_Private(win->map, v));
30
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];
34
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;
38
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];
41
42 vw->localupdate = w->localupdate;
43 if (vw->localupdate) PetscCall(PetscObjectReference((PetscObject)vw->localupdate));
44
45 vw->ghost = w->ghost;
46 if (vw->ghost) PetscCall(PetscObjectReference((PetscObject)vw->ghost));
47 }
48
49 /* New vector should inherit stashing property of parent */
50 (*v)->stash.donotstash = win->stash.donotstash;
51 (*v)->stash.ignorenegidx = win->stash.ignorenegidx;
52
53 PetscCall(PetscObjectListDuplicate(((PetscObject)win)->olist, &((PetscObject)*v)->olist));
54 PetscCall(PetscFunctionListDuplicate(((PetscObject)win)->qlist, &((PetscObject)*v)->qlist));
55
56 (*v)->stash.bs = win->stash.bs;
57 (*v)->bstash.bs = win->bstash.bs;
58 PetscFunctionReturn(PETSC_SUCCESS);
59 }
60
VecDuplicate_MPI(Vec win,Vec * v)61 PetscErrorCode VecDuplicate_MPI(Vec win, Vec *v)
62 {
63 PetscFunctionBegin;
64 PetscCall(VecDuplicateWithArray_MPI(win, NULL, v));
65 PetscFunctionReturn(PETSC_SUCCESS);
66 }
67
VecDuplicateVecs_MPI_GEMV(Vec w,PetscInt m,Vec * V[])68 static PetscErrorCode VecDuplicateVecs_MPI_GEMV(Vec w, PetscInt m, Vec *V[])
69 {
70 Vec_MPI *wmpi = (Vec_MPI *)w->data;
71
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"
81
82 PetscCall(PetscMalloc1(m, V));
83 VecGetLocalSizeAligned(w, 64, &lda); // get in lda the 64-bytes aligned local size
84
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 }
107
VecSetOption_MPI(Vec V,VecOption op,PetscBool flag)108 static PetscErrorCode VecSetOption_MPI(Vec V, VecOption op, PetscBool flag)
109 {
110 Vec_MPI *v = (Vec_MPI *)V->data;
111
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 }
130
VecResetArray_MPI(Vec vin)131 PetscErrorCode VecResetArray_MPI(Vec vin)
132 {
133 Vec_MPI *v = (Vec_MPI *)vin->data;
134
135 PetscFunctionBegin;
136 v->array = v->unplacedarray;
137 v->unplacedarray = NULL;
138 if (v->localrep) PetscCall(VecResetArray(v->localrep));
139 PetscFunctionReturn(PETSC_SUCCESS);
140 }
141
VecAssemblySend_MPI_Private(MPI_Comm comm,const PetscMPIInt tag[],PetscMPIInt rankid,PetscMPIInt rank,void * sdata,MPI_Request req[],void * ctx)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;
148
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 }
164
VecAssemblyRecv_MPI_Private(MPI_Comm comm,const PetscMPIInt tag[],PetscMPIInt rank,void * rdata,MPI_Request req[],void * ctx)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;
172
173 PetscFunctionBegin;
174 PetscCall(PetscSegBufferGet(x->segrecvframe, 1, &frame));
175
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 }
187
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 }
201
VecAssemblyBegin_MPI_BTS(Vec X)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;
208
209 PetscFunctionBegin;
210 if (X->stash.donotstash) PetscFunctionReturn(PETSC_SUCCESS);
211
212 PetscCall(PetscObjectGetComm((PetscObject)X, &comm));
213 PetscCall(VecGetBlockSize(X, &bs));
214 if (PetscDefined(USE_DEBUG)) {
215 InsertMode addv;
216
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 */
221
222 PetscCall(VecStashSortCompress_Private(&X->stash));
223 PetscCall(VecStashSortCompress_Private(&X->bstash));
224
225 if (!x->sendranks) {
226 PetscMPIInt nowners, bnowners, *owners, *bowners;
227 PetscInt ntmp;
228
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];
240
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 }
258
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 }
272
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;
277
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 }
287
VecAssemblyEnd_MPI_BTS(Vec X)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;
296
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 }
303
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);
323
324 npending--;
325 if (!blockmsg) { /* Scalar stash */
326 PetscMPIInt count;
327
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;
336
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 }
383
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 }
390
VecAssemblyReset_MPI(Vec X)391 PetscErrorCode VecAssemblyReset_MPI(Vec X)
392 {
393 Vec_MPI *x = (Vec_MPI *)X->data;
394
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 }
408
VecSetFromOptions_MPI(Vec X,PetscOptionItems PetscOptionsObject)409 static PetscErrorCode VecSetFromOptions_MPI(Vec X, PetscOptionItems PetscOptionsObject)
410 {
411 #if !defined(PETSC_HAVE_MPIUNI)
412 PetscBool flg = PETSC_FALSE, set;
413
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 }
429
VecGetLocalToGlobalMapping_MPI_VecGhost(Vec X,ISLocalToGlobalMapping * ltg)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;
435
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 }
454
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 };
547
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()
552
553 If alloc is true and array is NULL then this routine allocates the space, otherwise
554 no space is allocated.
555 */
VecCreate_MPI_Private(Vec v,PetscBool alloc,PetscInt nghost,const PetscScalar array[])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.
561
562 PetscFunctionBegin;
563 PetscCall(PetscNew(&s));
564 v->data = (void *)s;
565 v->ops[0] = DvOps;
566
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));
569
570 // allocate multiple vectors together
571 if (mdot_use_gemv || maxpy_use_gemv) v->ops[0].duplicatevecs = VecDuplicateVecs_MPI_GEMV;
572
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;
580
581 s->nghost = nghost;
582 v->petscnative = PETSC_TRUE;
583 if (array) v->offloadmask = PETSC_OFFLOAD_CPU;
584
585 PetscCall(PetscLayoutSetUp(v->map));
586
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 }
597
598 /* By default parallel vectors do not have local representation */
599 s->localrep = NULL;
600 s->localupdate = NULL;
601 s->ghost = NULL;
602
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));
610
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 }
618
619 /*
620 Create a VECMPI with the given layout and array
621
622 Collective
623
624 Input Parameter:
625 + map - the layout
626 - array - the array on host
627
628 Output Parameter:
629 . V - The vector object
630 */
VecCreateMPIWithLayoutAndArray_Private(PetscLayout map,const PetscScalar array[],Vec * V)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 }
638
639 /*MC
640 VECMPI - VECMPI = "mpi" - The basic parallel vector
641
642 Options Database Key:
643 . -vec_type mpi - sets the vector type to `VECMPI` during a call to `VecSetFromOptions()`
644
645 Level: beginner
646
647 .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecSetType()`, `VecSetFromOptions()`, `VecCreateMPIWithArray()`, `VECMPI`, `VecType`, `VecCreateMPI()`, `VecCreateMPI()`
648 M*/
649
VecCreate_MPI(Vec vv)650 PetscErrorCode VecCreate_MPI(Vec vv)
651 {
652 PetscFunctionBegin;
653 PetscCall(VecCreate_MPI_Private(vv, PETSC_TRUE, 0, NULL));
654 PetscFunctionReturn(PETSC_SUCCESS);
655 }
656
657 /*MC
658 VECSTANDARD = "standard" - A `VECSEQ` on one process and `VECMPI` on more than one process
659
660 Options Database Key:
661 . -vec_type standard - sets a vector type to standard on calls to `VecSetFromOptions()`
662
663 Level: beginner
664
665 .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeq()`, `VecCreateMPI()`
666 M*/
667
VecCreate_Standard(Vec v)668 PETSC_EXTERN PetscErrorCode VecCreate_Standard(Vec v)
669 {
670 PetscMPIInt size;
671
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 }
681
682 /*@
683 VecCreateMPIWithArray - Creates a parallel, array-style vector,
684 where the user provides the array space to store the vector values.
685
686 Collective
687
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
694
695 Output Parameter:
696 . vv - the vector
697
698 Level: intermediate
699
700 Notes:
701 Use `VecDuplicate()` or `VecDuplicateVecs()` to form additional vectors of the
702 same type as an existing vector.
703
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.
706
707 PETSc does NOT free `array` when the vector is destroyed via `VecDestroy()`.
708
709 The user should not free `array` until the vector is destroyed.
710
711 .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeqWithArray()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateGhost()`,
712 `VecCreateMPI()`, `VecCreateGhostWithArray()`, `VecPlaceArray()`
713 @*/
VecCreateMPIWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,const PetscScalar array[],Vec * vv)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 }
725
726 /*@
727 VecCreateGhostWithArray - Creates a parallel vector with ghost padding on each processor;
728 the caller allocates the array space.
729
730 Collective
731
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)
739
740 Output Parameter:
741 . vv - the global vector representation (without ghost points as part of vector)
742
743 Level: advanced
744
745 Notes:
746 Use `VecGhostGetLocalForm()` to access the local, ghosted representation
747 of the vector.
748
749 This also automatically sets the `ISLocalToGlobalMapping()` for this vector.
750
751 .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`,
752 `VecCreateGhost()`, `VecCreateSeqWithArray()`, `VecCreateMPIWithArray()`,
753 `VecCreateGhostBlock()`, `VecCreateGhostBlockWithArray()`, `VecMPISetGhost()`, `VecGhostUpdateBegin()`, `VecGhostUpdateEnd()`
754 @*/
VecCreateGhostWithArray(MPI_Comm comm,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],const PetscScalar array[],Vec * vv)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;
760
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));
775
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));
783
784 w->ghost = from;
785 (*vv)->ops->getlocaltoglobalmapping = VecGetLocalToGlobalMapping_MPI_VecGhost;
786 PetscFunctionReturn(PETSC_SUCCESS);
787 }
788
789 /*@
790 VecGhostGetGhostIS - Return ghosting indices of a ghost vector
791
792 Input Parameters:
793 . X - ghost vector
794
795 Output Parameter:
796 . ghost - ghosting indices
797
798 Level: beginner
799
800 .seealso: `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`
801 @*/
VecGhostGetGhostIS(Vec X,IS * ghost)802 PetscErrorCode VecGhostGetGhostIS(Vec X, IS *ghost)
803 {
804 Vec_MPI *w;
805 PetscBool flg;
806
807 PetscFunctionBegin;
808 PetscValidType(X, 1);
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 }
816
817 /*@
818 VecCreateGhost - Creates a parallel vector with ghost padding on each processor.
819
820 Collective
821
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)
828
829 Output Parameter:
830 . vv - the global vector representation (without ghost points as part of vector)
831
832 Level: advanced
833
834 Notes:
835 Use `VecGhostGetLocalForm()` to access the local, ghosted representation
836 of the vector.
837
838 This also automatically sets the `ISLocalToGlobalMapping()` for this vector.
839
840 .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeq()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateMPI()`,
841 `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`, `VecGhostUpdateBegin()`,
842 `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`, `VecGhostUpdateEnd()`,
843 `VecCreateGhostBlock()`, `VecCreateGhostBlockWithArray()`, `VecMPISetGhost()`
844
845 @*/
VecCreateGhost(MPI_Comm comm,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],Vec * vv)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 }
852
853 /*@
854 VecMPISetGhost - Sets the ghost points for an MPI ghost vector
855
856 Collective
857
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)
862
863 Level: advanced
864
865 Notes:
866 Use `VecGhostGetLocalForm()` to access the local, ghosted representation
867 of the vector.
868
869 This also automatically sets the `ISLocalToGlobalMapping()` for this vector.
870
871 You must call this AFTER you have set the type of the vector (with` VecSetType()`) and the size (with `VecSetSizes()`).
872
873 .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeq()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateMPI()`,
874 `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`, `VecGhostUpdateBegin()`,
875 `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`, `VecGhostUpdateEnd()`,
876 `VecCreateGhostBlock()`, `VecCreateGhostBlockWithArray()`
877 @*/
VecMPISetGhost(Vec vv,PetscInt nghost,const PetscInt ghosts[])878 PetscErrorCode VecMPISetGhost(Vec vv, PetscInt nghost, const PetscInt ghosts[])
879 {
880 PetscBool flg;
881
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;
891
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));
903
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));
911
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 }
920
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.
924
925 Collective
926
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$)
935
936 Output Parameter:
937 . vv - the global vector representation (without ghost points as part of vector)
938
939 Level: advanced
940
941 Notes:
942 Use `VecGhostGetLocalForm()` to access the local, ghosted representation
943 of the vector.
944
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
948
949 .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`,
950 `VecCreateGhost()`, `VecCreateSeqWithArray()`, `VecCreateMPIWithArray()`,
951 `VecCreateGhostWithArray()`, `VecCreateGhostBlock()`, `VecGhostUpdateBegin()`, `VecGhostUpdateEnd()`
952 @*/
VecCreateGhostBlockWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],const PetscScalar array[],Vec * vv)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;
960
961 PetscFunctionBegin;
962 *vv = NULL;
963
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));
979
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));
988
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;
994
995 for (i = 0; i < nb; i++) indices[i] = rstart + i;
996 for (i = 0; i < nghost; i++) indices[nb + i] = ghosts[i];
997
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 }
1003
1004 /*@
1005 VecCreateGhostBlock - Creates a parallel vector with ghost padding on each processor.
1006 The indicing of the ghost points is done with blocks.
1007
1008 Collective
1009
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)
1017
1018 Output Parameter:
1019 . vv - the global vector representation (without ghost points as part of vector)
1020
1021 Level: advanced
1022
1023 Notes:
1024 Use `VecGhostGetLocalForm()` to access the local, ghosted representation
1025 of the vector.
1026
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`
1030
1031 .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeq()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateMPI()`,
1032 `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`, `VecGhostUpdateBegin()`, `VecGhostUpdateEnd()`
1033 `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`, `VecCreateGhostBlockWithArray()`
1034 @*/
VecCreateGhostBlock(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],Vec * vv)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 }
1041