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, &ltog);
670:   VecSetLocalToGlobalMapping(*vv, ltog);
671:   ISLocalToGlobalMappingDestroy(&ltog);
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, &ltog);
780:     VecSetLocalToGlobalMapping(vv, ltog);
781:     ISLocalToGlobalMappingDestroy(&ltog);
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, &ltog);
868:   VecSetLocalToGlobalMapping(*vv, ltog);
869:   ISLocalToGlobalMappingDestroy(&ltog);
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: }