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;

 14:   VecDot_Seq(xin,yin,&work);
 15:   MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));
 16:   *z   = sum;
 17:   return(0);
 18: }

 20: PetscErrorCode VecTDot_MPI(Vec xin,Vec yin,PetscScalar *z)
 21: {
 22:   PetscScalar    sum,work;

 26:   VecTDot_Seq(xin,yin,&work);
 27:   MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));
 28:   *z   = sum;
 29:   return(0);
 30: }

 32: extern PetscErrorCode VecView_MPI_Draw(Vec,PetscViewer);

 34: static PetscErrorCode VecPlaceArray_MPI(Vec vin,const PetscScalar *a)
 35: {
 37:   Vec_MPI        *v = (Vec_MPI*)vin->data;

 40:   if (v->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"VecPlaceArray() was already called on this vector, without a call to VecResetArray()");
 41:   v->unplacedarray = v->array;  /* save previous array so reset can bring it back */
 42:   v->array         = (PetscScalar*)a;
 43:   if (v->localrep) {
 44:     VecPlaceArray(v->localrep,a);
 45:   }
 46:   return(0);
 47: }

 49: PetscErrorCode VecDuplicate_MPI(Vec win,Vec *v)
 50: {
 52:   Vec_MPI        *vw,*w = (Vec_MPI*)win->data;
 53:   PetscScalar    *array;

 56:   VecCreate(PetscObjectComm((PetscObject)win),v);
 57:   PetscLayoutReference(win->map,&(*v)->map);

 59:   VecCreate_MPI_Private(*v,PETSC_TRUE,w->nghost,NULL);
 60:   vw   = (Vec_MPI*)(*v)->data;
 61:   PetscMemcpy((*v)->ops,win->ops,sizeof(struct _VecOps));

 63:   /* save local representation of the parallel vector (and scatter) if it exists */
 64:   if (w->localrep) {
 65:     VecGetArray(*v,&array);
 66:     VecCreateSeqWithArray(PETSC_COMM_SELF,PetscAbs(win->map->bs),win->map->n+w->nghost,array,&vw->localrep);
 67:     PetscMemcpy(vw->localrep->ops,w->localrep->ops,sizeof(struct _VecOps));
 68:     VecRestoreArray(*v,&array);
 69:     PetscLogObjectParent((PetscObject)*v,(PetscObject)vw->localrep);

 71:     vw->localupdate = w->localupdate;
 72:     if (vw->localupdate) {
 73:       PetscObjectReference((PetscObject)vw->localupdate);
 74:     }
 75:   }

 77:   /* New vector should inherit stashing property of parent */
 78:   (*v)->stash.donotstash   = win->stash.donotstash;
 79:   (*v)->stash.ignorenegidx = win->stash.ignorenegidx;

 81:   PetscObjectListDuplicate(((PetscObject)win)->olist,&((PetscObject)(*v))->olist);
 82:   PetscFunctionListDuplicate(((PetscObject)win)->qlist,&((PetscObject)(*v))->qlist);

 84:   (*v)->map->bs   = PetscAbs(win->map->bs);
 85:   (*v)->bstash.bs = win->bstash.bs;
 86:   return(0);
 87: }

 89: static PetscErrorCode VecSetOption_MPI(Vec V,VecOption op,PetscBool flag)
 90: {
 91:   Vec_MPI        *v = (Vec_MPI*)V->data;

 95:   switch (op) {
 96:   case VEC_IGNORE_OFF_PROC_ENTRIES: V->stash.donotstash = flag;
 97:     break;
 98:   case VEC_IGNORE_NEGATIVE_INDICES: V->stash.ignorenegidx = flag;
 99:     break;
100:   case VEC_SUBSET_OFF_PROC_ENTRIES:
101:     v->assembly_subset = flag; /* See the same logic in MatAssembly wrt MAT_SUBSET_OFF_PROC_ENTRIES */
102:     if (!v->assembly_subset) { /* User indicates "do not reuse the communication pattern" */
103:       VecAssemblyReset_MPI(V); /* Reset existing pattern to free memory */
104:       v->first_assembly_done = PETSC_FALSE; /* Mark the first assembly is not done */
105:     }
106:     break;
107:   }

109:   return(0);
110: }

112: static PetscErrorCode VecResetArray_MPI(Vec vin)
113: {
114:   Vec_MPI        *v = (Vec_MPI*)vin->data;

118:   v->array         = v->unplacedarray;
119:   v->unplacedarray = NULL;
120:   if (v->localrep) {
121:     VecResetArray(v->localrep);
122:   }
123:   return(0);
124: }

126: static PetscErrorCode VecAssemblySend_MPI_Private(MPI_Comm comm,const PetscMPIInt tag[],PetscMPIInt rankid,PetscMPIInt rank,void *sdata,MPI_Request req[],void *ctx)
127: {
128:   Vec X = (Vec)ctx;
129:   Vec_MPI *x = (Vec_MPI*)X->data;
130:   VecAssemblyHeader *hdr = (VecAssemblyHeader*)sdata;
131:   PetscInt bs = X->map->bs;

135:   /* x->first_assembly_done indicates we are reusing a communication network. In that case, some
136:      messages can be empty, but we have to send them this time if we sent them before because the
137:      receiver is expecting them.
138:    */
139:   if (hdr->count || (x->first_assembly_done && x->sendptrs[rankid].ints)) {
140:     MPI_Isend(x->sendptrs[rankid].ints,hdr->count,MPIU_INT,rank,tag[0],comm,&req[0]);
141:     MPI_Isend(x->sendptrs[rankid].scalars,hdr->count,MPIU_SCALAR,rank,tag[1],comm,&req[1]);
142:   }
143:   if (hdr->bcount || (x->first_assembly_done && x->sendptrs[rankid].intb)) {
144:     MPI_Isend(x->sendptrs[rankid].intb,hdr->bcount,MPIU_INT,rank,tag[2],comm,&req[2]);
145:     MPI_Isend(x->sendptrs[rankid].scalarb,hdr->bcount*bs,MPIU_SCALAR,rank,tag[3],comm,&req[3]);
146:   }
147:   return(0);
148: }

