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) {
38: VecPlaceArray(v->localrep,a);
39: }
40: return 0;
41: }
43: PetscErrorCode VecDuplicate_MPI(Vec win,Vec *v)
44: {
45: Vec_MPI *vw,*w = (Vec_MPI*)win->data;
46: PetscScalar *array;
48: VecCreate(PetscObjectComm((PetscObject)win),v);
49: PetscLayoutReference(win->map,&(*v)->map);
51: VecCreate_MPI_Private(*v,PETSC_TRUE,w->nghost,NULL);
52: vw = (Vec_MPI*)(*v)->data;
53: PetscMemcpy((*v)->ops,win->ops,sizeof(struct _VecOps));
55: /* save local representation of the parallel vector (and scatter) if it exists */
56: if (w->localrep) {
57: VecGetArray(*v,&array);
58: VecCreateSeqWithArray(PETSC_COMM_SELF,PetscAbs(win->map->bs),win->map->n+w->nghost,array,&vw->localrep);
59: PetscMemcpy(vw->localrep->ops,w->localrep->ops,sizeof(struct _VecOps));
60: VecRestoreArray(*v,&array);
61: PetscLogObjectParent((PetscObject)*v,(PetscObject)vw->localrep);
63: vw->localupdate = w->localupdate;
64: if (vw->localupdate) {
65: PetscObjectReference((PetscObject)vw->localupdate);
66: }
67: }
69: /* New vector should inherit stashing property of parent */
70: (*v)->stash.donotstash = win->stash.donotstash;
71: (*v)->stash.ignorenegidx = win->stash.ignorenegidx;
73: PetscObjectListDuplicate(((PetscObject)win)->olist,&((PetscObject)(*v))->olist);
74: PetscFunctionListDuplicate(((PetscObject)win)->qlist,&((PetscObject)(*v))->qlist);
76: (*v)->map->bs = PetscAbs(win->map->bs);
77: (*v)->bstash.bs = win->bstash.bs;
78: return 0;
79: }
81: static PetscErrorCode VecSetOption_MPI(Vec V,VecOption op,PetscBool flag)
82: {
83: Vec_MPI *v = (Vec_MPI*)V->data;
85: switch (op) {
86: case VEC_IGNORE_OFF_PROC_ENTRIES: V->stash.donotstash = flag;
87: break;
88: case VEC_IGNORE_NEGATIVE_INDICES: V->stash.ignorenegidx = flag;
89: break;
90: case VEC_SUBSET_OFF_PROC_ENTRIES:
91: v->assembly_subset = flag; /* See the same logic in MatAssembly wrt MAT_SUBSET_OFF_PROC_ENTRIES */
92: if (!v->assembly_subset) { /* User indicates "do not reuse the communication pattern" */
93: VecAssemblyReset_MPI(V); /* Reset existing pattern to free memory */
94: v->first_assembly_done = PETSC_FALSE; /* Mark the first assembly is not done */
95: }
96: break;
97: }
99: return 0;
100: }
102: static PetscErrorCode VecResetArray_MPI(Vec vin)
103: {
104: Vec_MPI *v = (Vec_MPI*)vin->data;
106: v->array = v->unplacedarray;
107: v->unplacedarray = NULL;
108: if (v->localrep) {
109: VecResetArray(v->localrep);
110: }
111: return 0;
112: }
114: static PetscErrorCode VecAssemblySend_MPI_Private(MPI_Comm comm,const PetscMPIInt tag[],PetscMPIInt rankid,PetscMPIInt rank,void *sdata,MPI_Request req[],void *ctx)
115: {
116: Vec X = (Vec)ctx;
117: Vec_MPI *x = (Vec_MPI*)X->data;
118: VecAssemblyHeader *hdr = (VecAssemblyHeader*)sdata;
119: PetscInt bs = X->map->bs;
121: /* x->first_assembly_done indicates we are reusing a communication network. In that case, some
122: messages can be empty, but we have to send them this time if we sent them before because the
123: receiver is expecting them.
124: */
125: if (hdr->count || (x->first_assembly_done && x->sendptrs[rankid].ints)) {
126: MPI_Isend(x->sendptrs[rankid].ints,hdr->count,MPIU_INT,rank,tag[0],comm,&req[0]);
127: MPI_Isend(x->sendptrs[rankid].scalars,hdr->count,MPIU_SCALAR,rank,tag[1],comm,&req[1]);
128: }
129: if (hdr->bcount || (x->first_assembly_done && x->sendptrs[rankid].intb)) {
130: MPI_Isend(x->sendptrs[rankid].intb,hdr->bcount,MPIU_INT,rank,tag[2],comm,&req[2]);
131: MPI_Isend(x->sendptrs[rankid].scalarb,hdr->bcount*bs,MPIU_SCALAR,rank,tag[3],comm,&req[3]);
132: }
133: return 0;
134: }
136: static PetscErrorCode VecAssemblyRecv_MPI_Private(MPI_Comm comm,const PetscMPIInt tag[],PetscMPIInt rank,void *rdata,MPI_Request req[],void *ctx)
137: {
138: Vec X = (Vec)ctx;
139: Vec_MPI *x = (Vec_MPI*)X->data;
140: VecAssemblyHeader *hdr = (VecAssemblyHeader*)rdata;
141: PetscInt bs = X->map->bs;
142: VecAssemblyFrame *frame;
144: PetscSegBufferGet(x->segrecvframe,1,&frame);
146: if (hdr->count) {
147: PetscSegBufferGet(x->segrecvint,hdr->count,&frame->ints);
148: MPI_Irecv(frame->ints,hdr->count,MPIU_INT,rank,tag[0],comm,&req[0]);
149: PetscSegBufferGet(x->segrecvscalar,hdr->count,&frame->scalars);
150: MPI_Irecv(frame->scalars,hdr->count,MPIU_SCALAR,rank,tag[1],comm,&req[1]);
151: frame->pendings = 2;
152: } else {
153: frame->ints = NULL;
154: frame->scalars = NULL;
155: frame->pendings = 0;
156: }
158: if (hdr->bcount) {
159: PetscSegBufferGet(x->segrecvint,hdr->bcount,&frame->intb);
160: MPI_Irecv(frame->intb,hdr->bcount,MPIU_INT,rank,tag[2],comm,&req[2]);
161: PetscSegBufferGet(x->segrecvscalar,hdr->bcount*bs,&frame->scalarb);
162: MPI_Irecv(frame->scalarb,hdr->bcount*bs,MPIU_SCALAR,rank,tag[3],comm,&req[3]);
163: frame->pendingb = 2;
164: } else {
165: frame->intb = NULL;
166: frame->scalarb = NULL;
167: frame->pendingb = 0;
168: }
169: return 0;
170: }
172: static PetscErrorCode VecAssemblyBegin_MPI_BTS(Vec X)
173: {
174: Vec_MPI *x = (Vec_MPI*)X->data;
175: MPI_Comm comm;
176: PetscInt i,j,jb,bs;
178: if (X->stash.donotstash) return 0;
180: PetscObjectGetComm((PetscObject)X,&comm);
181: VecGetBlockSize(X,&bs);
182: if (PetscDefined(USE_DEBUG)) {
183: InsertMode addv;
184: MPIU_Allreduce((PetscEnum*)&X->stash.insertmode,(PetscEnum*)&addv,1,MPIU_ENUM,MPI_BOR,comm);
186: }
187: X->bstash.insertmode = X->stash.insertmode; /* Block stash implicitly tracks InsertMode of scalar stash */
189: VecStashSortCompress_Private(&X->stash);
190: VecStashSortCompress_Private(&X->bstash);
192: if (!x->sendranks) {
193: PetscMPIInt nowners,bnowners,*owners,*bowners;
194: PetscInt ntmp;
195: VecStashGetOwnerList_Private(&X->stash,X->map,&nowners,&owners);
196: VecStashGetOwnerList_Private(&X->bstash,X->map,&bnowners,&bowners);
197: PetscMergeMPIIntArray(nowners,owners,bnowners,bowners,&ntmp,&x->sendranks);
198: x->nsendranks = ntmp;
199: PetscFree(owners);
200: PetscFree(bowners);
201: PetscMalloc1(x->nsendranks,&x->sendhdr);
202: PetscCalloc1(x->nsendranks,&x->sendptrs);
203: }
204: for (i=0,j=0,jb=0; i<x->nsendranks; i++) {
205: PetscMPIInt rank = x->sendranks[i];
206: x->sendhdr[i].insertmode = X->stash.insertmode;
207: /* Initialize pointers for non-empty stashes the first time around. Subsequent assemblies with
208: * VEC_SUBSET_OFF_PROC_ENTRIES will leave the old pointers (dangling because the stash has been collected) when
209: * there is nothing new to send, so that size-zero messages get sent instead. */
210: x->sendhdr[i].count = 0;
211: if (X->stash.n) {
212: x->sendptrs[i].ints = &X->stash.idx[j];
213: x->sendptrs[i].scalars = &X->stash.array[j];
214: for (; j<X->stash.n && X->stash.idx[j] < X->map->range[rank+1]; j++) x->sendhdr[i].count++;
215: }
216: x->sendhdr[i].bcount = 0;
217: if (X->bstash.n) {
218: x->sendptrs[i].intb = &X->bstash.idx[jb];
219: x->sendptrs[i].scalarb = &X->bstash.array[jb*bs];
220: for (; jb<X->bstash.n && X->bstash.idx[jb]*bs < X->map->range[rank+1]; jb++) x->sendhdr[i].bcount++;
221: }
222: }
224: if (!x->segrecvint) PetscSegBufferCreate(sizeof(PetscInt),1000,&x->segrecvint);
225: if (!x->segrecvscalar) PetscSegBufferCreate(sizeof(PetscScalar),1000,&x->segrecvscalar);
226: if (!x->segrecvframe) PetscSegBufferCreate(sizeof(VecAssemblyFrame),50,&x->segrecvframe);
227: if (x->first_assembly_done) { /* this is not the first assembly */
228: PetscMPIInt tag[4];
229: for (i=0; i<4; i++) PetscCommGetNewTag(comm,&tag[i]);
230: for (i=0; i<x->nsendranks; i++) {
231: VecAssemblySend_MPI_Private(comm,tag,i,x->sendranks[i],x->sendhdr+i,x->sendreqs+4*i,X);
232: }
233: for (i=0; i<x->nrecvranks; i++) {
234: VecAssemblyRecv_MPI_Private(comm,tag,x->recvranks[i],x->recvhdr+i,x->recvreqs+4*i,X);
235: }
236: x->use_status = PETSC_TRUE;
237: } else { /* First time assembly */
238: 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);
239: x->use_status = PETSC_FALSE;
240: }
242: /* The first_assembly_done flag is only meaningful when x->assembly_subset is set.
243: This line says when assembly_subset is set, then we mark that the first assembly is done.
244: */
245: x->first_assembly_done = x->assembly_subset;
247: {
248: PetscInt nstash,reallocs;
249: VecStashGetInfo_Private(&X->stash,&nstash,&reallocs);
250: PetscInfo(X,"Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n",nstash,reallocs);
251: VecStashGetInfo_Private(&X->bstash,&nstash,&reallocs);
252: PetscInfo(X,"Block-Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n",nstash,reallocs);
253: }
254: return 0;
255: }
257: static PetscErrorCode VecAssemblyEnd_MPI_BTS(Vec X)
258: {
259: Vec_MPI *x = (Vec_MPI*)X->data;
260: PetscInt bs = X->map->bs;
261: PetscMPIInt npending,*some_indices,r;
262: MPI_Status *some_statuses;
263: PetscScalar *xarray;
264: VecAssemblyFrame *frame;
266: if (X->stash.donotstash) {
267: X->stash.insertmode = NOT_SET_VALUES;
268: X->bstash.insertmode = NOT_SET_VALUES;
269: return 0;
270: }
273: VecGetArray(X,&xarray);
274: PetscSegBufferExtractInPlace(x->segrecvframe,&frame);
275: PetscMalloc2(4*x->nrecvranks,&some_indices,x->use_status?4*x->nrecvranks:0,&some_statuses);
276: for (r=0,npending=0; r<x->nrecvranks; r++) npending += frame[r].pendings + frame[r].pendingb;
277: while (npending>0) {
278: PetscMPIInt ndone=0,ii;
279: /* Filling MPI_Status fields requires some resources from the MPI library. We skip it on the first assembly, or
280: * when VEC_SUBSET_OFF_PROC_ENTRIES has not been set, because we could exchange exact sizes in the initial
281: * rendezvous. When the rendezvous is elided, however, we use MPI_Status to get actual message lengths, so that
282: * subsequent assembly can set a proper subset of the values. */
283: MPI_Waitsome(4*x->nrecvranks,x->recvreqs,&ndone,some_indices,x->use_status?some_statuses:MPI_STATUSES_IGNORE);
284: for (ii=0; ii<ndone; ii++) {
285: PetscInt i = some_indices[ii]/4,j,k;
286: InsertMode imode = (InsertMode)x->recvhdr[i].insertmode;
287: PetscInt *recvint;
288: PetscScalar *recvscalar;
289: PetscBool intmsg = (PetscBool)(some_indices[ii]%2 == 0);
290: PetscBool blockmsg = (PetscBool)((some_indices[ii]%4)/2 == 1);
291: npending--;
292: if (!blockmsg) { /* Scalar stash */
293: PetscMPIInt count;
294: if (--frame[i].pendings > 0) continue;
295: if (x->use_status) {
296: MPI_Get_count(&some_statuses[ii],intmsg ? MPIU_INT : MPIU_SCALAR,&count);
297: } else count = x->recvhdr[i].count;
298: for (j=0,recvint=frame[i].ints,recvscalar=frame[i].scalars; j<count; j++,recvint++) {
299: PetscInt loc = *recvint - X->map->rstart;
301: switch (imode) {
302: case ADD_VALUES:
303: xarray[loc] += *recvscalar++;
304: break;
305: case INSERT_VALUES:
306: xarray[loc] = *recvscalar++;
307: break;
308: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Insert mode not supported 0x%x",imode);
309: }
310: }
311: } else { /* Block stash */
312: PetscMPIInt count;
313: if (--frame[i].pendingb > 0) continue;
314: if (x->use_status) {
315: MPI_Get_count(&some_statuses[ii],intmsg ? MPIU_INT : MPIU_SCALAR,&count);
316: if (!intmsg) count /= bs; /* Convert from number of scalars to number of blocks */
317: } else count = x->recvhdr[i].bcount;
318: for (j=0,recvint=frame[i].intb,recvscalar=frame[i].scalarb; j<count; j++,recvint++) {
319: PetscInt loc = (*recvint)*bs - X->map->rstart;
320: switch (imode) {
321: case ADD_VALUES:
322: for (k=loc; k<loc+bs; k++) xarray[k] += *recvscalar++;
323: break;
324: case INSERT_VALUES:
325: for (k=loc; k<loc+bs; k++) xarray[k] = *recvscalar++;
326: break;
327: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Insert mode not supported 0x%x",imode);
328: }
329: }
330: }
331: }
332: }
333: VecRestoreArray(X,&xarray);
334: MPI_Waitall(4*x->nsendranks,x->sendreqs,MPI_STATUSES_IGNORE);
335: PetscFree2(some_indices,some_statuses);
336: if (x->assembly_subset) {
337: void *dummy; /* reset segbuffers */
338: PetscSegBufferExtractInPlace(x->segrecvint,&dummy);
339: PetscSegBufferExtractInPlace(x->segrecvscalar,&dummy);
340: } else {
341: VecAssemblyReset_MPI(X);
342: }
344: X->stash.insertmode = NOT_SET_VALUES;
345: X->bstash.insertmode = NOT_SET_VALUES;
346: VecStashScatterEnd_Private(&X->stash);
347: VecStashScatterEnd_Private(&X->bstash);
348: return 0;
349: }
351: PetscErrorCode VecAssemblyReset_MPI(Vec X)
352: {
353: Vec_MPI *x = (Vec_MPI*)X->data;
355: PetscFree(x->sendreqs);
356: PetscFree(x->recvreqs);
357: PetscFree(x->sendranks);
358: PetscFree(x->recvranks);
359: PetscFree(x->sendhdr);
360: PetscFree(x->recvhdr);
361: PetscFree(x->sendptrs);
362: PetscSegBufferDestroy(&x->segrecvint);
363: PetscSegBufferDestroy(&x->segrecvscalar);
364: PetscSegBufferDestroy(&x->segrecvframe);
365: return 0;
366: }
368: static PetscErrorCode VecSetFromOptions_MPI(PetscOptionItems *PetscOptionsObject,Vec X)
369: {
370: #if !defined(PETSC_HAVE_MPIUNI)
371: PetscBool flg = PETSC_FALSE,set;
373: PetscOptionsHead(PetscOptionsObject,"VecMPI Options");
374: PetscOptionsBool("-vec_assembly_legacy","Use MPI 1 version of assembly","",flg,&flg,&set);
375: if (set) {
376: X->ops->assemblybegin = flg ? VecAssemblyBegin_MPI : VecAssemblyBegin_MPI_BTS;
377: X->ops->assemblyend = flg ? VecAssemblyEnd_MPI : VecAssemblyEnd_MPI_BTS;
378: }
379: PetscOptionsTail();
380: #else
381: X->ops->assemblybegin = VecAssemblyBegin_MPI;
382: X->ops->assemblyend = VecAssemblyEnd_MPI;
383: #endif
384: return 0;
385: }
387: static struct _VecOps DvOps = {
388: PetscDesignatedInitializer(duplicate,VecDuplicate_MPI), /* 1 */
389: PetscDesignatedInitializer(duplicatevecs,VecDuplicateVecs_Default),
390: PetscDesignatedInitializer(destroyvecs,VecDestroyVecs_Default),
391: PetscDesignatedInitializer(dot,VecDot_MPI),
392: PetscDesignatedInitializer(mdot,VecMDot_MPI),
393: PetscDesignatedInitializer(norm,VecNorm_MPI),
394: PetscDesignatedInitializer(tdot,VecTDot_MPI),
395: PetscDesignatedInitializer(mtdot,VecMTDot_MPI),
396: PetscDesignatedInitializer(scale,VecScale_Seq),
397: PetscDesignatedInitializer(copy,VecCopy_Seq), /* 10 */
398: PetscDesignatedInitializer(set,VecSet_Seq),
399: PetscDesignatedInitializer(swap,VecSwap_Seq),
400: PetscDesignatedInitializer(axpy,VecAXPY_Seq),
401: PetscDesignatedInitializer(axpby,VecAXPBY_Seq),
402: PetscDesignatedInitializer(maxpy,VecMAXPY_Seq),
403: PetscDesignatedInitializer(aypx,VecAYPX_Seq),
404: PetscDesignatedInitializer(waxpy,VecWAXPY_Seq),
405: PetscDesignatedInitializer(axpbypcz,VecAXPBYPCZ_Seq),
406: PetscDesignatedInitializer(pointwisemult,VecPointwiseMult_Seq),
407: PetscDesignatedInitializer(pointwisedivide,VecPointwiseDivide_Seq),
408: PetscDesignatedInitializer(setvalues,VecSetValues_MPI), /* 20 */
409: PetscDesignatedInitializer(assemblybegin,VecAssemblyBegin_MPI_BTS),
410: PetscDesignatedInitializer(assemblyend,VecAssemblyEnd_MPI_BTS),
411: PetscDesignatedInitializer(getarray,NULL),
412: PetscDesignatedInitializer(getsize,VecGetSize_MPI),
413: PetscDesignatedInitializer(getlocalsize,VecGetSize_Seq),
414: PetscDesignatedInitializer(restorearray,NULL),
415: PetscDesignatedInitializer(max,VecMax_MPI),
416: PetscDesignatedInitializer(min,VecMin_MPI),
417: PetscDesignatedInitializer(setrandom,VecSetRandom_Seq),
418: PetscDesignatedInitializer(setoption,VecSetOption_MPI),
419: PetscDesignatedInitializer(setvaluesblocked,VecSetValuesBlocked_MPI),
420: PetscDesignatedInitializer(destroy,VecDestroy_MPI),
421: PetscDesignatedInitializer(view,VecView_MPI),
422: PetscDesignatedInitializer(placearray,VecPlaceArray_MPI),
423: PetscDesignatedInitializer(replacearray,VecReplaceArray_Seq),
424: PetscDesignatedInitializer(dot_local,VecDot_Seq),
425: PetscDesignatedInitializer(tdot_local,VecTDot_Seq),
426: PetscDesignatedInitializer(norm_local,VecNorm_Seq),
427: PetscDesignatedInitializer(mdot_local,VecMDot_Seq),
428: PetscDesignatedInitializer(mtdot_local,VecMTDot_Seq),
429: PetscDesignatedInitializer(load,VecLoad_Default),
430: PetscDesignatedInitializer(reciprocal,VecReciprocal_Default),
431: PetscDesignatedInitializer(conjugate,VecConjugate_Seq),
432: PetscDesignatedInitializer(setlocaltoglobalmapping,NULL),
433: PetscDesignatedInitializer(setvalueslocal,NULL),
434: PetscDesignatedInitializer(resetarray,VecResetArray_MPI),
435: PetscDesignatedInitializer(setfromoptions,VecSetFromOptions_MPI),/*set from options */
436: PetscDesignatedInitializer(maxpointwisedivide,VecMaxPointwiseDivide_Seq),
437: PetscDesignatedInitializer(pointwisemax,VecPointwiseMax_Seq),
438: PetscDesignatedInitializer(pointwisemaxabs,VecPointwiseMaxAbs_Seq),
439: PetscDesignatedInitializer(pointwisemin,VecPointwiseMin_Seq),
440: PetscDesignatedInitializer(getvalues,VecGetValues_MPI),
441: PetscDesignatedInitializer(sqrt,NULL),
442: PetscDesignatedInitializer(abs,NULL),
443: PetscDesignatedInitializer(exp,NULL),
444: PetscDesignatedInitializer(log,NULL),
445: PetscDesignatedInitializer(shift,NULL),
446: PetscDesignatedInitializer(create,NULL), /* really? */
447: PetscDesignatedInitializer(stridegather,VecStrideGather_Default),
448: PetscDesignatedInitializer(stridescatter,VecStrideScatter_Default),
449: PetscDesignatedInitializer(dotnorm2,NULL),
450: PetscDesignatedInitializer(getsubvector,NULL),
451: PetscDesignatedInitializer(restoresubvector,NULL),
452: PetscDesignatedInitializer(getarrayread,NULL),
453: PetscDesignatedInitializer(restorearrayread,NULL),
454: PetscDesignatedInitializer(stridesubsetgather,VecStrideSubSetGather_Default),
455: PetscDesignatedInitializer(stridesubsetscatter,VecStrideSubSetScatter_Default),
456: PetscDesignatedInitializer(viewnative,VecView_MPI),
457: PetscDesignatedInitializer(loadnative,NULL),
458: PetscDesignatedInitializer(getlocalvector,NULL),
459: };
461: /*
462: VecCreate_MPI_Private - Basic create routine called by VecCreate_MPI() (i.e. VecCreateMPI()),
463: VecCreateMPIWithArray(), VecCreate_Shared() (i.e. VecCreateShared()), VecCreateGhost(),
464: VecDuplicate_MPI(), VecCreateGhostWithArray(), VecDuplicate_MPI(), and VecDuplicate_Shared()
466: If alloc is true and array is NULL then this routine allocates the space, otherwise
467: no space is allocated.
468: */
469: PetscErrorCode VecCreate_MPI_Private(Vec v,PetscBool alloc,PetscInt nghost,const PetscScalar array[])
470: {
471: Vec_MPI *s;
473: PetscNewLog(v,&s);
474: v->data = (void*)s;
475: PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));
476: s->nghost = nghost;
477: v->petscnative = PETSC_TRUE;
478: if (array) v->offloadmask = PETSC_OFFLOAD_CPU;
480: PetscLayoutSetUp(v->map);
482: s->array = (PetscScalar*)array;
483: s->array_allocated = NULL;
484: if (alloc && !array) {
485: PetscInt n = v->map->n+nghost;
486: PetscCalloc1(n,&s->array);
487: PetscLogObjectMemory((PetscObject)v,n*sizeof(PetscScalar));
488: s->array_allocated = s->array;
489: }
491: /* By default parallel vectors do not have local representation */
492: s->localrep = NULL;
493: s->localupdate = NULL;
495: v->stash.insertmode = NOT_SET_VALUES;
496: v->bstash.insertmode = NOT_SET_VALUES;
497: /* create the stashes. The block-size for bstash is set later when
498: VecSetValuesBlocked is called.
499: */
500: VecStashCreate_Private(PetscObjectComm((PetscObject)v),1,&v->stash);
501: VecStashCreate_Private(PetscObjectComm((PetscObject)v),PetscAbs(v->map->bs),&v->bstash);
503: #if defined(PETSC_HAVE_MATLAB_ENGINE)
504: PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);
505: PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);
506: #endif
507: PetscObjectChangeTypeName((PetscObject)v,VECMPI);
508: return 0;
509: }
511: /*MC
512: VECMPI - VECMPI = "mpi" - The basic parallel vector
514: Options Database Keys:
515: . -vec_type mpi - sets the vector type to VECMPI during a call to VecSetFromOptions()
517: Level: beginner
519: .seealso: VecCreate(), VecSetType(), VecSetFromOptions(), VecCreateMPIWithArray(), VECMPI, VecType, VecCreateMPI(), VecCreateMPI()
520: M*/
522: PetscErrorCode VecCreate_MPI(Vec vv)
523: {
524: VecCreate_MPI_Private(vv,PETSC_TRUE,0,NULL);
525: return 0;
526: }
528: /*MC
529: VECSTANDARD = "standard" - A VECSEQ on one process and VECMPI on more than one process
531: Options Database Keys:
532: . -vec_type standard - sets a vector type to standard on calls to VecSetFromOptions()
534: Level: beginner
536: .seealso: VecCreateSeq(), VecCreateMPI()
537: M*/
539: PETSC_EXTERN PetscErrorCode VecCreate_Standard(Vec v)
540: {
541: PetscMPIInt size;
543: MPI_Comm_size(PetscObjectComm((PetscObject)v),&size);
544: if (size == 1) {
545: VecSetType(v,VECSEQ);
546: } else {
547: VecSetType(v,VECMPI);
548: }
549: return 0;
550: }
552: /*@C
553: VecCreateMPIWithArray - Creates a parallel, array-style vector,
554: where the user provides the array space to store the vector values.
556: Collective
558: Input Parameters:
559: + comm - the MPI communicator to use
560: . bs - block size, same meaning as VecSetBlockSize()
561: . n - local vector length, cannot be PETSC_DECIDE
562: . N - global vector length (or PETSC_DECIDE to have calculated)
563: - array - the user provided array to store the vector values
565: Output Parameter:
566: . vv - the vector
568: Notes:
569: Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
570: same type as an existing vector.
572: If the user-provided array is NULL, then VecPlaceArray() can be used
573: at a later stage to SET the array for storing the vector values.
575: PETSc does NOT free the array when the vector is destroyed via VecDestroy().
576: The user should not free the array until the vector is destroyed.
578: Level: intermediate
580: .seealso: VecCreateSeqWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateGhost(),
581: VecCreateMPI(), VecCreateGhostWithArray(), VecPlaceArray()
583: @*/
584: PetscErrorCode VecCreateMPIWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,const PetscScalar array[],Vec *vv)
585: {
587: PetscSplitOwnership(comm,&n,&N);
588: VecCreate(comm,vv);
589: VecSetSizes(*vv,n,N);
590: VecSetBlockSize(*vv,bs);
591: VecCreate_MPI_Private(*vv,PETSC_FALSE,0,array);
592: return 0;
593: }
595: /*@C
596: VecCreateGhostWithArray - Creates a parallel vector with ghost padding on each processor;
597: the caller allocates the array space.
599: Collective
601: Input Parameters:
602: + comm - the MPI communicator to use
603: . n - local vector length
604: . N - global vector length (or PETSC_DECIDE to have calculated if n is given)
605: . nghost - number of local ghost points
606: . ghosts - global indices of ghost points (or NULL if not needed), these do not need to be in increasing order (sorted)
607: - array - the space to store the vector values (as long as n + nghost)
609: Output Parameter:
610: . vv - the global vector representation (without ghost points as part of vector)
612: Notes:
613: Use VecGhostGetLocalForm() to access the local, ghosted representation
614: of the vector.
616: This also automatically sets the ISLocalToGlobalMapping() for this vector.
618: Level: advanced
620: .seealso: VecCreate(), VecGhostGetLocalForm(), VecGhostRestoreLocalForm(),
621: VecCreateGhost(), VecCreateSeqWithArray(), VecCreateMPIWithArray(),
622: VecCreateGhostBlock(), VecCreateGhostBlockWithArray(), VecMPISetGhost()
624: @*/
625: PetscErrorCode VecCreateGhostWithArray(MPI_Comm comm,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],const PetscScalar array[],Vec *vv)
626: {
627: Vec_MPI *w;
628: PetscScalar *larray;
629: IS from,to;
630: ISLocalToGlobalMapping ltog;
631: PetscInt rstart,i,*indices;
633: *vv = NULL;
638: PetscSplitOwnership(comm,&n,&N);
639: /* Create global representation */
640: VecCreate(comm,vv);
641: VecSetSizes(*vv,n,N);
642: VecCreate_MPI_Private(*vv,PETSC_TRUE,nghost,array);
643: w = (Vec_MPI*)(*vv)->data;
644: /* Create local representation */
645: VecGetArray(*vv,&larray);
646: VecCreateSeqWithArray(PETSC_COMM_SELF,1,n+nghost,larray,&w->localrep);
647: PetscLogObjectParent((PetscObject)*vv,(PetscObject)w->localrep);
648: VecRestoreArray(*vv,&larray);
650: /*
651: Create scatter context for scattering (updating) ghost values
652: */
653: ISCreateGeneral(comm,nghost,ghosts,PETSC_COPY_VALUES,&from);
654: ISCreateStride(PETSC_COMM_SELF,nghost,n,1,&to);
655: VecScatterCreate(*vv,from,w->localrep,to,&w->localupdate);
656: PetscLogObjectParent((PetscObject)*vv,(PetscObject)w->localupdate);
657: ISDestroy(&to);
658: ISDestroy(&from);
660: /* set local to global mapping for ghosted vector */
661: PetscMalloc1(n+nghost,&indices);
662: VecGetOwnershipRange(*vv,&rstart,NULL);
663: for (i=0; i<n; i++) {
664: indices[i] = rstart + i;
665: }
666: for (i=0; i<nghost; i++) {
667: indices[n+i] = ghosts[i];
668: }
669: ISLocalToGlobalMappingCreate(comm,1,n+nghost,indices,PETSC_OWN_POINTER,<og);
670: VecSetLocalToGlobalMapping(*vv,ltog);
671: ISLocalToGlobalMappingDestroy(<og);
672: return 0;
673: }
675: /*@
676: VecCreateGhost - Creates a parallel vector with ghost padding on each processor.
678: Collective
680: Input Parameters:
681: + comm - the MPI communicator to use
682: . n - local vector length
683: . N - global vector length (or PETSC_DECIDE to have calculated if n is given)
684: . nghost - number of local ghost points
685: - ghosts - global indices of ghost points, these do not need to be in increasing order (sorted)
687: Output Parameter:
688: . vv - the global vector representation (without ghost points as part of vector)
690: Notes:
691: Use VecGhostGetLocalForm() to access the local, ghosted representation
692: of the vector.
694: This also automatically sets the ISLocalToGlobalMapping() for this vector.
696: Level: advanced
698: .seealso: VecCreateSeq(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateMPI(),
699: VecGhostGetLocalForm(), VecGhostRestoreLocalForm(), VecGhostUpdateBegin(),
700: VecCreateGhostWithArray(), VecCreateMPIWithArray(), VecGhostUpdateEnd(),
701: VecCreateGhostBlock(), VecCreateGhostBlockWithArray(), VecMPISetGhost()
703: @*/
704: PetscErrorCode VecCreateGhost(MPI_Comm comm,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],Vec *vv)
705: {
706: VecCreateGhostWithArray(comm,n,N,nghost,ghosts,NULL,vv);
707: return 0;
708: }
710: /*@
711: VecMPISetGhost - Sets the ghost points for an MPI ghost vector
713: Collective on Vec
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 existant 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: (*vv->ops->destroy)(vv);
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: PetscLogObjectParent((PetscObject)vv,(PetscObject)w->localrep);
762: VecRestoreArray(vv,&larray);
764: /*
765: Create scatter context for scattering (updating) ghost values
766: */
767: ISCreateGeneral(comm,nghost,ghosts,PETSC_COPY_VALUES,&from);
768: ISCreateStride(PETSC_COMM_SELF,nghost,n,1,&to);
769: VecScatterCreate(vv,from,w->localrep,to,&w->localupdate);
770: PetscLogObjectParent((PetscObject)vv,(PetscObject)w->localupdate);
771: ISDestroy(&to);
772: ISDestroy(&from);
774: /* set local to global mapping for ghosted vector */
775: PetscMalloc1(n+nghost,&indices);
776: VecGetOwnershipRange(vv,&rstart,NULL);
778: for (i=0; i<n; i++) indices[i] = rstart + i;
779: for (i=0; i<nghost; i++) indices[n+i] = ghosts[i];
781: ISLocalToGlobalMappingCreate(comm,1,n+nghost,indices,PETSC_OWN_POINTER,<og);
782: VecSetLocalToGlobalMapping(vv,ltog);
783: ISLocalToGlobalMappingDestroy(<og);
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: PetscLogObjectParent((PetscObject)*vv,(PetscObject)w->localrep);
848: VecRestoreArray(*vv,&larray);
850: /*
851: Create scatter context for scattering (updating) ghost values
852: */
853: ISCreateBlock(comm,bs,nghost,ghosts,PETSC_COPY_VALUES,&from);
854: ISCreateStride(PETSC_COMM_SELF,bs*nghost,n,1,&to);
855: VecScatterCreate(*vv,from,w->localrep,to,&w->localupdate);
856: PetscLogObjectParent((PetscObject)*vv,(PetscObject)w->localupdate);
857: ISDestroy(&to);
858: ISDestroy(&from);
860: /* set local to global mapping for ghosted vector */
861: nb = n/bs;
862: PetscMalloc1(nb+nghost,&indices);
863: VecGetOwnershipRange(*vv,&rstart,NULL);
864: rstart = rstart/bs;
866: for (i=0; i<nb; i++) indices[i] = rstart + i;
867: for (i=0; i<nghost; i++) indices[nb+i] = ghosts[i];
869: ISLocalToGlobalMappingCreate(comm,bs,nb+nghost,indices,PETSC_OWN_POINTER,<og);
870: VecSetLocalToGlobalMapping(*vv,ltog);
871: ISLocalToGlobalMappingDestroy(<og);
872: return 0;
873: }
875: /*@
876: VecCreateGhostBlock - Creates a parallel vector with ghost padding on each processor.
877: The indicing of the ghost points is done with blocks.
879: Collective
881: Input Parameters:
882: + comm - the MPI communicator to use
883: . bs - the block size
884: . n - local vector length
885: . N - global vector length (or PETSC_DECIDE to have calculated if n is given)
886: . nghost - number of local ghost blocks
887: - ghosts - global indices of ghost blocks, counts are by block, not by individual index, these do not need to be in increasing order (sorted)
889: Output Parameter:
890: . vv - the global vector representation (without ghost points as part of vector)
892: Notes:
893: Use VecGhostGetLocalForm() to access the local, ghosted representation
894: of the vector.
896: n is the local vector size (total local size not the number of blocks) while nghost
897: is the number of blocks in the ghost portion, i.e. the number of elements in the ghost
898: portion is bs*nghost
900: Level: advanced
902: .seealso: VecCreateSeq(), VecCreate(), VecDuplicate(), VecDuplicateVecs(), VecCreateMPI(),
903: VecGhostGetLocalForm(), VecGhostRestoreLocalForm(),
904: VecCreateGhostWithArray(), VecCreateMPIWithArray(), VecCreateGhostBlockWithArray()
906: @*/
907: PetscErrorCode VecCreateGhostBlock(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,PetscInt nghost,const PetscInt ghosts[],Vec *vv)
908: {
909: VecCreateGhostBlockWithArray(comm,bs,n,N,nghost,ghosts,NULL,vv);
910: return 0;
911: }