Actual source code: pbvec.c
2: /*
3: This file contains routines for Parallel vector operations.
4: */
5: #include <petscsys.h>
6: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
8: PetscErrorCode VecDot_MPI(Vec xin, Vec yin, PetscScalar *z)
9: {
10: PetscScalar sum, work;
12: VecDot_Seq(xin, yin, &work);
13: MPIU_Allreduce(&work, &sum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)xin));
14: *z = sum;
15: return 0;
16: }
18: PetscErrorCode VecTDot_MPI(Vec xin, Vec yin, PetscScalar *z)
19: {
20: PetscScalar sum, work;
22: VecTDot_Seq(xin, yin, &work);
23: MPIU_Allreduce(&work, &sum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)xin));
24: *z = sum;
25: return 0;
26: }
28: extern PetscErrorCode VecView_MPI_Draw(Vec, PetscViewer);
30: static PetscErrorCode VecPlaceArray_MPI(Vec vin, const PetscScalar *a)
31: {
32: Vec_MPI *v = (Vec_MPI *)vin->data;
35: v->unplacedarray = v->array; /* save previous array so reset can bring it back */
36: v->array = (PetscScalar *)a;
37: if (v->localrep) VecPlaceArray(v->localrep, a);
38: return 0;
39: }
41: PetscErrorCode VecDuplicate_MPI(Vec win, Vec *v)
42: {
43: Vec_MPI *vw, *w = (Vec_MPI *)win->data;
44: PetscScalar *array;
46: VecCreate(PetscObjectComm((PetscObject)win), v);
47: PetscLayoutReference(win->map, &(*v)->map);
49: VecCreate_MPI_Private(*v, PETSC_TRUE, w->nghost, NULL);
50: vw = (Vec_MPI *)(*v)->data;
51: PetscMemcpy((*v)->ops, win->ops, sizeof(struct _VecOps));
53: /* save local representation of the parallel vector (and scatter) if it exists */
54: if (w->localrep) {
55: VecGetArray(*v, &array);
56: VecCreateSeqWithArray(PETSC_COMM_SELF, PetscAbs(win->map->bs), win->map->n + w->nghost, array, &vw->localrep);
57: PetscMemcpy(vw->localrep->ops, w->localrep->ops, sizeof(struct _VecOps));
58: VecRestoreArray(*v, &array);
60: vw->localupdate = w->localupdate;
61: if (vw->localupdate) PetscObjectReference((PetscObject)vw->localupdate);
62: }
64: /* New vector should inherit stashing property of parent */
65: (*v)->stash.donotstash = win->stash.donotstash;
66: (*v)->stash.ignorenegidx = win->stash.ignorenegidx;
68: PetscObjectListDuplicate(((PetscObject)win)->olist, &((PetscObject)(*v))->olist);
69: PetscFunctionListDuplicate(((PetscObject)win)->qlist, &((PetscObject)(*v))->qlist);
71: (*v)->map->bs = PetscAbs(win->map->bs);
72: (*v)->bstash.bs = win->bstash.bs;
73: return 0;
74: }
76: static PetscErrorCode VecSetOption_MPI(Vec V, VecOption op, PetscBool flag)
77: {
78: Vec_MPI *v = (Vec_MPI *)V->data;
80: switch (op) {
81: case VEC_IGNORE_OFF_PROC_ENTRIES:
82: V->stash.donotstash = flag;
83: break;
84: case VEC_IGNORE_NEGATIVE_INDICES:
85: V->stash.ignorenegidx = flag;
86: break;
87: case VEC_SUBSET_OFF_PROC_ENTRIES:
88: v->assembly_subset = flag; /* See the same logic in MatAssembly wrt MAT_SUBSET_OFF_PROC_ENTRIES */
89: if (!v->assembly_subset) { /* User indicates "do not reuse the communication pattern" */
90: VecAssemblyReset_MPI(V); /* Reset existing pattern to free memory */
91: v->first_assembly_done = PETSC_FALSE; /* Mark the first assembly is not done */
92: }
93: break;
94: }
96: return 0;
97: }
99: static PetscErrorCode VecResetArray_MPI(Vec vin)
100: {
101: Vec_MPI *v = (Vec_MPI *)vin->data;
103: v->array = v->unplacedarray;
104: v->unplacedarray = NULL;
105: if (v->localrep) VecResetArray(v->localrep);
106: return 0;
107: }
109: static PetscErrorCode VecAssemblySend_MPI_Private(MPI_Comm comm, const PetscMPIInt tag[], PetscMPIInt rankid, PetscMPIInt rank, void *sdata, MPI_Request req[], void *ctx)
110: {
111: Vec X = (Vec)ctx;
112: Vec_MPI *x = (Vec_MPI *)X->data;
113: VecAssemblyHeader *hdr = (VecAssemblyHeader *)sdata;
114: PetscInt bs = X->map->bs;
116: /* x->first_assembly_done indicates we are reusing a communication network. In that case, some
117: messages can be empty, but we have to send them this time if we sent them before because the
118: receiver is expecting them.
119: */
120: if (hdr->count || (x->first_assembly_done && x->sendptrs[rankid].ints)) {
121: MPI_Isend(x->sendptrs[rankid].ints, hdr->count, MPIU_INT, rank, tag[0], comm, &req[0]);
122: MPI_Isend(x->sendptrs[rankid].scalars, hdr->count, MPIU_SCALAR, rank, tag[1], comm, &req[1]);
123: }
124: if (hdr->bcount || (x->first_assembly_done && x->sendptrs[rankid].intb)) {
125: MPI_Isend(x->sendptrs[rankid].intb, hdr->bcount, MPIU_INT, rank, tag[2], comm, &req[2]);
126: MPI_Isend(x->sendptrs[rankid].scalarb, hdr->bcount * bs, MPIU_SCALAR, rank, tag[3], comm, &req[3]);
127: }
128: return 0;
129: }
131: static PetscErrorCode VecAssemblyRecv_MPI_Private(MPI_Comm comm, const PetscMPIInt tag[], PetscMPIInt rank, void *rdata, MPI_Request req[], void *ctx)
132: {
133: Vec X = (Vec)ctx;
134: Vec_MPI *x = (Vec_MPI *)X->data;
135: VecAssemblyHeader *hdr = (VecAssemblyHeader *)rdata;
136: PetscInt bs = X->map->bs;
137: VecAssemblyFrame *frame;
139: PetscSegBufferGet(x->segrecvframe, 1, &frame);
141: if (hdr->count) {
142: PetscSegBufferGet(x->segrecvint, hdr->count, &frame->ints);
143: MPI_Irecv(frame->ints, hdr->count, MPIU_INT, rank, tag[0], comm, &req[0]);
144: PetscSegBufferGet(x->segrecvscalar, hdr->count, &frame->scalars);
145: MPI_Irecv(frame->scalars, hdr->count, MPIU_SCALAR, rank, tag[1], comm, &req[1]);
146: frame->pendings = 2;
147: } else {
148: frame->ints = NULL;
149: frame->scalars = NULL;
150: frame->pendings = 0;
151: }
153: if (hdr->bcount) {
154: PetscSegBufferGet(x->segrecvint, hdr->bcount, &frame->intb);
155: MPI_Irecv(frame->intb, hdr->bcount, MPIU_INT, rank, tag[2], comm, &req[2]);
156: PetscSegBufferGet(x->segrecvscalar, hdr->bcount * bs, &frame->scalarb);
157: MPI_Irecv(frame->scalarb, hdr->bcount * bs, MPIU_SCALAR, rank, tag[3], comm, &req[3]);
158: frame->pendingb = 2;
159: } else {
160: frame->intb = NULL;
161: frame->scalarb = NULL;
162: frame->pendingb = 0;
163: }
164: return 0;
165: }
167: static PetscErrorCode VecAssemblyBegin_MPI_BTS(Vec X)
168: {
169: Vec_MPI *x = (Vec_MPI *)X->data;
170: MPI_Comm comm;
171: PetscInt i, j, jb, bs;
173: if (X->stash.donotstash) return 0;
175: PetscObjectGetComm((PetscObject)X, &comm);
176: VecGetBlockSize(X, &bs);
177: if (PetscDefined(USE_DEBUG)) {
178: InsertMode addv;
179: MPIU_Allreduce((PetscEnum *)&X->stash.insertmode, (PetscEnum *)&addv, 1, MPIU_ENUM, MPI_BOR, comm);
181: }
182: X->bstash.insertmode = X->stash.insertmode; /* Block stash implicitly tracks InsertMode of scalar stash */
184: VecStashSortCompress_Private(&X->stash);
185: VecStashSortCompress_Private(&X->bstash);
187: if (!x->sendranks) {
188: PetscMPIInt nowners, bnowners, *owners, *bowners;
189: PetscInt ntmp;
190: VecStashGetOwnerList_Private(&X->stash, X->map, &nowners, &owners);
191: VecStashGetOwnerList_Private(&X->bstash, X->map, &bnowners, &bowners);
192: PetscMergeMPIIntArray(nowners, owners, bnowners, bowners, &ntmp, &x->sendranks);
193: x->nsendranks = ntmp;
194: PetscFree(owners);
195: PetscFree(bowners);
196: PetscMalloc1(x->nsendranks, &x->sendhdr);
197: PetscCalloc1(x->nsendranks, &x->sendptrs);
198: }
199: for (i = 0, j = 0, jb = 0; i < x->nsendranks; i++) {
200: PetscMPIInt rank = x->sendranks[i];
201: x->sendhdr[i].insertmode = X->stash.insertmode;
202: /* Initialize pointers for non-empty stashes the first time around. Subsequent assemblies with
203: * VEC_SUBSET_OFF_PROC_ENTRIES will leave the old pointers (dangling because the stash has been collected) when
204: * there is nothing new to send, so that size-zero messages get sent instead. */
205: x->sendhdr[i].count = 0;
206: if (X->stash.n) {
207: x->sendptrs[i].ints = &X->stash.idx[j];
208: x->sendptrs[i].scalars = &X->stash.array[j];
209: for (; j < X->stash.n && X->stash.idx[j] < X->map->range[rank + 1]; j++) x->sendhdr[i].count++;
210: }
211: x->sendhdr[i].bcount = 0;
212: if (X->bstash.n) {
213: x->sendptrs[i].intb = &X->bstash.idx[jb];
214: x->sendptrs[i].scalarb = &X->bstash.array[jb * bs];
215: for (; jb < X->bstash.n && X->bstash.idx[jb] * bs < X->map->range[rank + 1]; jb++) x->sendhdr[i].bcount++;
216: }
217: }
219: if (!x->segrecvint) PetscSegBufferCreate(sizeof(PetscInt), 1000, &x->segrecvint);
220: if (!x->segrecvscalar) PetscSegBufferCreate(sizeof(PetscScalar), 1000, &x->segrecvscalar);
221: if (!x->segrecvframe) PetscSegBufferCreate(sizeof(VecAssemblyFrame), 50, &x->segrecvframe);
222: if (x->first_assembly_done) { /* this is not the first assembly */
223: PetscMPIInt tag[4];
224: for (i = 0; i < 4; i++) PetscCommGetNewTag(comm, &tag[i]);
225: for (i = 0; i < x->nsendranks; i++) VecAssemblySend_MPI_Private(comm, tag, i, x->sendranks[i], x->sendhdr + i, x->sendreqs + 4 * i, X);
226: for (i = 0; i < x->nrecvranks; i++) VecAssemblyRecv_MPI_Private(comm, tag, x->recvranks[i], x->recvhdr + i, x->recvreqs + 4 * i, X);
227: x->use_status = PETSC_TRUE;
228: } else { /* First time assembly */
229: 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);
230: x->use_status = PETSC_FALSE;
231: }
233: /* The first_assembly_done flag is only meaningful when x->assembly_subset is set.
234: This line says when assembly_subset is set, then we mark that the first assembly is done.
235: */
236: x->first_assembly_done = x->assembly_subset;
238: {
239: PetscInt nstash, reallocs;
240: VecStashGetInfo_Private(&X->stash, &nstash, &reallocs);
241: PetscInfo(X, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs);
242: VecStashGetInfo_Private(&X->bstash, &nstash, &reallocs);
243: PetscInfo(X, "Block-Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs);
244: }
245: return 0;
246: }
248: static PetscErrorCode VecAssemblyEnd_MPI_BTS(Vec X)
249: {
250: Vec_MPI *x = (Vec_MPI *)X->data;
251: PetscInt bs = X->map->bs;
252: PetscMPIInt npending, *some_indices, r;
253: MPI_Status *some_statuses;
254: PetscScalar *xarray;
255: VecAssemblyFrame *frame;
257: if (X->stash.donotstash) {
258: X->stash.insertmode = NOT_SET_VALUES;
259: X->bstash.insertmode = NOT_SET_VALUES;
260: return 0;
261: }
264: VecGetArray(X, &xarray);
265: PetscSegBufferExtractInPlace(x->segrecvframe, &frame);
266: PetscMalloc2(4 * x->nrecvranks, &some_indices, x->use_status ? 4 * x->nrecvranks : 0, &some_statuses);
267: for (r = 0, npending = 0; r < x->nrecvranks; r++) npending += frame[r].pendings + frame[r].pendingb;
268: while (npending > 0) {
269: PetscMPIInt ndone = 0, ii;
270: /* Filling MPI_Status fields requires some resources from the MPI library. We skip it on the first assembly, or
271: * when VEC_SUBSET_OFF_PROC_ENTRIES has not been set, because we could exchange exact sizes in the initial
272: * rendezvous. When the rendezvous is elided, however, we use MPI_Status to get actual message lengths, so that
273: * subsequent assembly can set a proper subset of the values. */
274: MPI_Waitsome(4 * x->nrecvranks, x->recvreqs, &ndone, some_indices, x->use_status ? some_statuses : MPI_STATUSES_IGNORE);
275: for (ii = 0; ii < ndone; ii++) {
276: PetscInt i = some_indices[ii] / 4, j, k;
277: InsertMode imode = (InsertMode)x->recvhdr[i].insertmode;
278: PetscInt *recvint;
279: PetscScalar *recvscalar;
280: PetscBool intmsg = (PetscBool)(some_indices[ii] % 2 == 0);
281: PetscBool blockmsg = (PetscBool)((some_indices[ii] % 4) / 2 == 1);
282: npending--;
283: if (!blockmsg) { /* Scalar stash */
284: PetscMPIInt count;
285: if (--frame[i].pendings > 0) continue;
286: if (x->use_status) {
287: MPI_Get_count(&some_statuses[ii], intmsg ? MPIU_INT : MPIU_SCALAR, &count);
288: } else count = x->recvhdr[i].count;
289: for (j = 0, recvint = frame[i].ints, recvscalar = frame[i].scalars; j < count; j++, recvint++) {
290: PetscInt loc = *recvint - X->map->rstart;
292: switch (imode) {
293: case ADD_VALUES:
294: xarray[loc] += *recvscalar++;
295: break;
296: case INSERT_VALUES:
297: xarray[loc] = *recvscalar++;
298: break;
299: default:
300: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Insert mode not supported 0x%x", imode);
301: }
302: }
303: } else { /* Block stash */
304: PetscMPIInt count;
305: if (--frame[i].pendingb > 0) continue;
306: if (x->use_status) {
307: MPI_Get_count(&some_statuses[ii], intmsg ? MPIU_INT : MPIU_SCALAR, &count);
308: if (!intmsg) count /= bs; /* Convert from number of scalars to number of blocks */
309: } else count = x->recvhdr[i].bcount;
310: for (j = 0, recvint = frame[i].intb, recvscalar = frame[i].scalarb; j < count; j++, recvint++) {
311: PetscInt loc = (*recvint) * bs - X->map->rstart;
312: switch (imode) {
313: case ADD_VALUES:
314: for (k = loc; k < loc + bs; k++) xarray[k] += *recvscalar++;
315: break;
316: case INSERT_VALUES:
317: for (k = loc; k < loc + bs; k++) xarray[k] = *recvscalar++;
318: break;
319: default:
320: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Insert mode not supported 0x%x", imode);
321: }
322: }
323: }
324: }
325: }
326: VecRestoreArray(X, &xarray);
327: MPI_Waitall(4 * x->nsendranks, x->sendreqs, MPI_STATUSES_IGNORE);
328: PetscFree2(some_indices, some_statuses);
329: if (x->assembly_subset) {
330: PetscSegBufferExtractInPlace(x->segrecvint, NULL);
331: PetscSegBufferExtractInPlace(x->segrecvscalar, NULL);
332: } else {
333: VecAssemblyReset_MPI(X);
334: }
336: X->stash.insertmode = NOT_SET_VALUES;
337: X->bstash.insertmode = NOT_SET_VALUES;
338: VecStashScatterEnd_Private(&X->stash);
339: VecStashScatterEnd_Private(&X->bstash);
340: return 0;
341: }
343: PetscErrorCode VecAssemblyReset_MPI(Vec X)
344: {
345: Vec_MPI *x = (Vec_MPI *)X->data;
347: PetscFree(x->sendreqs);
348: PetscFree(x->recvreqs);
349: PetscFree(x->sendranks);
350: PetscFree(x->recvranks);
351: PetscFree(x->sendhdr);
352: PetscFree(x->recvhdr);
353: PetscFree(x->sendptrs);
354: PetscSegBufferDestroy(&x->segrecvint);
355: PetscSegBufferDestroy(&x->segrecvscalar);
356: PetscSegBufferDestroy(&x->segrecvframe);
357: return 0;
358: }
360: static PetscErrorCode VecSetFromOptions_MPI(Vec X, PetscOptionItems *PetscOptionsObject)
361: {
362: #if !defined(PETSC_HAVE_MPIUNI)
363: PetscBool flg = PETSC_FALSE, set;
365: PetscOptionsHeadBegin(PetscOptionsObject, "VecMPI Options");
366: PetscOptionsBool("-vec_assembly_legacy", "Use MPI 1 version of assembly", "", flg, &flg, &set);
367: if (set) {
368: X->ops->assemblybegin = flg ? VecAssemblyBegin_MPI : VecAssemblyBegin_MPI_BTS;
369: X->ops->assemblyend = flg ? VecAssemblyEnd_MPI : VecAssemblyEnd_MPI_BTS;
370: }
371: PetscOptionsHeadEnd();
372: #else
373: X->ops->assemblybegin = VecAssemblyBegin_MPI;
374: X->ops->assemblyend = VecAssemblyEnd_MPI;
375: #endif
376: return 0;
377: }
379: static struct _VecOps DvOps = {PetscDesignatedInitializer(duplicate, VecDuplicate_MPI), /* 1 */
380: PetscDesignatedInitializer(duplicatevecs, VecDuplicateVecs_Default),
381: PetscDesignatedInitializer(destroyvecs, VecDestroyVecs_Default),
382: PetscDesignatedInitializer(dot, VecDot_MPI),
383: PetscDesignatedInitializer(mdot, VecMDot_MPI),
384: PetscDesignatedInitializer(norm, VecNorm_MPI),
385: PetscDesignatedInitializer(tdot, VecTDot_MPI),
386: PetscDesignatedInitializer(mtdot, VecMTDot_MPI),
387: PetscDesignatedInitializer(scale, VecScale_Seq),
388: PetscDesignatedInitializer(copy, VecCopy_Seq), /* 10 */
389: PetscDesignatedInitializer(set, VecSet_Seq),
390: PetscDesignatedInitializer(swap, VecSwap_Seq),
391: PetscDesignatedInitializer(axpy, VecAXPY_Seq),
392: PetscDesignatedInitializer(axpby, VecAXPBY_Seq),
393: PetscDesignatedInitializer(maxpy, VecMAXPY_Seq),
394: PetscDesignatedInitializer(aypx, VecAYPX_Seq),
395: PetscDesignatedInitializer(waxpy, VecWAXPY_Seq),
396: PetscDesignatedInitializer(axpbypcz, VecAXPBYPCZ_Seq),
397: PetscDesignatedInitializer(pointwisemult, VecPointwiseMult_Seq),
398: PetscDesignatedInitializer(pointwisedivide, VecPointwiseDivide_Seq),
399: PetscDesignatedInitializer(setvalues, VecSetValues_MPI), /* 20 */
400: PetscDesignatedInitializer(assemblybegin, VecAssemblyBegin_MPI_BTS),
401: PetscDesignatedInitializer(assemblyend, VecAssemblyEnd_MPI_BTS),
402: PetscDesignatedInitializer(getarray, NULL),
403: PetscDesignatedInitializer(getsize, VecGetSize_MPI),
404: PetscDesignatedInitializer(getlocalsize, VecGetSize_Seq),
405: PetscDesignatedInitializer(restorearray, NULL),
406: PetscDesignatedInitializer(max, VecMax_MPI),
407: PetscDesignatedInitializer(min, VecMin_MPI),
408: PetscDesignatedInitializer(setrandom, VecSetRandom_Seq),
409: PetscDesignatedInitializer(setoption, VecSetOption_MPI),
410: PetscDesignatedInitializer(setvaluesblocked, VecSetValuesBlocked_MPI),
411: PetscDesignatedInitializer(destroy, VecDestroy_MPI),
412: PetscDesignatedInitializer(view, VecView_MPI),
413: PetscDesignatedInitializer(placearray, VecPlaceArray_MPI),
414: PetscDesignatedInitializer(replacearray, VecReplaceArray_Seq),
415: PetscDesignatedInitializer(dot_local, VecDot_Seq),
416: PetscDesignatedInitializer(tdot_local, VecTDot_Seq),
417: PetscDesignatedInitializer(norm_local, VecNorm_Seq),
418: PetscDesignatedInitializer(mdot_local, VecMDot_Seq),
419: PetscDesignatedInitializer(mtdot_local, VecMTDot_Seq),
420: PetscDesignatedInitializer(load, VecLoad_Default),
421: PetscDesignatedInitializer(reciprocal, VecReciprocal_Default),
422: PetscDesignatedInitializer(conjugate, VecConjugate_Seq),
423: PetscDesignatedInitializer(setlocaltoglobalmapping, NULL),
424: PetscDesignatedInitializer(setvalueslocal, NULL),
425: PetscDesignatedInitializer(resetarray, VecResetArray_MPI),
426: PetscDesignatedInitializer(setfromoptions, VecSetFromOptions_MPI), /*set from options */
427: PetscDesignatedInitializer(maxpointwisedivide, VecMaxPointwiseDivide_Seq),
428: PetscDesignatedInitializer(pointwisemax, VecPointwiseMax_Seq),
429: PetscDesignatedInitializer(pointwisemaxabs, VecPointwiseMaxAbs_Seq),
430: PetscDesignatedInitializer(pointwisemin, VecPointwiseMin_Seq),
431: PetscDesignatedInitializer(getvalues, VecGetValues_MPI),
432: PetscDesignatedInitializer(sqrt, NULL),
433: PetscDesignatedInitializer(abs, NULL),
434: PetscDesignatedInitializer(exp, NULL),
435: PetscDesignatedInitializer(log, NULL),
436: PetscDesignatedInitializer(shift, NULL),
437: PetscDesignatedInitializer(create, NULL), /* really? */
438: PetscDesignatedInitializer(stridegather, VecStrideGather_Default),
439: PetscDesignatedInitializer(stridescatter, VecStrideScatter_Default),
440: PetscDesignatedInitializer(dotnorm2, NULL),
441: PetscDesignatedInitializer(getsubvector, NULL),
442: PetscDesignatedInitializer(restoresubvector, NULL),
443: PetscDesignatedInitializer(getarrayread, NULL),
444: PetscDesignatedInitializer(restorearrayread, NULL),
445: PetscDesignatedInitializer(stridesubsetgather, VecStrideSubSetGather_Default),
446: PetscDesignatedInitializer(stridesubsetscatter, VecStrideSubSetScatter_Default),
447: PetscDesignatedInitializer(viewnative, VecView_MPI),
448: PetscDesignatedInitializer(loadnative, NULL),
449: PetscDesignatedInitializer(createlocalvector, NULL),
450: PetscDesignatedInitializer(getlocalvector, NULL),
451: PetscDesignatedInitializer(restorelocalvector, NULL),
452: PetscDesignatedInitializer(getlocalvectorread, NULL),
453: PetscDesignatedInitializer(restorelocalvectorread, NULL),
454: PetscDesignatedInitializer(bindtocpu, NULL),
455: PetscDesignatedInitializer(getarraywrite, NULL),
456: PetscDesignatedInitializer(restorearraywrite, NULL),
457: PetscDesignatedInitializer(getarrayandmemtype, NULL),
458: PetscDesignatedInitializer(restorearrayandmemtype, NULL),
459: PetscDesignatedInitializer(getarrayreadandmemtype, NULL),
460: PetscDesignatedInitializer(restorearrayreadandmemtype, NULL),
461: PetscDesignatedInitializer(getarraywriteandmemtype, NULL),
462: PetscDesignatedInitializer(restorearraywriteandmemtype, NULL),
463: PetscDesignatedInitializer(concatenate, NULL),
464: PetscDesignatedInitializer(sum, NULL),
465: PetscDesignatedInitializer(setpreallocationcoo, VecSetPreallocationCOO_MPI),
466: PetscDesignatedInitializer(setvaluescoo, VecSetValuesCOO_MPI)};
468: /*
469: VecCreate_MPI_Private - Basic create routine called by VecCreate_MPI() (i.e. VecCreateMPI()),
470: VecCreateMPIWithArray(), VecCreate_Shared() (i.e. VecCreateShared()), VecCreateGhost(),
471: VecDuplicate_MPI(), VecCreateGhostWithArray(), VecDuplicate_MPI(), and VecDuplicate_Shared()
473: If alloc is true and array is NULL then this routine allocates the space, otherwise
474: no space is allocated.
475: */
476: PetscErrorCode VecCreate_MPI_Private(Vec v, PetscBool alloc, PetscInt nghost, const PetscScalar array[])
477: {
478: Vec_MPI *s;
480: PetscNew(&s);
481: v->data = (void *)s;
482: PetscMemcpy(v->ops, &DvOps, sizeof(DvOps));
483: s->nghost = nghost;
484: v->petscnative = PETSC_TRUE;
485: if (array) v->offloadmask = PETSC_OFFLOAD_CPU;
487: PetscLayoutSetUp(v->map);
489: s->array = (PetscScalar *)array;
490: s->array_allocated = NULL;
491: if (alloc && !array) {
492: PetscInt n = v->map->n + nghost;
493: PetscCalloc1(n, &s->array);
494: s->array_allocated = s->array;
495: }
497: /* By default parallel vectors do not have local representation */
498: s->localrep = NULL;
499: s->localupdate = NULL;
501: v->stash.insertmode = NOT_SET_VALUES;
502: v->bstash.insertmode = NOT_SET_VALUES;
503: /* create the stashes. The block-size for bstash is set later when
504: VecSetValuesBlocked is called.
505: */
506: VecStashCreate_Private(PetscObjectComm((PetscObject)v), 1, &v->stash);
507: VecStashCreate_Private(PetscObjectComm((PetscObject)v), PetscAbs(v->map->bs), &v->bstash);
509: #if defined(PETSC_HAVE_MATLAB)
510: PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", VecMatlabEnginePut_Default);
511: PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", VecMatlabEngineGet_Default);
512: #endif
513: PetscObjectChangeTypeName((PetscObject)v, VECMPI);
514: return 0;
515: }
517: /*MC
518: VECMPI - VECMPI = "mpi" - The basic parallel vector
520: Options Database Keys:
521: . -vec_type mpi - sets the vector type to VECMPI during a call to VecSetFromOptions()
523: Level: beginner
525: .seealso: `VecCreate()`, `VecSetType()`, `VecSetFromOptions()`, `VecCreateMPIWithArray()`, `VECMPI`, `VecType`, `VecCreateMPI()`, `VecCreateMPI()`
526: M*/
528: PetscErrorCode VecCreate_MPI(Vec vv)
529: {
530: VecCreate_MPI_Private(vv, PETSC_TRUE, 0, NULL);
531: return 0;
532: }
534: /*MC
535: VECSTANDARD = "standard" - A VECSEQ on one process and VECMPI on more than one process
537: Options Database Keys:
538: . -vec_type standard - sets a vector type to standard on calls to VecSetFromOptions()
540: Level: beginner
542: .seealso: `VecCreateSeq()`, `VecCreateMPI()`
543: M*/
545: PETSC_EXTERN PetscErrorCode VecCreate_Standard(Vec v)
546: {
547: PetscMPIInt size;
549: MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
550: if (size == 1) {
551: VecSetType(v, VECSEQ);
552: } else {
553: VecSetType(v, VECMPI);
554: }
555: return 0;
556: }
558: /*@C
559: VecCreateMPIWithArray - Creates a parallel, array-style vector,
560: where the user provides the array space to store the vector values.
562: Collective
564: Input Parameters:
565: + comm - the MPI communicator to use
566: . bs - block size, same meaning as VecSetBlockSize()
567: . n - local vector length, cannot be PETSC_DECIDE
568: . N - global vector length (or PETSC_DECIDE to have calculated)
569: - array - the user provided array to store the vector values
571: Output Parameter:
572: . vv - the vector
574: Notes:
575: Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
576: same type as an existing vector.
578: If the user-provided array is NULL, then VecPlaceArray() can be used
579: at a later stage to SET the array for storing the vector values.
581: PETSc does NOT free the array when the vector is destroyed via VecDestroy().
582: The user should not free the array until the vector is destroyed.
584: Level: intermediate
586: .seealso: `VecCreateSeqWithArray()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateGhost()`,
587: `VecCreateMPI()`, `VecCreateGhostWithArray()`, `VecPlaceArray()`
589: @*/
590: PetscErrorCode VecCreateMPIWithArray(MPI_Comm comm, PetscInt bs, PetscInt n, PetscInt N, const PetscScalar array[], Vec *vv)
591: {
593: PetscSplitOwnership(comm, &n, &N);
594: VecCreate(comm, vv);
595: VecSetSizes(*vv, n, N);
596: VecSetBlockSize(*vv, bs);
597: VecCreate_MPI_Private(*vv, PETSC_FALSE, 0, array);
598: return 0;
599: }
601: /*@C
602: VecCreateGhostWithArray - Creates a parallel vector with ghost padding on each processor;
603: the caller allocates the array space.
605: Collective
607: Input Parameters:
608: + comm - the MPI communicator to use
609: . n - local vector length
610: . N - global vector length (or PETSC_DECIDE to have calculated if n is given)
611: . nghost - number of local ghost points
612: . ghosts - global indices of ghost points (or NULL if not needed), these do not need to be in increasing order (sorted)
613: - array - the space to store the vector values (as long as n + nghost)
615: Output Parameter:
616: . vv - the global vector representation (without ghost points as part of vector)
618: Notes:
619: Use VecGhostGetLocalForm() to access the local, ghosted representation
620: of the vector.
622: This also automatically sets the ISLocalToGlobalMapping() for this vector.
624: Level: advanced
626: .seealso: `VecCreate()`, `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`,
627: `VecCreateGhost()`, `VecCreateSeqWithArray()`, `VecCreateMPIWithArray()`,
628: `VecCreateGhostBlock()`, `VecCreateGhostBlockWithArray()`, `VecMPISetGhost()`
630: @*/
631: PetscErrorCode VecCreateGhostWithArray(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[], const PetscScalar array[], Vec *vv)
632: {
633: Vec_MPI *w;
634: PetscScalar *larray;
635: IS from, to;
636: ISLocalToGlobalMapping ltog;
637: PetscInt rstart, i, *indices;
639: *vv = NULL;
644: PetscSplitOwnership(comm, &n, &N);
645: /* Create global representation */
646: VecCreate(comm, vv);
647: VecSetSizes(*vv, n, N);
648: VecCreate_MPI_Private(*vv, PETSC_TRUE, nghost, array);
649: w = (Vec_MPI *)(*vv)->data;
650: /* Create local representation */
651: VecGetArray(*vv, &larray);
652: VecCreateSeqWithArray(PETSC_COMM_SELF, 1, n + nghost, larray, &w->localrep);
653: VecRestoreArray(*vv, &larray);
655: /*
656: Create scatter context for scattering (updating) ghost values
657: */
658: ISCreateGeneral(comm, nghost, ghosts, PETSC_COPY_VALUES, &from);
659: ISCreateStride(PETSC_COMM_SELF, nghost, n, 1, &to);
660: VecScatterCreate(*vv, from, w->localrep, to, &w->localupdate);
661: ISDestroy(&to);
662: ISDestroy(&from);
664: /* set local to global mapping for ghosted vector */
665: PetscMalloc1(n + nghost, &indices);
666: VecGetOwnershipRange(*vv, &rstart, NULL);
667: for (i = 0; i < n; i++) indices[i] = rstart + i;
668: for (i = 0; i < nghost; i++) indices[n + i] = ghosts[i];
669: ISLocalToGlobalMappingCreate(comm, 1, n + nghost, indices, PETSC_OWN_POINTER, <og);
670: VecSetLocalToGlobalMapping(*vv, ltog);
671: ISLocalToGlobalMappingDestroy(<og);
672: return 0;
673: }
675: /*@
676: VecCreateGhost - Creates a parallel vector with ghost padding on each processor.
678: Collective
680: Input Parameters:
681: + comm - the MPI communicator to use
682: . n - local vector length
683: . N - global vector length (or PETSC_DECIDE to have calculated if n is given)
684: . nghost - number of local ghost points
685: - ghosts - global indices of ghost points, these do not need to be in increasing order (sorted)
687: Output Parameter:
688: . vv - the global vector representation (without ghost points as part of vector)
690: Notes:
691: Use VecGhostGetLocalForm() to access the local, ghosted representation
692: of the vector.
694: This also automatically sets the ISLocalToGlobalMapping() for this vector.
696: Level: advanced
698: .seealso: `VecCreateSeq()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateMPI()`,
699: `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`, `VecGhostUpdateBegin()`,
700: `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`, `VecGhostUpdateEnd()`,
701: `VecCreateGhostBlock()`, `VecCreateGhostBlockWithArray()`, `VecMPISetGhost()`
703: @*/
704: PetscErrorCode VecCreateGhost(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[], Vec *vv)
705: {
706: VecCreateGhostWithArray(comm, n, N, nghost, ghosts, NULL, vv);
707: return 0;
708: }
710: /*@
711: VecMPISetGhost - Sets the ghost points for an MPI ghost vector
713: Collective
715: Input Parameters:
716: + vv - the MPI vector
717: . nghost - number of local ghost points
718: - ghosts - global indices of ghost points, these do not need to be in increasing order (sorted)
720: Notes:
721: Use VecGhostGetLocalForm() to access the local, ghosted representation
722: of the vector.
724: This also automatically sets the ISLocalToGlobalMapping() for this vector.
726: You must call this AFTER you have set the type of the vector (with VecSetType()) and the size (with VecSetSizes()).
728: Level: advanced
730: .seealso: `VecCreateSeq()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateMPI()`,
731: `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`, `VecGhostUpdateBegin()`,
732: `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`, `VecGhostUpdateEnd()`,
733: `VecCreateGhostBlock()`, `VecCreateGhostBlockWithArray()`
735: @*/
736: PetscErrorCode VecMPISetGhost(Vec vv, PetscInt nghost, const PetscInt ghosts[])
737: {
738: PetscBool flg;
740: PetscObjectTypeCompare((PetscObject)vv, VECMPI, &flg);
741: /* if already fully existent VECMPI then basically destroy it and rebuild with ghosting */
742: if (flg) {
743: PetscInt n, N;
744: Vec_MPI *w;
745: PetscScalar *larray;
746: IS from, to;
747: ISLocalToGlobalMapping ltog;
748: PetscInt rstart, i, *indices;
749: MPI_Comm comm;
751: PetscObjectGetComm((PetscObject)vv, &comm);
752: n = vv->map->n;
753: N = vv->map->N;
754: PetscUseTypeMethod(vv, destroy);
755: VecSetSizes(vv, n, N);
756: VecCreate_MPI_Private(vv, PETSC_TRUE, nghost, NULL);
757: w = (Vec_MPI *)(vv)->data;
758: /* Create local representation */
759: VecGetArray(vv, &larray);
760: VecCreateSeqWithArray(PETSC_COMM_SELF, 1, n + nghost, larray, &w->localrep);
761: VecRestoreArray(vv, &larray);
763: /*
764: Create scatter context for scattering (updating) ghost values
765: */
766: ISCreateGeneral(comm, nghost, ghosts, PETSC_COPY_VALUES, &from);
767: ISCreateStride(PETSC_COMM_SELF, nghost, n, 1, &to);
768: VecScatterCreate(vv, from, w->localrep, to, &w->localupdate);
769: ISDestroy(&to);
770: ISDestroy(&from);
772: /* set local to global mapping for ghosted vector */
773: PetscMalloc1(n + nghost, &indices);
774: VecGetOwnershipRange(vv, &rstart, NULL);
776: for (i = 0; i < n; i++) indices[i] = rstart + i;
777: for (i = 0; i < nghost; i++) indices[n + i] = ghosts[i];
779: ISLocalToGlobalMappingCreate(comm, 1, n + nghost, indices, PETSC_OWN_POINTER, <og);
780: VecSetLocalToGlobalMapping(vv, ltog);
781: ISLocalToGlobalMappingDestroy(<og);
782: } else {
785: }
786: return 0;
787: }
789: /* ------------------------------------------------------------------------------------------*/
790: /*@C
791: VecCreateGhostBlockWithArray - Creates a parallel vector with ghost padding on each processor;
792: the caller allocates the array space. Indices in the ghost region are based on blocks.
794: Collective
796: Input Parameters:
797: + comm - the MPI communicator to use
798: . bs - block size
799: . n - local vector length
800: . N - global vector length (or PETSC_DECIDE to have calculated if n is given)
801: . nghost - number of local ghost blocks
802: . 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)
803: - array - the space to store the vector values (as long as n + nghost*bs)
805: Output Parameter:
806: . vv - the global vector representation (without ghost points as part of vector)
808: Notes:
809: Use VecGhostGetLocalForm() to access the local, ghosted representation
810: of the vector.
812: n is the local vector size (total local size not the number of blocks) while nghost
813: is the number of blocks in the ghost portion, i.e. the number of elements in the ghost
814: portion is bs*nghost
816: Level: advanced
818: .seealso: `VecCreate()`, `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`,
819: `VecCreateGhost()`, `VecCreateSeqWithArray()`, `VecCreateMPIWithArray()`,
820: `VecCreateGhostWithArray()`, `VecCreateGhostBlock()`
822: @*/
823: PetscErrorCode VecCreateGhostBlockWithArray(MPI_Comm comm, PetscInt bs, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[], const PetscScalar array[], Vec *vv)
824: {
825: Vec_MPI *w;
826: PetscScalar *larray;
827: IS from, to;
828: ISLocalToGlobalMapping ltog;
829: PetscInt rstart, i, nb, *indices;
831: *vv = NULL;
837: PetscSplitOwnership(comm, &n, &N);
838: /* Create global representation */
839: VecCreate(comm, vv);
840: VecSetSizes(*vv, n, N);
841: VecSetBlockSize(*vv, bs);
842: VecCreate_MPI_Private(*vv, PETSC_TRUE, nghost * bs, array);
843: w = (Vec_MPI *)(*vv)->data;
844: /* Create local representation */
845: VecGetArray(*vv, &larray);
846: VecCreateSeqWithArray(PETSC_COMM_SELF, bs, n + bs * nghost, larray, &w->localrep);
847: VecRestoreArray(*vv, &larray);
849: /*
850: Create scatter context for scattering (updating) ghost values
851: */
852: ISCreateBlock(comm, bs, nghost, ghosts, PETSC_COPY_VALUES, &from);
853: ISCreateStride(PETSC_COMM_SELF, bs * nghost, n, 1, &to);
854: VecScatterCreate(*vv, from, w->localrep, to, &w->localupdate);
855: ISDestroy(&to);
856: ISDestroy(&from);
858: /* set local to global mapping for ghosted vector */
859: nb = n / bs;
860: PetscMalloc1(nb + nghost, &indices);
861: VecGetOwnershipRange(*vv, &rstart, NULL);
862: rstart = rstart / bs;
864: for (i = 0; i < nb; i++) indices[i] = rstart + i;
865: for (i = 0; i < nghost; i++) indices[nb + i] = ghosts[i];
867: ISLocalToGlobalMappingCreate(comm, bs, nb + nghost, indices, PETSC_OWN_POINTER, <og);
868: VecSetLocalToGlobalMapping(*vv, ltog);
869: ISLocalToGlobalMappingDestroy(<og);
870: return 0;
871: }
873: /*@
874: VecCreateGhostBlock - Creates a parallel vector with ghost padding on each processor.
875: The indicing of the ghost points is done with blocks.
877: Collective
879: Input Parameters:
880: + comm - the MPI communicator to use
881: . bs - the block size
882: . n - local vector length
883: . N - global vector length (or PETSC_DECIDE to have calculated if n is given)
884: . nghost - number of local ghost blocks
885: - ghosts - global indices of ghost blocks, counts are by block, not by individual index, these do not need to be in increasing order (sorted)
887: Output Parameter:
888: . vv - the global vector representation (without ghost points as part of vector)
890: Notes:
891: Use VecGhostGetLocalForm() to access the local, ghosted representation
892: of the vector.
894: n is the local vector size (total local size not the number of blocks) while nghost
895: is the number of blocks in the ghost portion, i.e. the number of elements in the ghost
896: portion is bs*nghost
898: Level: advanced
900: .seealso: `VecCreateSeq()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateMPI()`,
901: `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`,
902: `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`, `VecCreateGhostBlockWithArray()`
904: @*/
905: PetscErrorCode VecCreateGhostBlock(MPI_Comm comm, PetscInt bs, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[], Vec *vv)
906: {
907: VecCreateGhostBlockWithArray(comm, bs, n, N, nghost, ghosts, NULL, vv);
908: return 0;
909: }