150: static PetscErrorCode VecAssemblyRecv_MPI_Private(MPI_Comm comm,const PetscMPIInt tag[],PetscMPIInt rank,void *rdata,MPI_Request req[],void *ctx)
151: {
152:   Vec X = (Vec)ctx;
153:   Vec_MPI *x = (Vec_MPI*)X->data;
154:   VecAssemblyHeader *hdr = (VecAssemblyHeader*)rdata;
156:   PetscInt bs = X->map->bs;
157:   VecAssemblyFrame *frame;

160:   PetscSegBufferGet(x->segrecvframe,1,&frame);

162:   if (hdr->count) {
163:     PetscSegBufferGet(x->segrecvint,hdr->count,&frame->ints);
164:     MPI_Irecv(frame->ints,hdr->count,MPIU_INT,rank,tag[0],comm,&req[0]);
165:     PetscSegBufferGet(x->segrecvscalar,hdr->count,&frame->scalars);
166:     MPI_Irecv(frame->scalars,hdr->count,MPIU_SCALAR,rank,tag[1],comm,&req[1]);
167:     frame->pendings = 2;
168:   } else {
169:     frame->ints = NULL;
170:     frame->scalars = NULL;
171:     frame->pendings = 0;
172:   }

174:   if (hdr->bcount) {
175:     PetscSegBufferGet(x->segrecvint,hdr->bcount,&frame->intb);
176:     MPI_Irecv(frame->intb,hdr->bcount,MPIU_INT,rank,tag[2],comm,&req[2]);
177:     PetscSegBufferGet(x->segrecvscalar,hdr->bcount*bs,&frame->scalarb);
178:     MPI_Irecv(frame->scalarb,hdr->bcount*bs,MPIU_SCALAR,rank,tag[3],comm,&req[3]);
179:     frame->pendingb = 2;
180:   } else {
181:     frame->intb = NULL;
182:     frame->scalarb = NULL;
183:     frame->pendingb = 0;
184:   }
185:   return(0);
186: }

