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,<og);
706: VecSetLocalToGlobalMapping(*vv,ltog);
707: ISLocalToGlobalMappingDestroy(<og);
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,<og);
823: VecSetLocalToGlobalMapping(vv,ltog);
824: ISLocalToGlobalMappingDestroy(<og);
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,<og);
913: VecSetLocalToGlobalMapping(*vv,ltog);
914: ISLocalToGlobalMappingDestroy(<og);
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: }