188: static PetscErrorCode VecAssemblyBegin_MPI_BTS(Vec X)
189: {
190:   Vec_MPI        *x = (Vec_MPI*)X->data;
192:   MPI_Comm       comm;
193:   PetscInt       i,j,jb,bs;

196:   if (X->stash.donotstash) return(0);

198:   PetscObjectGetComm((PetscObject)X,&comm);
199:   VecGetBlockSize(X,&bs);
200:   if (PetscDefined(USE_DEBUG)) {
201:     InsertMode addv;
202:     MPIU_Allreduce((PetscEnum*)&X->stash.insertmode,(PetscEnum*)&addv,1,MPIU_ENUM,MPI_BOR,comm);
203:     if (addv == (ADD_VALUES|INSERT_VALUES)) SETERRQ(comm,PETSC_ERR_ARG_NOTSAMETYPE,"Some processors inserted values while others added");
204:   }
205:   X->bstash.insertmode = X->stash.insertmode; /* Block stash implicitly tracks InsertMode of scalar stash */

207:   VecStashSortCompress_Private(&X->stash);
208:   VecStashSortCompress_Private(&X->bstash);

210:   if (!x->sendranks) {
211:     PetscMPIInt nowners,bnowners,*owners,*bowners;
212:     PetscInt ntmp;
213:     VecStashGetOwnerList_Private(&X->stash,X->map,&nowners,&owners);
214:     VecStashGetOwnerList_Private(&X->bstash,X->map,&bnowners,&bowners);
215:     PetscMergeMPIIntArray(nowners,owners,bnowners,bowners,&ntmp,&x->sendranks);
216:     x->nsendranks = ntmp;
217:     PetscFree(owners);
218:     PetscFree(bowners);
219:     PetscMalloc1(x->nsendranks,&x->sendhdr);
220:     PetscCalloc1(x->nsendranks,&x->sendptrs);
221:   }
222:   for (i=0,j=0,jb=0; i<x->nsendranks; i++) {
223:     PetscMPIInt rank = x->sendranks[i];
224:     x->sendhdr[i].insertmode = X->stash.insertmode;
225:     /* Initialize pointers for non-empty stashes the first time around.  Subsequent assemblies with
226:      * VEC_SUBSET_OFF_PROC_ENTRIES will leave the old pointers (dangling because the stash has been collected) when
227:      * there is nothing new to send, so that size-zero messages get sent instead. */
228:     x->sendhdr[i].count = 0;
229:     if (X->stash.n) {
230:       x->sendptrs[i].ints    = &X->stash.idx[j];
231:       x->sendptrs[i].scalars = &X->stash.array[j];
232:       for (; j<X->stash.n && X->stash.idx[j] < X->map->range[rank+1]; j++) x->sendhdr[i].count++;
233:     }
234:     x->sendhdr[i].bcount = 0;
235:     if (X->bstash.n) {
236:       x->sendptrs[i].intb    = &X->bstash.idx[jb];
237:       x->sendptrs[i].scalarb = &X->bstash.array[jb*bs];
238:       for (; jb<X->bstash.n && X->bstash.idx[jb]*bs < X->map->range[rank+1]; jb++) x->sendhdr[i].bcount++;
239:     }
240:   }

242:   if (!x->segrecvint) {PetscSegBufferCreate(sizeof(PetscInt),1000,&x->segrecvint);}
243:   if (!x->segrecvscalar) {PetscSegBufferCreate(sizeof(PetscScalar),1000,&x->segrecvscalar);}
244:   if (!x->segrecvframe) {PetscSegBufferCreate(sizeof(VecAssemblyFrame),50,&x->segrecvframe);}
245:   if (x->first_assembly_done) { /* this is not the first assembly */
246:     PetscMPIInt tag[4];
247:     for (i=0; i<4; i++) {PetscCommGetNewTag(comm,&tag[i]);}
248:     for (i=0; i<x->nsendranks; i++) {
249:       VecAssemblySend_MPI_Private(comm,tag,i,x->sendranks[i],x->sendhdr+i,x->sendreqs+4*i,X);
250:     }
251:     for (i=0; i<x->nrecvranks; i++) {
252:       VecAssemblyRecv_MPI_Private(comm,tag,x->recvranks[i],x->recvhdr+i,x->recvreqs+4*i,X);
253:     }
254:     x->use_status = PETSC_TRUE;
255:   } else { /* First time assembly */
256:     PetscCommBuildTwoSidedFReq(comm,3,MPIU_INT,x->nsendranks,x->sendranks,(PetscInt*)x->sendhdr,&x->nrecvranks,&x->recvranks,&x->recvhdr,4,&x->sendreqs,&x->recvreqs,VecAssemblySend_MPI_Private,VecAssemblyRecv_MPI_Private,X);
257:     x->use_status = PETSC_FALSE;
258:   }

260:   /* The first_assembly_done flag is only meaningful when x->assembly_subset is set.
261:      This line says when assembly_subset is set, then we mark that the first assembly is done.
262:    */
263:   x->first_assembly_done = x->assembly_subset;

265:   {
266:     PetscInt nstash,reallocs;
267:     VecStashGetInfo_Private(&X->stash,&nstash,&reallocs);
268:     PetscInfo2(X,"Stash has %D entries, uses %D mallocs.\n",nstash,reallocs);
269:     VecStashGetInfo_Private(&X->bstash,&nstash,&reallocs);
270:     PetscInfo2(X,"Block-Stash has %D entries, uses %D mallocs.\n",nstash,reallocs);
271:   }
272:   return(0);
273: }

275: static PetscErrorCode VecAssemblyEnd_MPI_BTS(Vec X)
276: {
277:   Vec_MPI *x = (Vec_MPI*)X->data;
278:   PetscInt bs = X->map->bs;
279:   PetscMPIInt npending,*some_indices,r;
280:   MPI_Status  *some_statuses;
281:   PetscScalar *xarray;
283:   VecAssemblyFrame *frame;

286:   if (X->stash.donotstash) {
287:     X->stash.insertmode = NOT_SET_VALUES;
288:     X->bstash.insertmode = NOT_SET_VALUES;
289:     return(0);
290:   }

292:   if (!x->segrecvframe) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing segrecvframe! Probably you forgot to call VecAssemblyBegin first");
293:   VecGetArray(X,&xarray);
294:   PetscSegBufferExtractInPlace(x->segrecvframe,&frame);
295:   PetscMalloc2(4*x->nrecvranks,&some_indices,x->use_status?4*x->nrecvranks:0,&some_statuses);
296:   for (r=0,npending=0; r<x->nrecvranks; r++) npending += frame[r].pendings + frame[r].pendingb;
297:   while (npending>0) {
298:     PetscMPIInt ndone=0,ii;
299:     /* Filling MPI_Status fields requires some resources from the MPI library.  We skip it on the first assembly, or
300:      * when VEC_SUBSET_OFF_PROC_ENTRIES has not been set, because we could exchange exact sizes in the initial
301:      * rendezvous.  When the rendezvous is elided, however, we use MPI_Status to get actual message lengths, so that
302:      * subsequent assembly can set a proper subset of the values. */
303:     MPI_Waitsome(4*x->nrecvranks,x->recvreqs,&ndone,some_indices,x->use_status?some_statuses:MPI_STATUSES_IGNORE);
304:     for (ii=0; ii<ndone; ii++) {
305:       PetscInt i = some_indices[ii]/4,j,k;
306:       InsertMode imode = (InsertMode)x->recvhdr[i].insertmode;
307:       PetscInt *recvint;
308:       PetscScalar *recvscalar;
309:       PetscBool intmsg = (PetscBool)(some_indices[ii]%2 == 0);
310:       PetscBool blockmsg = (PetscBool)((some_indices[ii]%4)/2 == 1);
311:       npending--;
312:       if (!blockmsg) { /* Scalar stash */
313:         PetscMPIInt count;
314:         if (--frame[i].pendings > 0) continue;
315:         if (x->use_status) {
316:           MPI_Get_count(&some_statuses[ii],intmsg ? MPIU_INT : MPIU_SCALAR,&count);
317:         } else count = x->recvhdr[i].count;
318:         for (j=0,recvint=frame[i].ints,recvscalar=frame[i].scalars; j<count; j++,recvint++) {
319:           PetscInt loc = *recvint - X->map->rstart;
320:           if (*recvint < X->map->rstart || X->map->rend <= *recvint) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Received vector entry %D out of local range [%D,%D)]",*recvint,X->map->rstart,X->map->rend);
321:           switch (imode) {
322:           case ADD_VALUES:
323:             xarray[loc] += *recvscalar++;
324:             break;
325:           case INSERT_VALUES:
326:             xarray[loc] = *recvscalar++;
327:             break;
328:           default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Insert mode not supported 0x%x",imode);
329:           }
330:         }
331:       } else {                  /* Block stash */
332:         PetscMPIInt count;
333:         if (--frame[i].pendingb > 0) continue;
334:         if (x->use_status) {
335:           MPI_Get_count(&some_statuses[ii],intmsg ? MPIU_INT : MPIU_SCALAR,&count);
336:           if (!intmsg) count /= bs; /* Convert from number of scalars to number of blocks */
337:         } else count = x->recvhdr[i].bcount;
338:         for (j=0,recvint=frame[i].intb,recvscalar=frame[i].scalarb; j<count; j++,recvint++) {
339:           PetscInt loc = (*recvint)*bs - X->map->rstart;
340:           switch (imode) {
341:           case ADD_VALUES:
342:             for (k=loc; k<loc+bs; k++) xarray[k] += *recvscalar++;
343:             break;
344:           case INSERT_VALUES:
345:             for (k=loc; k<loc+bs; k++) xarray[k] = *recvscalar++;
346:             break;
347:           default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Insert mode not supported 0x%x",imode);
348:           }
349:         }
350:       }
351:     }
352:   }
353:   VecRestoreArray(X,&xarray);
354:   MPI_Waitall(4*x->nsendranks,x->sendreqs,MPI_STATUSES_IGNORE);
355:   PetscFree2(some_indices,some_statuses);
356:   if (x->assembly_subset) {
357:     void *dummy;                /* reset segbuffers */
358:     PetscSegBufferExtractInPlace(x->segrecvint,&dummy);
359:     PetscSegBufferExtractInPlace(x->segrecvscalar,&dummy);
360:   } else {
361:     VecAssemblyReset_MPI(X);
362:   }

364:   X->stash.insertmode = NOT_SET_VALUES;
365:   X->bstash.insertmode = NOT_SET_VALUES;
366:   VecStashScatterEnd_Private(&X->stash);
367:   VecStashScatterEnd_Private(&X->bstash);
368:   return(0);
369: }

371: PetscErrorCode VecAssemblyReset_MPI(Vec X)
372: {
373:   Vec_MPI *x = (Vec_MPI*)X->data;

377:   PetscFree(x->sendreqs);
378:   PetscFree(x->recvreqs);
379:   PetscFree(x->sendranks);
380:   PetscFree(x->recvranks);
381:   PetscFree(x->sendhdr);
382:   PetscFree(x->recvhdr);
383:   PetscFree(x->sendptrs);
384:   PetscSegBufferDestroy(&x->segrecvint);
385:   PetscSegBufferDestroy(&x->segrecvscalar);
386:   PetscSegBufferDestroy(&x->segrecvframe);
387:   return(0);
388: }

390: static PetscErrorCode VecSetFromOptions_MPI(PetscOptionItems *PetscOptionsObject,Vec X)
391: {
392: #if !defined(PETSC_HAVE_MPIUNI)
394:   PetscBool      flg = PETSC_FALSE,set;

397:   PetscOptionsHead(PetscOptionsObject,"VecMPI Options");
398:   PetscOptionsBool("-vec_assembly_legacy","Use MPI 1 version of assembly","",flg,&flg,&set);
399:   if (set) {
400:     X->ops->assemblybegin = flg ? VecAssemblyBegin_MPI : VecAssemblyBegin_MPI_BTS;
401:     X->ops->assemblyend   = flg ? VecAssemblyEnd_MPI   : VecAssemblyEnd_MPI_BTS;
402:   }
403:   PetscOptionsTail();
404: #else
406:   X->ops->assemblybegin = VecAssemblyBegin_MPI;
407:   X->ops->assemblyend   = VecAssemblyEnd_MPI;
408: #endif
409:   return(0);
410: }

412: static struct _VecOps DvOps = { VecDuplicate_MPI, /* 1 */
413:                                 VecDuplicateVecs_Default,
414:                                 VecDestroyVecs_Default,
415:                                 VecDot_MPI,
416:                                 VecMDot_MPI,
417:                                 VecNorm_MPI,
418:                                 VecTDot_MPI,
419:                                 VecMTDot_MPI,
420:                                 VecScale_Seq,
421:                                 VecCopy_Seq, /* 10 */
422:                                 VecSet_Seq,
423:                                 VecSwap_Seq,
424:                                 VecAXPY_Seq,
425:                                 VecAXPBY_Seq,
426:                                 VecMAXPY_Seq,
427:                                 VecAYPX_Seq,
428:                                 VecWAXPY_Seq,
429:                                 VecAXPBYPCZ_Seq,
430:                                 VecPointwiseMult_Seq,
431:                                 VecPointwiseDivide_Seq,
432:                                 VecSetValues_MPI, /* 20 */
433:                                 VecAssemblyBegin_MPI_BTS,
434:                                 VecAssemblyEnd_MPI_BTS,
435:                                 NULL,
436:                                 VecGetSize_MPI,
437:                                 VecGetSize_Seq,
438:                                 NULL,
439:                                 VecMax_MPI,
440:                                 VecMin_MPI,
441:                                 VecSetRandom_Seq,
442:                                 VecSetOption_MPI,
443:                                 VecSetValuesBlocked_MPI,
444:                                 VecDestroy_MPI,
445:                                 VecView_MPI,
446:                                 VecPlaceArray_MPI,
447:                                 VecReplaceArray_Seq,
448:                                 VecDot_Seq,
449:                                 VecTDot_Seq,
450:                                 VecNorm_Seq,
451:                                 VecMDot_Seq,
452:                                 VecMTDot_Seq,
453:                                 VecLoad_Default,
454:                                 VecReciprocal_Default,
455:                                 VecConjugate_Seq,
456:                                 NULL,
457:                                 NULL,
458:                                 VecResetArray_MPI,
459:                                 VecSetFromOptions_MPI,/*set from options */
460:                                 VecMaxPointwiseDivide_Seq,
461:                                 VecPointwiseMax_Seq,
462:                                 VecPointwiseMaxAbs_Seq,
463:                                 VecPointwiseMin_Seq,
464:                                 VecGetValues_MPI,
465:                                 NULL,
466:                                 NULL,
467:                                 NULL,
468:                                 NULL,
469:                                 NULL,
470:                                 NULL,
471:                                 VecStrideGather_Default,
472:                                 VecStrideScatter_Default,
473:                                 NULL,
474:                                 NULL,
475:                                 NULL,
476:                                 NULL,
477:                                 NULL,
478:                                 VecStrideSubSetGather_Default,
479:                                 VecStrideSubSetScatter_Default,
480:                                 NULL,
481:                                 NULL,
482:                                 NULL
483: };

485: /*
486:     VecCreate_MPI_Private - Basic create routine called by VecCreate_MPI() (i.e. VecCreateMPI()),
487:     VecCreateMPIWithArray(), VecCreate_Shared() (i.e. VecCreateShared()), VecCreateGhost(),
488:     VecDuplicate_MPI(), VecCreateGhostWithArray(), VecDuplicate_MPI(), and VecDuplicate_Shared()

490:     If alloc is true and array is NULL then this routine allocates the space, otherwise
491:     no space is allocated.
492: */
493: PetscErrorCode VecCreate_MPI_Private(Vec v,PetscBool alloc,PetscInt nghost,const PetscScalar array[])
494: {
495:   Vec_MPI        *s;

499:   PetscNewLog(v,&s);
500:   v->data        = (void*)s;
501:   PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));
502:   s->nghost      = nghost;
503:   v->petscnative = PETSC_TRUE;
504:   if (array) v->offloadmask = PETSC_OFFLOAD_CPU;

506:   PetscLayoutSetUp(v->map);

508:   s->array           = (PetscScalar*)array;
509:   s->array_allocated = NULL;
510:   if (alloc && !array) {
511:     PetscInt n = v->map->n+nghost;
512:     PetscCalloc1(n,&s->array);
513:     PetscLogObjectMemory((PetscObject)v,n*sizeof(PetscScalar));
514:     s->array_allocated = s->array;
515:   }

517:   /* By default parallel vectors do not have local representation */
518:   s->localrep    = NULL;
519:   s->localupdate = NULL;

521:   v->stash.insertmode = NOT_SET_VALUES;
522:   v->bstash.insertmode = NOT_SET_VALUES;
523:   /* create the stashes. The block-size for bstash is set later when
524:      VecSetValuesBlocked is called.
525:   */
526:   VecStashCreate_Private(PetscObjectComm((PetscObject)v),1,&v->stash);
527:   VecStashCreate_Private(PetscObjectComm((PetscObject)v),PetscAbs(v->map->bs),&v->bstash);

529: #if defined(PETSC_HAVE_MATLAB_ENGINE)
530:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);
531:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);
532: #endif
533:   PetscObjectChangeTypeName((PetscObject)v,VECMPI);
534:   return(0);
535: }

537: /*MC
538:    VECMPI - VECMPI = "mpi" - The basic parallel vector

540:    Options Database Keys:
541: . -vec_type mpi - sets the vector type to VECMPI during a call to VecSetFromOptions()

543:   Level: beginner

545: .seealso: VecCreate(), VecSetType(), VecSetFromOptions(), VecCreateMPIWithArray(), VECMPI, VecType, VecCreateMPI(), VecCreateMPI()
546: M*/

548: PetscErrorCode VecCreate_MPI(Vec vv)
549: {

553:   VecCreate_MPI_Private(vv,PETSC_TRUE,0,NULL);
554:   return(0);
555: }

557: /*MC
558:    VECSTANDARD = "standard" - A VECSEQ on one process and VECMPI on more than one process

560:    Options Database Keys:
561: . -vec_type standard - sets a vector type to standard on calls to VecSetFromOptions()

563:   Level: beginner

565: .seealso: VecCreateSeq(), VecCreateMPI()
566: M*/

568: PETSC_EXTERN PetscErrorCode VecCreate_Standard(Vec v)
569: {
571:   PetscMPIInt    size;

574:   MPI_Comm_size(PetscObjectComm((PetscObject)v),&size);
575:   if (size == 1) {
576:     VecSetType(v,VECSEQ);
577:   } else {
578:     VecSetType(v,VECMPI);
579:   }
580:   return(0);
581: }

583: /*@C
584:    VecCreateMPIWithArray - Creates a parallel, array-style vector,
585:    where the user provides the array space to store the vector values.

587:    Collective

589:    Input Parameters:
590: +  comm  - the MPI communicator to use
591: .  bs    - block size, same meaning as VecSetBlockSize()
592: .  n     - local vector length, cannot be PETSC_DECIDE
593: .  N     - global vector length (or PETSC_DECIDE to have calculated)
594: -  array - the user provided array to store the vector values

596:    Output Parameter:
597: .  vv - the vector

599:    Notes:
600:    Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
601:    same type as an existing vector.

603:    If the user-provided array is NULL, then VecPlaceArray() can be used
604:    at a later stage to SET the array for storing the vector values.

606:    PETSc does NOT free the array when the vector is destroyed via VecDestroy().
607:    The user should not free the array until the vector is destroyed.

609:    Level: intermediate

611: .seealso: VecCreateSeqWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateGhost(),
612:           VecCreateMPI(), VecCreateGhostWithArray(), VecPlaceArray()

614: @*/
615: PetscErrorCode  VecCreateMPIWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,const PetscScalar array[],Vec *vv)
616: {

620:   if (n == PETSC_DECIDE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must set local size of vector");
621:   PetscSplitOwnership(comm,&n,&N);
622:   VecCreate(comm,vv);
623:   VecSetSizes(*vv,n,N);
624:   VecSetBlockSize(*vv,bs);
625:   VecCreate_MPI_Private(*vv,PETSC_FALSE,0,array);
626:   return(0);
627: }

629: /*@C
630:    VecCreateGhostWithArray - Creates a parallel vector with ghost padding on each processor;
631:    the caller allocates the array space.

633:    Collective

635:    Input Parameters:
636: +  comm - the MPI communicator to use
637: .  n - local vector length
638: .  N - global vector length (or PETSC_DECIDE to have calculated if n is given)
639: .  nghost - number of local ghost points
640: .  ghosts - global indices of ghost points (or NULL if not needed), these do not need to be in increasing order (sorted)
641: -  array - the space to store the vector values (as long as n + nghost)

643:    Output Parameter:
644: .  vv - the global vector representation (without ghost points as part of vector)

646:    Notes:
647:    Use VecGhostGetLocalForm() to access the local, ghosted representation
648:    of the vector.

650:    This also automatically sets the ISLocalToGlobalMapping() for this vector.

652:    Level: advanced

654: .seealso: VecCreate(), VecGhostGetLocalForm(), VecGhostRestoreLocalForm(),
655:           VecCreateGhost(), VecCreateSeqWithArray(), VecCreateMPIWithArray(),
656:           VecCreateGhostBlock(), VecCreateGhostBlockWithArray(), VecMPISetGhost()

658: @*/
659: PetscErrorCode  VecCreateGhostWithArray(MPI_Comm comm,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],const PetscScalar array[],Vec *vv)
660: {
661:   PetscErrorCode         ierr;
662:   Vec_MPI                *w;
663:   PetscScalar            *larray;
664:   IS                     from,to;
665:   ISLocalToGlobalMapping ltog;
666:   PetscInt               rstart,i,*indices;

669:   *vv = NULL;

671:   if (n == PETSC_DECIDE)      SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must set local size");
672:   if (nghost == PETSC_DECIDE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must set local ghost size");
673:   if (nghost < 0)             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Ghost length must be >= 0");
674:   PetscSplitOwnership(comm,&n,&N);
675:   /* Create global representation */
676:   VecCreate(comm,vv);
677:   VecSetSizes(*vv,n,N);
678:   VecCreate_MPI_Private(*vv,PETSC_TRUE,nghost,array);
679:   w    = (Vec_MPI*)(*vv)->data;
680:   /* Create local representation */
681:   VecGetArray(*vv,&larray);
682:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,n+nghost,larray,&w->localrep);
683:   PetscLogObjectParent((PetscObject)*vv,(PetscObject)w->localrep);
684:   VecRestoreArray(*vv,&larray);

686:   /*
687:        Create scatter context for scattering (updating) ghost values
688:   */
689:   ISCreateGeneral(comm,nghost,ghosts,PETSC_COPY_VALUES,&from);
690:   ISCreateStride(PETSC_COMM_SELF,nghost,n,1,&to);
691:   VecScatterCreate(*vv,from,w->localrep,to,&w->localupdate);
692:   PetscLogObjectParent((PetscObject)*vv,(PetscObject)w->localupdate);
693:   ISDestroy(&to);
694:   ISDestroy(&from);

696:   /* set local to global mapping for ghosted vector */
697:   PetscMalloc1(n+nghost,&indices);
698:   VecGetOwnershipRange(*vv,&rstart,NULL);
699:   for (i=0; i<n; i++) {
700:     indices[i] = rstart + i;
701:   }
702:   for (i=0; i<nghost; i++) {
703:     indices[n+i] = ghosts[i];
704:   }
705:   ISLocalToGlobalMappingCreate(comm,1,n+nghost,indices,PETSC_OWN_POINTER,&ltog);
706:   VecSetLocalToGlobalMapping(*vv,ltog);
707:   ISLocalToGlobalMappingDestroy(&ltog);
708:   return(0);
709: }

711: /*@
712:    VecCreateGhost - Creates a parallel vector with ghost padding on each processor.

714:    Collective

716:    Input Parameters:
717: +  comm - the MPI communicator to use
718: .  n - local vector length
719: .  N - global vector length (or PETSC_DECIDE to have calculated if n is given)
720: .  nghost - number of local ghost points
721: -  ghosts - global indices of ghost points, these do not need to be in increasing order (sorted)

723:    Output Parameter:
724: .  vv - the global vector representation (without ghost points as part of vector)

726:    Notes:
727:    Use VecGhostGetLocalForm() to access the local, ghosted representation
728:    of the vector.

730:    This also automatically sets the ISLocalToGlobalMapping() for this vector.

732:    Level: advanced

734: .seealso: VecCreateSeq(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateMPI(),
735:           VecGhostGetLocalForm(), VecGhostRestoreLocalForm(), VecGhostUpdateBegin(),
736:           VecCreateGhostWithArray(), VecCreateMPIWithArray(), VecGhostUpdateEnd(),
737:           VecCreateGhostBlock(), VecCreateGhostBlockWithArray(), VecMPISetGhost()

739: @*/
740: PetscErrorCode  VecCreateGhost(MPI_Comm comm,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],Vec *vv)
741: {

745:   VecCreateGhostWithArray(comm,n,N,nghost,ghosts,NULL,vv);
746:   return(0);
747: }

749: /*@
750:    VecMPISetGhost - Sets the ghost points for an MPI ghost vector

752:    Collective on Vec

754:    Input Parameters:
755: +  vv - the MPI vector
756: .  nghost - number of local ghost points
757: -  ghosts - global indices of ghost points, these do not need to be in increasing order (sorted)

759:    Notes:
760:    Use VecGhostGetLocalForm() to access the local, ghosted representation
761:    of the vector.

763:    This also automatically sets the ISLocalToGlobalMapping() for this vector.

765:    You must call this AFTER you have set the type of the vector (with VecSetType()) and the size (with VecSetSizes()).

767:    Level: advanced

769: .seealso: VecCreateSeq(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateMPI(),
770:           VecGhostGetLocalForm(), VecGhostRestoreLocalForm(), VecGhostUpdateBegin(),
771:           VecCreateGhostWithArray(), VecCreateMPIWithArray(), VecGhostUpdateEnd(),
772:           VecCreateGhostBlock(), VecCreateGhostBlockWithArray()

774: @*/
775: PetscErrorCode  VecMPISetGhost(Vec vv,PetscInt nghost,const PetscInt ghosts[])
776: {
778:   PetscBool      flg;

781:   PetscObjectTypeCompare((PetscObject)vv,VECMPI,&flg);
782:   /* if already fully existant VECMPI then basically destroy it and rebuild with ghosting */
783:   if (flg) {
784:     PetscInt               n,N;
785:     Vec_MPI                *w;
786:     PetscScalar            *larray;
787:     IS                     from,to;
788:     ISLocalToGlobalMapping ltog;
789:     PetscInt               rstart,i,*indices;
790:     MPI_Comm               comm;

792:     PetscObjectGetComm((PetscObject)vv,&comm);
793:     n    = vv->map->n;
794:     N    = vv->map->N;
795:     (*vv->ops->destroy)(vv);
796:     VecSetSizes(vv,n,N);
797:     VecCreate_MPI_Private(vv,PETSC_TRUE,nghost,NULL);
798:     w    = (Vec_MPI*)(vv)->data;
799:     /* Create local representation */
800:     VecGetArray(vv,&larray);
801:     VecCreateSeqWithArray(PETSC_COMM_SELF,1,n+nghost,larray,&w->localrep);
802:     PetscLogObjectParent((PetscObject)vv,(PetscObject)w->localrep);
803:     VecRestoreArray(vv,&larray);

805:     /*
806:      Create scatter context for scattering (updating) ghost values
807:      */
808:     ISCreateGeneral(comm,nghost,ghosts,PETSC_COPY_VALUES,&from);
809:     ISCreateStride(PETSC_COMM_SELF,nghost,n,1,&to);
810:     VecScatterCreate(vv,from,w->localrep,to,&w->localupdate);
811:     PetscLogObjectParent((PetscObject)vv,(PetscObject)w->localupdate);
812:     ISDestroy(&to);
813:     ISDestroy(&from);

815:     /* set local to global mapping for ghosted vector */
816:     PetscMalloc1(n+nghost,&indices);
817:     VecGetOwnershipRange(vv,&rstart,NULL);

819:     for (i=0; i<n; i++)      indices[i]   = rstart + i;
820:     for (i=0; i<nghost; i++) indices[n+i] = ghosts[i];

822:     ISLocalToGlobalMappingCreate(comm,1,n+nghost,indices,PETSC_OWN_POINTER,&ltog);
823:     VecSetLocalToGlobalMapping(vv,ltog);
824:     ISLocalToGlobalMappingDestroy(&ltog);
825:   } else if (vv->ops->create == VecCreate_MPI) SETERRQ(PetscObjectComm((PetscObject)vv),PETSC_ERR_ARG_WRONGSTATE,"Must set local or global size before setting ghosting");
826:   else if (!((PetscObject)vv)->type_name) SETERRQ(PetscObjectComm((PetscObject)vv),PETSC_ERR_ARG_WRONGSTATE,"Must set type to VECMPI before ghosting");
827:   return(0);
828: }

830: /* ------------------------------------------------------------------------------------------*/
831: /*@C
832:    VecCreateGhostBlockWithArray - Creates a parallel vector with ghost padding on each processor;
833:    the caller allocates the array space. Indices in the ghost region are based on blocks.

835:    Collective

837:    Input Parameters:
838: +  comm - the MPI communicator to use
839: .  bs - block size
840: .  n - local vector length
841: .  N - global vector length (or PETSC_DECIDE to have calculated if n is given)
842: .  nghost - number of local ghost blocks
843: .  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)
844: -  array - the space to store the vector values (as long as n + nghost*bs)

846:    Output Parameter:
847: .  vv - the global vector representation (without ghost points as part of vector)

849:    Notes:
850:    Use VecGhostGetLocalForm() to access the local, ghosted representation
851:    of the vector.

853:    n is the local vector size (total local size not the number of blocks) while nghost
854:    is the number of blocks in the ghost portion, i.e. the number of elements in the ghost
855:    portion is bs*nghost

857:    Level: advanced

859: .seealso: VecCreate(), VecGhostGetLocalForm(), VecGhostRestoreLocalForm(),
860:           VecCreateGhost(), VecCreateSeqWithArray(), VecCreateMPIWithArray(),
861:           VecCreateGhostWithArray(), VecCreateGhostBlock()

863: @*/
864: PetscErrorCode  VecCreateGhostBlockWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],const PetscScalar array[],Vec *vv)
865: {
866:   PetscErrorCode         ierr;
867:   Vec_MPI                *w;
868:   PetscScalar            *larray;
869:   IS                     from,to;
870:   ISLocalToGlobalMapping ltog;
871:   PetscInt               rstart,i,nb,*indices;

874:   *vv = NULL;

876:   if (n == PETSC_DECIDE)      SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must set local size");
877:   if (nghost == PETSC_DECIDE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must set local ghost size");
878:   if (nghost < 0)             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Ghost length must be >= 0");
879:   if (n % bs)                 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local size must be a multiple of block size");
880:   PetscSplitOwnership(comm,&n,&N);
881:   /* Create global representation */
882:   VecCreate(comm,vv);
883:   VecSetSizes(*vv,n,N);
884:   VecSetBlockSize(*vv,bs);
885:   VecCreate_MPI_Private(*vv,PETSC_TRUE,nghost*bs,array);
886:   w    = (Vec_MPI*)(*vv)->data;
887:   /* Create local representation */
888:   VecGetArray(*vv,&larray);
889:   VecCreateSeqWithArray(PETSC_COMM_SELF,bs,n+bs*nghost,larray,&w->localrep);
890:   PetscLogObjectParent((PetscObject)*vv,(PetscObject)w->localrep);
891:   VecRestoreArray(*vv,&larray);

893:   /*
894:        Create scatter context for scattering (updating) ghost values
895:   */
896:   ISCreateBlock(comm,bs,nghost,ghosts,PETSC_COPY_VALUES,&from);
897:   ISCreateStride(PETSC_COMM_SELF,bs*nghost,n,1,&to);
898:   VecScatterCreate(*vv,from,w->localrep,to,&w->localupdate);
899:   PetscLogObjectParent((PetscObject)*vv,(PetscObject)w->localupdate);
900:   ISDestroy(&to);
901:   ISDestroy(&from);

903:   /* set local to global mapping for ghosted vector */
904:   nb     = n/bs;
905:   PetscMalloc1(nb+nghost,&indices);
906:   VecGetOwnershipRange(*vv,&rstart,NULL);
907:   rstart = rstart/bs;

909:   for (i=0; i<nb; i++)      indices[i]    = rstart + i;
910:   for (i=0; i<nghost; i++)  indices[nb+i] = ghosts[i];

912:   ISLocalToGlobalMappingCreate(comm,bs,nb+nghost,indices,PETSC_OWN_POINTER,&ltog);
913:   VecSetLocalToGlobalMapping(*vv,ltog);
914:   ISLocalToGlobalMappingDestroy(&ltog);
915:   return(0);
916: }

918: /*@
919:    VecCreateGhostBlock - Creates a parallel vector with ghost padding on each processor.
920:         The indicing of the ghost points is done with blocks.

922:    Collective

924:    Input Parameters:
925: +  comm - the MPI communicator to use
926: .  bs - the block size
927: .  n - local vector length
928: .  N - global vector length (or PETSC_DECIDE to have calculated if n is given)
929: .  nghost - number of local ghost blocks
930: -  ghosts - global indices of ghost blocks, counts are by block, not by individual index, these do not need to be in increasing order (sorted)

932:    Output Parameter:
933: .  vv - the global vector representation (without ghost points as part of vector)

935:    Notes:
936:    Use VecGhostGetLocalForm() to access the local, ghosted representation
937:    of the vector.

939:    n is the local vector size (total local size not the number of blocks) while nghost
940:    is the number of blocks in the ghost portion, i.e. the number of elements in the ghost
941:    portion is bs*nghost

943:    Level: advanced

945: .seealso: VecCreateSeq(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateMPI(),
946:           VecGhostGetLocalForm(), VecGhostRestoreLocalForm(),
947:           VecCreateGhostWithArray(), VecCreateMPIWithArray(), VecCreateGhostBlockWithArray()

949: @*/
950: PetscErrorCode  VecCreateGhostBlock(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],Vec *vv)
951: {

955:   VecCreateGhostBlockWithArray(comm,bs,n,N,nghost,ghosts,NULL,vv);
956:   return(0);
957: }