Actual source code: vscat.c
petsc-3.6.4 2016-04-12
2: /*
3: Code for creating scatters between vectors. This file
4: includes the code for scattering between sequential vectors and
5: some special cases for parallel scatters.
6: */
8: #include <petsc/private/vecimpl.h> /*I "petscvec.h" I*/
10: #if defined(PETSC_HAVE_CUSP)
11: PETSC_INTERN PetscErrorCode VecScatterCUSPIndicesCreate_StoS(PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt*,PetscInt*,PetscCUSPIndices*);
12: PETSC_INTERN PetscErrorCode VecScatterCUSPIndicesDestroy(PetscCUSPIndices*);
13: PETSC_INTERN PetscErrorCode VecScatterCUSP_StoS(Vec,Vec,PetscCUSPIndices,InsertMode,ScatterMode);
14: #endif
16: /* Logging support */
17: PetscClassId VEC_SCATTER_CLASSID;
19: #if defined(PETSC_USE_DEBUG)
20: /*
21: Checks if any indices are less than zero and generates an error
22: */
25: static PetscErrorCode VecScatterCheckIndices_Private(PetscInt nmax,PetscInt n,const PetscInt *idx)
26: {
27: PetscInt i;
30: for (i=0; i<n; i++) {
31: if (idx[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative index %D at %D location",idx[i],i);
32: if (idx[i] >= nmax) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D at %D location greater than max %D",idx[i],i,nmax);
33: }
34: return(0);
35: }
36: #endif
38: /*
39: This is special scatter code for when the entire parallel vector is copied to each processor.
41: This code was written by Cameron Cooper, Occidental College, Fall 1995,
42: will working at ANL as a SERS student.
43: */
46: PetscErrorCode VecScatterBegin_MPI_ToAll(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
47: {
49: PetscInt yy_n,xx_n;
50: PetscScalar *xv,*yv;
53: VecGetLocalSize(y,&yy_n);
54: VecGetLocalSize(x,&xx_n);
55: VecGetArrayPair(x,y,&xv,&yv);
57: if (mode & SCATTER_REVERSE) {
58: PetscScalar *xvt,*xvt2;
59: VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
60: PetscInt i;
61: PetscMPIInt *disply = scat->displx;
63: if (addv == INSERT_VALUES) {
64: PetscInt rstart,rend;
65: /*
66: copy the correct part of the local vector into the local storage of
67: the MPI one. Note: this operation only makes sense if all the local
68: vectors have the same values
69: */
70: VecGetOwnershipRange(y,&rstart,&rend);
71: PetscMemcpy(yv,xv+rstart,yy_n*sizeof(PetscScalar));
72: } else {
73: MPI_Comm comm;
74: PetscMPIInt rank;
75: PetscObjectGetComm((PetscObject)y,&comm);
76: MPI_Comm_rank(comm,&rank);
77: if (scat->work1) xvt = scat->work1;
78: else {
79: PetscMalloc1(xx_n,&xvt);
80: scat->work1 = xvt;
81: }
82: if (!rank) { /* I am the zeroth processor, values are accumulated here */
83: if (scat->work2) xvt2 = scat->work2;
84: else {
85: PetscMalloc1(xx_n,&xvt2);
86: scat->work2 = xvt2;
87: }
88: MPI_Gatherv(yv,yy_n,MPIU_SCALAR,xvt2,scat->count,disply,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
89: MPI_Reduce(xv,xvt,xx_n,MPIU_SCALAR,MPIU_SUM,0,PetscObjectComm((PetscObject)ctx));
90: if (addv == ADD_VALUES) {
91: for (i=0; i<xx_n; i++) xvt[i] += xvt2[i];
92: #if !defined(PETSC_USE_COMPLEX)
93: } else if (addv == MAX_VALUES) {
94: for (i=0; i<xx_n; i++) xvt[i] = PetscMax(xvt[i],xvt2[i]);
95: #endif
96: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
97: MPI_Scatterv(xvt,scat->count,disply,MPIU_SCALAR,yv,yy_n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
98: } else {
99: MPI_Gatherv(yv,yy_n,MPIU_SCALAR,0, 0,0,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
100: MPI_Reduce(xv,xvt,xx_n,MPIU_SCALAR,MPIU_SUM,0,PetscObjectComm((PetscObject)ctx));
101: MPI_Scatterv(0,scat->count,disply,MPIU_SCALAR,yv,yy_n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
102: }
103: }
104: } else {
105: PetscScalar *yvt;
106: VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
107: PetscInt i;
108: PetscMPIInt *displx = scat->displx;
110: if (addv == INSERT_VALUES) {
111: MPI_Allgatherv(xv,xx_n,MPIU_SCALAR,yv,scat->count,displx,MPIU_SCALAR,PetscObjectComm((PetscObject)ctx));
112: } else {
113: if (scat->work1) yvt = scat->work1;
114: else {
115: PetscMalloc1(yy_n,&yvt);
116: scat->work1 = yvt;
117: }
118: MPI_Allgatherv(xv,xx_n,MPIU_SCALAR,yvt,scat->count,displx,MPIU_SCALAR,PetscObjectComm((PetscObject)ctx));
119: if (addv == ADD_VALUES) {
120: for (i=0; i<yy_n; i++) yv[i] += yvt[i];
121: #if !defined(PETSC_USE_COMPLEX)
122: } else if (addv == MAX_VALUES) {
123: for (i=0; i<yy_n; i++) yv[i] = PetscMax(yv[i],yvt[i]);
124: #endif
125: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
126: }
127: }
128: VecRestoreArrayPair(x,y,&xv,&yv);
129: return(0);
130: }
134: PetscErrorCode VecScatterView_MPI_ToAll(VecScatter in,PetscViewer viewer)
135: {
137: PetscBool isascii;
140: PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
141: if (isascii) {
142: PetscViewerASCIIPrintf(viewer,"Entire parallel vector is copied to each process\n");
143: }
144: return(0);
145: }
147: /*
148: This is special scatter code for when the entire parallel vector is copied to processor 0.
150: */
153: PetscErrorCode VecScatterBegin_MPI_ToOne(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
154: {
155: PetscErrorCode ierr;
156: PetscMPIInt rank;
157: PetscInt yy_n,xx_n;
158: PetscScalar *yv;
159: const PetscScalar *xv;
160: MPI_Comm comm;
163: VecGetLocalSize(y,&yy_n);
164: VecGetLocalSize(x,&xx_n);
165: VecGetArrayRead(x,&xv);
166: VecGetArray(y,&yv);
168: PetscObjectGetComm((PetscObject)x,&comm);
169: MPI_Comm_rank(comm,&rank);
171: /* -------- Reverse scatter; spread from processor 0 to other processors */
172: if (mode & SCATTER_REVERSE) {
173: PetscScalar *yvt;
174: VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
175: PetscInt i;
176: PetscMPIInt *disply = scat->displx;
178: if (addv == INSERT_VALUES) {
179: MPI_Scatterv((PetscScalar*)xv,scat->count,disply,MPIU_SCALAR,yv,yy_n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
180: } else {
181: if (scat->work2) yvt = scat->work2;
182: else {
183: PetscMalloc1(xx_n,&yvt);
184: scat->work2 = yvt;
185: }
186: MPI_Scatterv((PetscScalar*)xv,scat->count,disply,MPIU_SCALAR,yvt,yy_n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
187: if (addv == ADD_VALUES) {
188: for (i=0; i<yy_n; i++) yv[i] += yvt[i];
189: #if !defined(PETSC_USE_COMPLEX)
190: } else if (addv == MAX_VALUES) {
191: for (i=0; i<yy_n; i++) yv[i] = PetscMax(yv[i],yvt[i]);
192: #endif
193: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
194: }
195: /* --------- Forward scatter; gather all values onto processor 0 */
196: } else {
197: PetscScalar *yvt = 0;
198: VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
199: PetscInt i;
200: PetscMPIInt *displx = scat->displx;
202: if (addv == INSERT_VALUES) {
203: MPI_Gatherv((PetscScalar*)xv,xx_n,MPIU_SCALAR,yv,scat->count,displx,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
204: } else {
205: if (!rank) {
206: if (scat->work1) yvt = scat->work1;
207: else {
208: PetscMalloc1(yy_n,&yvt);
209: scat->work1 = yvt;
210: }
211: }
212: MPI_Gatherv((PetscScalar*)xv,xx_n,MPIU_SCALAR,yvt,scat->count,displx,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
213: if (!rank) {
214: if (addv == ADD_VALUES) {
215: for (i=0; i<yy_n; i++) yv[i] += yvt[i];
216: #if !defined(PETSC_USE_COMPLEX)
217: } else if (addv == MAX_VALUES) {
218: for (i=0; i<yy_n; i++) yv[i] = PetscMax(yv[i],yvt[i]);
219: #endif
220: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
221: }
222: }
223: }
224: VecRestoreArrayRead(x,&xv);
225: VecRestoreArray(y,&yv);
226: return(0);
227: }
229: /*
230: The follow to are used for both VecScatterBegin_MPI_ToAll() and VecScatterBegin_MPI_ToOne()
231: */
234: PetscErrorCode VecScatterDestroy_MPI_ToAll(VecScatter ctx)
235: {
236: VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
237: PetscErrorCode ierr;
240: PetscFree(scat->work1);
241: PetscFree(scat->work2);
242: PetscFree3(ctx->todata,scat->count,scat->displx);
243: return(0);
244: }
248: PetscErrorCode VecScatterDestroy_SGToSG(VecScatter ctx)
249: {
253: PetscFree2(((VecScatter_Seq_General*)ctx->todata)->vslots,((VecScatter_Seq_General*)ctx->fromdata)->vslots);
254: PetscFree2(ctx->todata,ctx->fromdata);
255: return(0);
256: }
260: PetscErrorCode VecScatterDestroy_SGToSS(VecScatter ctx)
261: {
265: PetscFree(((VecScatter_Seq_General*)ctx->fromdata)->vslots);
266: PetscFree2(ctx->todata,ctx->fromdata);
267: return(0);
268: }
272: PetscErrorCode VecScatterDestroy_SSToSG(VecScatter ctx)
273: {
277: PetscFree(((VecScatter_Seq_General*)ctx->todata)->vslots);
278: PetscFree2(ctx->todata,ctx->fromdata);
279: return(0);
280: }
284: PetscErrorCode VecScatterDestroy_SSToSS(VecScatter ctx)
285: {
289: PetscFree2(ctx->todata,ctx->fromdata);
290: return(0);
291: }
293: /* -------------------------------------------------------------------------*/
296: PetscErrorCode VecScatterCopy_MPI_ToAll(VecScatter in,VecScatter out)
297: {
298: VecScatter_MPI_ToAll *in_to = (VecScatter_MPI_ToAll*)in->todata,*sto;
299: PetscErrorCode ierr;
300: PetscMPIInt size,*count,*displx;
301: PetscInt i;
304: out->ops->begin = in->ops->begin;
305: out->ops->end = in->ops->end;
306: out->ops->copy = in->ops->copy;
307: out->ops->destroy = in->ops->destroy;
308: out->ops->view = in->ops->view;
310: MPI_Comm_size(PetscObjectComm((PetscObject)out),&size);
311: PetscMalloc3(1,&sto,size,&count,size,&displx);
312: sto->type = in_to->type;
313: sto->count = count;
314: sto->displx = displx;
315: for (i=0; i<size; i++) {
316: sto->count[i] = in_to->count[i];
317: sto->displx[i] = in_to->displx[i];
318: }
319: sto->work1 = 0;
320: sto->work2 = 0;
321: out->todata = (void*)sto;
322: out->fromdata = (void*)0;
323: return(0);
324: }
326: /* --------------------------------------------------------------------------------------*/
327: /*
328: Scatter: sequential general to sequential general
329: */
332: PetscErrorCode VecScatterBegin_SGToSG(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
333: {
334: VecScatter_Seq_General *gen_to = (VecScatter_Seq_General*)ctx->todata;
335: VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
336: PetscErrorCode ierr;
337: PetscInt i,n = gen_from->n,*fslots,*tslots;
338: PetscScalar *xv,*yv;
341: #if defined(PETSC_HAVE_CUSP)
342: if (x->valid_GPU_array == PETSC_CUSP_GPU) {
343: /* create the scatter indices if not done already */
344: if (!ctx->spptr) {
345: PetscInt tofirst = 0,tostep = 0,fromfirst = 0,fromstep = 0;
346: fslots = gen_from->vslots;
347: tslots = gen_to->vslots;
348: VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
349: }
350: /* next do the scatter */
351: VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
352: return(0);
353: }
354: #endif
356: VecGetArrayPair(x,y,&xv,&yv);
357: if (mode & SCATTER_REVERSE) {
358: gen_to = (VecScatter_Seq_General*)ctx->fromdata;
359: gen_from = (VecScatter_Seq_General*)ctx->todata;
360: }
361: fslots = gen_from->vslots;
362: tslots = gen_to->vslots;
364: if (addv == INSERT_VALUES) {
365: for (i=0; i<n; i++) yv[tslots[i]] = xv[fslots[i]];
366: } else if (addv == ADD_VALUES) {
367: for (i=0; i<n; i++) yv[tslots[i]] += xv[fslots[i]];
368: #if !defined(PETSC_USE_COMPLEX)
369: } else if (addv == MAX_VALUES) {
370: for (i=0; i<n; i++) yv[tslots[i]] = PetscMax(yv[tslots[i]],xv[fslots[i]]);
371: #endif
372: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
373: VecRestoreArrayPair(x,y,&xv,&yv);
374: return(0);
375: }
377: /*
378: Scatter: sequential general to sequential stride 1
379: */
382: PetscErrorCode VecScatterBegin_SGToSS_Stride1(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
383: {
384: VecScatter_Seq_Stride *gen_to = (VecScatter_Seq_Stride*)ctx->todata;
385: VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
386: PetscInt i,n = gen_from->n,*fslots = gen_from->vslots;
387: PetscErrorCode ierr;
388: PetscInt first = gen_to->first;
389: PetscScalar *xv,*yv;
392: #if defined(PETSC_HAVE_CUSP)
393: if (x->valid_GPU_array == PETSC_CUSP_GPU) {
394: /* create the scatter indices if not done already */
395: if (!ctx->spptr) {
396: PetscInt tofirst = first,tostep = 1,fromfirst = 0,fromstep = 0;
397: PetscInt *tslots = 0;
398: VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
399: }
400: /* next do the scatter */
401: VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
402: return(0);
403: }
404: #endif
406: VecGetArrayPair(x,y,&xv,&yv);
407: if (mode & SCATTER_REVERSE) {
408: xv += first;
409: if (addv == INSERT_VALUES) {
410: for (i=0; i<n; i++) yv[fslots[i]] = xv[i];
411: } else if (addv == ADD_VALUES) {
412: for (i=0; i<n; i++) yv[fslots[i]] += xv[i];
413: #if !defined(PETSC_USE_COMPLEX)
414: } else if (addv == MAX_VALUES) {
415: for (i=0; i<n; i++) yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[i]);
416: #endif
417: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
418: } else {
419: yv += first;
420: if (addv == INSERT_VALUES) {
421: for (i=0; i<n; i++) yv[i] = xv[fslots[i]];
422: } else if (addv == ADD_VALUES) {
423: for (i=0; i<n; i++) yv[i] += xv[fslots[i]];
424: #if !defined(PETSC_USE_COMPLEX)
425: } else if (addv == MAX_VALUES) {
426: for (i=0; i<n; i++) yv[i] = PetscMax(yv[i],xv[fslots[i]]);
427: #endif
428: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
429: }
430: VecRestoreArrayPair(x,y,&xv,&yv);
431: return(0);
432: }
434: /*
435: Scatter: sequential general to sequential stride
436: */
439: PetscErrorCode VecScatterBegin_SGToSS(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
440: {
441: VecScatter_Seq_Stride *gen_to = (VecScatter_Seq_Stride*)ctx->todata;
442: VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
443: PetscInt i,n = gen_from->n,*fslots = gen_from->vslots;
444: PetscErrorCode ierr;
445: PetscInt first = gen_to->first,step = gen_to->step;
446: PetscScalar *xv,*yv;
449: #if defined(PETSC_HAVE_CUSP)
450: if (x->valid_GPU_array == PETSC_CUSP_GPU) {
451: /* create the scatter indices if not done already */
452: if (!ctx->spptr) {
453: PetscInt tofirst = first,tostep = step,fromfirst = 0,fromstep = 0;
454: PetscInt * tslots = 0;
455: VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
456: }
457: /* next do the scatter */
458: VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
459: return(0);
460: }
461: #endif
463: VecGetArrayPair(x,y,&xv,&yv);
464: if (mode & SCATTER_REVERSE) {
465: if (addv == INSERT_VALUES) {
466: for (i=0; i<n; i++) yv[fslots[i]] = xv[first + i*step];
467: } else if (addv == ADD_VALUES) {
468: for (i=0; i<n; i++) yv[fslots[i]] += xv[first + i*step];
469: #if !defined(PETSC_USE_COMPLEX)
470: } else if (addv == MAX_VALUES) {
471: for (i=0; i<n; i++) yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[first + i*step]);
472: #endif
473: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
474: } else {
475: if (addv == INSERT_VALUES) {
476: for (i=0; i<n; i++) yv[first + i*step] = xv[fslots[i]];
477: } else if (addv == ADD_VALUES) {
478: for (i=0; i<n; i++) yv[first + i*step] += xv[fslots[i]];
479: #if !defined(PETSC_USE_COMPLEX)
480: } else if (addv == MAX_VALUES) {
481: for (i=0; i<n; i++) yv[first + i*step] = PetscMax(yv[first + i*step],xv[fslots[i]]);
482: #endif
483: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
484: }
485: VecRestoreArrayPair(x,y,&xv,&yv);
486: return(0);
487: }
489: /*
490: Scatter: sequential stride 1 to sequential general
491: */
494: PetscErrorCode VecScatterBegin_SSToSG_Stride1(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
495: {
496: VecScatter_Seq_Stride *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
497: VecScatter_Seq_General *gen_to = (VecScatter_Seq_General*)ctx->todata;
498: PetscInt i,n = gen_from->n,*tslots = gen_to->vslots;
499: PetscErrorCode ierr;
500: PetscInt first = gen_from->first;
501: PetscScalar *xv,*yv;
504: #if defined(PETSC_HAVE_CUSP)
505: if (x->valid_GPU_array == PETSC_CUSP_GPU) {
506: /* create the scatter indices if not done already */
507: if (!ctx->spptr) {
508: PetscInt tofirst = 0,tostep = 0,fromfirst = first,fromstep = 1;
509: PetscInt *fslots = 0;
510: VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
511: }
512: /* next do the scatter */
513: VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
514: return(0);
515: }
516: #endif
518: VecGetArrayPair(x,y,&xv,&yv);
519: if (mode & SCATTER_REVERSE) {
520: yv += first;
521: if (addv == INSERT_VALUES) {
522: for (i=0; i<n; i++) yv[i] = xv[tslots[i]];
523: } else if (addv == ADD_VALUES) {
524: for (i=0; i<n; i++) yv[i] += xv[tslots[i]];
525: #if !defined(PETSC_USE_COMPLEX)
526: } else if (addv == MAX_VALUES) {
527: for (i=0; i<n; i++) yv[i] = PetscMax(yv[i],xv[tslots[i]]);
528: #endif
529: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
530: } else {
531: xv += first;
532: if (addv == INSERT_VALUES) {
533: for (i=0; i<n; i++) yv[tslots[i]] = xv[i];
534: } else if (addv == ADD_VALUES) {
535: for (i=0; i<n; i++) yv[tslots[i]] += xv[i];
536: #if !defined(PETSC_USE_COMPLEX)
537: } else if (addv == MAX_VALUES) {
538: for (i=0; i<n; i++) yv[tslots[i]] = PetscMax(yv[tslots[i]],xv[i]);
539: #endif
540: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
541: }
542: VecRestoreArrayPair(x,y,&xv,&yv);
543: return(0);
544: }
548: /*
549: Scatter: sequential stride to sequential general
550: */
551: PetscErrorCode VecScatterBegin_SSToSG(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
552: {
553: VecScatter_Seq_Stride *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
554: VecScatter_Seq_General *gen_to = (VecScatter_Seq_General*)ctx->todata;
555: PetscInt i,n = gen_from->n,*tslots = gen_to->vslots;
556: PetscErrorCode ierr;
557: PetscInt first = gen_from->first,step = gen_from->step;
558: PetscScalar *xv,*yv;
561: #if defined(PETSC_HAVE_CUSP)
562: if (x->valid_GPU_array == PETSC_CUSP_GPU) {
563: /* create the scatter indices if not done already */
564: if (!ctx->spptr) {
565: PetscInt tofirst = 0,tostep = 0,fromfirst = first,fromstep = step;
566: PetscInt *fslots = 0;
567: VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
568: }
569: /* next do the scatter */
570: VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
571: return(0);
572: }
573: #endif
575: VecGetArrayPair(x,y,&xv,&yv);
576: if (mode & SCATTER_REVERSE) {
577: if (addv == INSERT_VALUES) {
578: for (i=0; i<n; i++) yv[first + i*step] = xv[tslots[i]];
579: } else if (addv == ADD_VALUES) {
580: for (i=0; i<n; i++) yv[first + i*step] += xv[tslots[i]];
581: #if !defined(PETSC_USE_COMPLEX)
582: } else if (addv == MAX_VALUES) {
583: for (i=0; i<n; i++) yv[first + i*step] = PetscMax(yv[first + i*step],xv[tslots[i]]);
584: #endif
585: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
586: } else {
587: if (addv == INSERT_VALUES) {
588: for (i=0; i<n; i++) yv[tslots[i]] = xv[first + i*step];
589: } else if (addv == ADD_VALUES) {
590: for (i=0; i<n; i++) yv[tslots[i]] += xv[first + i*step];
591: #if !defined(PETSC_USE_COMPLEX)
592: } else if (addv == MAX_VALUES) {
593: for (i=0; i<n; i++) yv[tslots[i]] = PetscMax(yv[tslots[i]],xv[first + i*step]);
594: #endif
595: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
596: }
597: VecRestoreArrayPair(x,y,&xv,&yv);
598: return(0);
599: }
603: PetscErrorCode VecScatterView_SSToSG(VecScatter in,PetscViewer viewer)
604: {
605: PetscErrorCode ierr;
606: VecScatter_Seq_Stride *in_from = (VecScatter_Seq_Stride*)in->todata;
607: VecScatter_Seq_General *in_to = (VecScatter_Seq_General*)in->fromdata;
608: PetscInt i;
609: PetscBool isascii;
612: PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
613: if (isascii) {
614: PetscViewerASCIIPrintf(viewer,"Sequential stride to general scatter\n");
615: for (i=0; i<in_to->n; i++) {
616: PetscViewerASCIIPrintf(viewer,"%D to %D\n",in_from->first + in_from->step*i,in_to->vslots[i]);
617: }
618: }
619: return(0);
620: }
621: /*
622: Scatter: sequential stride to sequential stride
623: */
626: PetscErrorCode VecScatterBegin_SSToSS(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
627: {
628: VecScatter_Seq_Stride *gen_to = (VecScatter_Seq_Stride*)ctx->todata;
629: VecScatter_Seq_Stride *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
630: PetscInt i,n = gen_from->n,to_first = gen_to->first,to_step = gen_to->step;
631: PetscErrorCode ierr;
632: PetscInt from_first = gen_from->first,from_step = gen_from->step;
633: PetscScalar *xv,*yv;
636: #if defined(PETSC_HAVE_CUSP)
637: if (x->valid_GPU_array == PETSC_CUSP_GPU) {
638: /* create the scatter indices if not done already */
639: if (!ctx->spptr) {
640: PetscInt *tslots = 0,*fslots = 0;
641: VecScatterCUSPIndicesCreate_StoS(n,to_first,from_first,to_step,from_step,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
642: }
643: /* next do the scatter */
644: VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
645: return(0);
646: }
647: #endif
649: VecGetArrayPair(x,y,&xv,&yv);
650: if (mode & SCATTER_REVERSE) {
651: from_first = gen_to->first;
652: to_first = gen_from->first;
653: from_step = gen_to->step;
654: to_step = gen_from->step;
655: }
657: if (addv == INSERT_VALUES) {
658: if (to_step == 1 && from_step == 1) {
659: PetscMemcpy(yv+to_first,xv+from_first,n*sizeof(PetscScalar));
660: } else {
661: for (i=0; i<n; i++) yv[to_first + i*to_step] = xv[from_first+i*from_step];
662: }
663: } else if (addv == ADD_VALUES) {
664: if (to_step == 1 && from_step == 1) {
665: yv += to_first; xv += from_first;
666: for (i=0; i<n; i++) yv[i] += xv[i];
667: } else {
668: for (i=0; i<n; i++) yv[to_first + i*to_step] += xv[from_first+i*from_step];
669: }
670: #if !defined(PETSC_USE_COMPLEX)
671: } else if (addv == MAX_VALUES) {
672: if (to_step == 1 && from_step == 1) {
673: yv += to_first; xv += from_first;
674: for (i=0; i<n; i++) yv[i] = PetscMax(yv[i],xv[i]);
675: } else {
676: for (i=0; i<n; i++) yv[to_first + i*to_step] = PetscMax(yv[to_first + i*to_step],xv[from_first+i*from_step]);
677: }
678: #endif
679: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
680: VecRestoreArrayPair(x,y,&xv,&yv);
681: return(0);
682: }
684: /* --------------------------------------------------------------------------*/
689: PetscErrorCode VecScatterCopy_SGToSG(VecScatter in,VecScatter out)
690: {
691: PetscErrorCode ierr;
692: VecScatter_Seq_General *in_to = (VecScatter_Seq_General*)in->todata,*out_to = NULL;
693: VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from = NULL;
696: out->ops->begin = in->ops->begin;
697: out->ops->end = in->ops->end;
698: out->ops->copy = in->ops->copy;
699: out->ops->destroy = in->ops->destroy;
700: out->ops->view = in->ops->view;
702: PetscMalloc2(1,&out_to,1,&out_from);
703: PetscMalloc2(in_to->n,&out_to->vslots,in_from->n,&out_from->vslots);
704: out_to->n = in_to->n;
705: out_to->type = in_to->type;
706: out_to->nonmatching_computed = PETSC_FALSE;
707: out_to->slots_nonmatching = 0;
708: out_to->is_copy = PETSC_FALSE;
709: PetscMemcpy(out_to->vslots,in_to->vslots,(out_to->n)*sizeof(PetscInt));
711: out_from->n = in_from->n;
712: out_from->type = in_from->type;
713: out_from->nonmatching_computed = PETSC_FALSE;
714: out_from->slots_nonmatching = 0;
715: out_from->is_copy = PETSC_FALSE;
716: PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));
718: out->todata = (void*)out_to;
719: out->fromdata = (void*)out_from;
720: return(0);
721: }
725: PetscErrorCode VecScatterView_SGToSG(VecScatter in,PetscViewer viewer)
726: {
727: PetscErrorCode ierr;
728: VecScatter_Seq_General *in_to = (VecScatter_Seq_General*)in->todata;
729: VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata;
730: PetscInt i;
731: PetscBool isascii;
734: PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
735: if (isascii) {
736: PetscViewerASCIIPrintf(viewer,"Sequential general scatter\n");
737: for (i=0; i<in_to->n; i++) {
738: PetscViewerASCIIPrintf(viewer,"%D to %D\n",in_from->vslots[i],in_to->vslots[i]);
739: }
740: }
741: return(0);
742: }
747: PetscErrorCode VecScatterCopy_SGToSS(VecScatter in,VecScatter out)
748: {
749: PetscErrorCode ierr;
750: VecScatter_Seq_Stride *in_to = (VecScatter_Seq_Stride*)in->todata,*out_to = NULL;
751: VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from = NULL;
754: out->ops->begin = in->ops->begin;
755: out->ops->end = in->ops->end;
756: out->ops->copy = in->ops->copy;
757: out->ops->destroy = in->ops->destroy;
758: out->ops->view = in->ops->view;
760: PetscMalloc2(1,&out_to,1,&out_from);
761: PetscMalloc1(in_from->n,&out_from->vslots);
762: out_to->n = in_to->n;
763: out_to->type = in_to->type;
764: out_to->first = in_to->first;
765: out_to->step = in_to->step;
766: out_to->type = in_to->type;
768: out_from->n = in_from->n;
769: out_from->type = in_from->type;
770: out_from->nonmatching_computed = PETSC_FALSE;
771: out_from->slots_nonmatching = 0;
772: out_from->is_copy = PETSC_FALSE;
773: PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));
775: out->todata = (void*)out_to;
776: out->fromdata = (void*)out_from;
777: return(0);
778: }
782: PetscErrorCode VecScatterView_SGToSS(VecScatter in,PetscViewer viewer)
783: {
784: PetscErrorCode ierr;
785: VecScatter_Seq_Stride *in_to = (VecScatter_Seq_Stride*)in->todata;
786: VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata;
787: PetscInt i;
788: PetscBool isascii;
791: PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
792: if (isascii) {
793: PetscViewerASCIIPrintf(viewer,"Sequential general scatter to stride\n");
794: for (i=0; i<in_to->n; i++) {
795: PetscViewerASCIIPrintf(viewer,"%D to %D\n",in_from->vslots[i],in_to->first + in_to->step*i);
796: }
797: }
798: return(0);
799: }
801: /* --------------------------------------------------------------------------*/
802: /*
803: Scatter: parallel to sequential vector, sequential strides for both.
804: */
807: PetscErrorCode VecScatterCopy_SSToSS(VecScatter in,VecScatter out)
808: {
809: VecScatter_Seq_Stride *in_to = (VecScatter_Seq_Stride*)in->todata,*out_to = NULL;
810: VecScatter_Seq_Stride *in_from = (VecScatter_Seq_Stride*)in->fromdata,*out_from = NULL;
811: PetscErrorCode ierr;
814: out->ops->begin = in->ops->begin;
815: out->ops->end = in->ops->end;
816: out->ops->copy = in->ops->copy;
817: out->ops->destroy = in->ops->destroy;
818: out->ops->view = in->ops->view;
820: PetscMalloc2(1,&out_to,1,&out_from);
821: out_to->n = in_to->n;
822: out_to->type = in_to->type;
823: out_to->first = in_to->first;
824: out_to->step = in_to->step;
825: out_to->type = in_to->type;
826: out_from->n = in_from->n;
827: out_from->type = in_from->type;
828: out_from->first = in_from->first;
829: out_from->step = in_from->step;
830: out_from->type = in_from->type;
831: out->todata = (void*)out_to;
832: out->fromdata = (void*)out_from;
833: return(0);
834: }
838: PetscErrorCode VecScatterView_SSToSS(VecScatter in,PetscViewer viewer)
839: {
840: VecScatter_Seq_Stride *in_to = (VecScatter_Seq_Stride*)in->todata;
841: VecScatter_Seq_Stride *in_from = (VecScatter_Seq_Stride*)in->fromdata;
842: PetscErrorCode ierr;
843: PetscBool isascii;
846: PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
847: if (isascii) {
848: PetscViewerASCIIPrintf(viewer,"Sequential stride count %D start %D step to start %D stride %D\n",in_to->n,in_to->first,in_to->step,in_from->first,in_from->step);
849: }
850: return(0);
851: }
854: extern PetscErrorCode VecScatterCreate_PtoS(PetscInt,const PetscInt*,PetscInt,const PetscInt*,Vec,Vec,PetscInt,VecScatter);
855: extern PetscErrorCode VecScatterCreate_PtoP(PetscInt,const PetscInt*,PetscInt,const PetscInt*,Vec,Vec,PetscInt,VecScatter);
856: extern PetscErrorCode VecScatterCreate_StoP(PetscInt,const PetscInt*,PetscInt,const PetscInt*,Vec,Vec,PetscInt,VecScatter);
858: /* =======================================================================*/
859: #define VEC_SEQ_ID 0
860: #define VEC_MPI_ID 1
861: #define IS_GENERAL_ID 0
862: #define IS_STRIDE_ID 1
863: #define IS_BLOCK_ID 2
865: /*
866: Blocksizes we have optimized scatters for
867: */
869: #define VecScatterOptimizedBS(mbs) (2 <= mbs)
871: PetscErrorCode VecScatterCreateEmpty(MPI_Comm comm,VecScatter *newctx)
872: {
873: VecScatter ctx;
877: PetscHeaderCreate(ctx,VEC_SCATTER_CLASSID,"VecScatter","VecScatter","Vec",comm,VecScatterDestroy,VecScatterView);
878: ctx->inuse = PETSC_FALSE;
879: ctx->beginandendtogether = PETSC_FALSE;
880: PetscOptionsGetBool(NULL,"-vecscatter_merge",&ctx->beginandendtogether,NULL);
881: if (ctx->beginandendtogether) {
882: PetscInfo(ctx,"Using combined (merged) vector scatter begin and end\n");
883: }
884: ctx->packtogether = PETSC_FALSE;
885: PetscOptionsGetBool(NULL,"-vecscatter_packtogether",&ctx->packtogether,NULL);
886: if (ctx->packtogether) {
887: PetscInfo(ctx,"Pack all messages before sending\n");
888: }
889: *newctx = ctx;
890: return(0);
891: }
895: /*@C
896: VecScatterCreate - Creates a vector scatter context.
898: Collective on Vec
900: Input Parameters:
901: + xin - a vector that defines the shape (parallel data layout of the vector)
902: of vectors from which we scatter
903: . yin - a vector that defines the shape (parallel data layout of the vector)
904: of vectors to which we scatter
905: . ix - the indices of xin to scatter (if NULL scatters all values)
906: - iy - the indices of yin to hold results (if NULL fills entire vector yin)
908: Output Parameter:
909: . newctx - location to store the new scatter context
911: Options Database Keys: (uses regular MPI_Sends by default)
912: + -vecscatter_view - Prints detail of communications
913: . -vecscatter_view ::ascii_info - Print less details about communication
914: . -vecscatter_ssend - Uses MPI_Ssend_init() instead of MPI_Send_init()
915: . -vecscatter_rsend - use ready receiver mode for MPI sends
916: . -vecscatter_merge - VecScatterBegin() handles all of the communication, VecScatterEnd() is a nop
917: eliminates the chance for overlap of computation and communication
918: . -vecscatter_sendfirst - Posts sends before receives
919: . -vecscatter_packtogether - Pack all messages before sending, receive all messages before unpacking
920: . -vecscatter_alltoall - Uses MPI all to all communication for scatter
921: . -vecscatter_window - Use MPI 2 window operations to move data
922: . -vecscatter_nopack - Avoid packing to work vector when possible (if used with -vecscatter_alltoall then will use MPI_Alltoallw()
923: - -vecscatter_reproduce - insure that the order of the communications are done the same for each scatter, this under certain circumstances
924: will make the results of scatters deterministic when otherwise they are not (it may be slower also).
926: $
927: $ --When packing is used--
928: $ MPI Datatypes (no packing) sendfirst merge packtogether persistent*
929: $ _nopack _sendfirst _merge _packtogether -vecscatter_
930: $ ----------------------------------------------------------------------------------------------------------------------------
931: $ Message passing Send p X X X always
932: $ Ssend p X X X always _ssend
933: $ Rsend p nonsense X X always _rsend
934: $ AlltoAll v or w X nonsense always X nonsense _alltoall
935: $ MPI_Win p nonsense p p nonsense _window
936: $
937: $ Since persistent sends and receives require a constant memory address they can only be used when data is packed into the work vector
938: $ because the in and out array may be different for each call to VecScatterBegin/End().
939: $
940: $ p indicates possible, but not implemented. X indicates implemented
941: $
943: Level: intermediate
945: Notes:
946: In calls to VecScatter() you can use different vectors than the xin and
947: yin you used above; BUT they must have the same parallel data layout, for example,
948: they could be obtained from VecDuplicate().
949: A VecScatter context CANNOT be used in two or more simultaneous scatters;
950: that is you cannot call a second VecScatterBegin() with the same scatter
951: context until the VecScatterEnd() has been called on the first VecScatterBegin().
952: In this case a separate VecScatter is needed for each concurrent scatter.
954: Currently the MPI_Send(), MPI_Ssend() and MPI_Rsend() all use PERSISTENT versions.
955: (this unfortunately requires that the same in and out arrays be used for each use, this
956: is why when not using MPI_alltoallw() we always need to pack the input into the work array before sending
957: and unpack upon receeving instead of using MPI datatypes to avoid the packing/unpacking).
959: Both ix and iy cannot be NULL at the same time.
961: Concepts: scatter^between vectors
962: Concepts: gather^between vectors
964: .seealso: VecScatterDestroy(), VecScatterCreateToAll(), VecScatterCreateToZero()
965: @*/
966: PetscErrorCode VecScatterCreate(Vec xin,IS ix,Vec yin,IS iy,VecScatter *newctx)
967: {
968: VecScatter ctx;
969: PetscErrorCode ierr;
970: PetscMPIInt size;
971: PetscInt xin_type = VEC_SEQ_ID,yin_type = VEC_SEQ_ID,*range;
972: PetscInt ix_type = IS_GENERAL_ID,iy_type = IS_GENERAL_ID;
973: MPI_Comm comm,ycomm;
974: PetscBool totalv,ixblock,iyblock,iystride,islocal,cando,flag;
975: IS tix = 0,tiy = 0;
978: if (!ix && !iy) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_SUP,"Cannot pass default in for both input and output indices");
980: /*
981: Determine if the vectors are "parallel", ie. it shares a comm with other processors, or
982: sequential (it does not share a comm). The difference is that parallel vectors treat the
983: index set as providing indices in the global parallel numbering of the vector, with
984: sequential vectors treat the index set as providing indices in the local sequential
985: numbering
986: */
987: PetscObjectGetComm((PetscObject)xin,&comm);
988: MPI_Comm_size(comm,&size);
989: if (size > 1) xin_type = VEC_MPI_ID;
991: PetscObjectGetComm((PetscObject)yin,&ycomm);
992: MPI_Comm_size(ycomm,&size);
993: if (size > 1) {comm = ycomm; yin_type = VEC_MPI_ID;}
995: /* generate the Scatter context */
996: PetscHeaderCreate(ctx,VEC_SCATTER_CLASSID,"VecScatter","VecScatter","Vec",comm,VecScatterDestroy,VecScatterView);
997: ctx->inuse = PETSC_FALSE;
999: ctx->beginandendtogether = PETSC_FALSE;
1000: PetscOptionsGetBool(NULL,"-vecscatter_merge",&ctx->beginandendtogether,NULL);
1001: if (ctx->beginandendtogether) {
1002: PetscInfo(ctx,"Using combined (merged) vector scatter begin and end\n");
1003: }
1004: ctx->packtogether = PETSC_FALSE;
1005: PetscOptionsGetBool(NULL,"-vecscatter_packtogether",&ctx->packtogether,NULL);
1006: if (ctx->packtogether) {
1007: PetscInfo(ctx,"Pack all messages before sending\n");
1008: }
1010: VecGetLocalSize(xin,&ctx->from_n);
1011: VecGetLocalSize(yin,&ctx->to_n);
1013: /*
1014: if ix or iy is not included; assume just grabbing entire vector
1015: */
1016: if (!ix && xin_type == VEC_SEQ_ID) {
1017: ISCreateStride(comm,ctx->from_n,0,1,&ix);
1018: tix = ix;
1019: } else if (!ix && xin_type == VEC_MPI_ID) {
1020: if (yin_type == VEC_MPI_ID) {
1021: PetscInt ntmp, low;
1022: VecGetLocalSize(xin,&ntmp);
1023: VecGetOwnershipRange(xin,&low,NULL);
1024: ISCreateStride(comm,ntmp,low,1,&ix);
1025: } else {
1026: PetscInt Ntmp;
1027: VecGetSize(xin,&Ntmp);
1028: ISCreateStride(comm,Ntmp,0,1,&ix);
1029: }
1030: tix = ix;
1031: } else if (!ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ix not given, but not Seq or MPI vector");
1033: if (!iy && yin_type == VEC_SEQ_ID) {
1034: ISCreateStride(comm,ctx->to_n,0,1,&iy);
1035: tiy = iy;
1036: } else if (!iy && yin_type == VEC_MPI_ID) {
1037: if (xin_type == VEC_MPI_ID) {
1038: PetscInt ntmp, low;
1039: VecGetLocalSize(yin,&ntmp);
1040: VecGetOwnershipRange(yin,&low,NULL);
1041: ISCreateStride(comm,ntmp,low,1,&iy);
1042: } else {
1043: PetscInt Ntmp;
1044: VecGetSize(yin,&Ntmp);
1045: ISCreateStride(comm,Ntmp,0,1,&iy);
1046: }
1047: tiy = iy;
1048: } else if (!iy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"iy not given, but not Seq or MPI vector");
1050: /*
1051: Determine types of index sets
1052: */
1053: PetscObjectTypeCompare((PetscObject)ix,ISBLOCK,&flag);
1054: if (flag) ix_type = IS_BLOCK_ID;
1055: PetscObjectTypeCompare((PetscObject)iy,ISBLOCK,&flag);
1056: if (flag) iy_type = IS_BLOCK_ID;
1057: PetscObjectTypeCompare((PetscObject)ix,ISSTRIDE,&flag);
1058: if (flag) ix_type = IS_STRIDE_ID;
1059: PetscObjectTypeCompare((PetscObject)iy,ISSTRIDE,&flag);
1060: if (flag) iy_type = IS_STRIDE_ID;
1062: /* ===========================================================================================================
1063: Check for special cases
1064: ==========================================================================================================*/
1065: /* ---------------------------------------------------------------------------*/
1066: if (xin_type == VEC_SEQ_ID && yin_type == VEC_SEQ_ID) {
1067: if (ix_type == IS_GENERAL_ID && iy_type == IS_GENERAL_ID) {
1068: PetscInt nx,ny;
1069: const PetscInt *idx,*idy;
1070: VecScatter_Seq_General *to = NULL,*from = NULL;
1072: ISGetLocalSize(ix,&nx);
1073: ISGetLocalSize(iy,&ny);
1074: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1075: ISGetIndices(ix,&idx);
1076: ISGetIndices(iy,&idy);
1077: PetscMalloc2(1,&to,1,&from);
1078: PetscMalloc2(nx,&to->vslots,nx,&from->vslots);
1079: to->n = nx;
1080: #if defined(PETSC_USE_DEBUG)
1081: VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1082: #endif
1083: PetscMemcpy(to->vslots,idy,nx*sizeof(PetscInt));
1084: from->n = nx;
1085: #if defined(PETSC_USE_DEBUG)
1086: VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1087: #endif
1088: PetscMemcpy(from->vslots,idx,nx*sizeof(PetscInt));
1089: to->type = VEC_SCATTER_SEQ_GENERAL;
1090: from->type = VEC_SCATTER_SEQ_GENERAL;
1091: ctx->todata = (void*)to;
1092: ctx->fromdata = (void*)from;
1093: ctx->ops->begin = VecScatterBegin_SGToSG;
1094: ctx->ops->end = 0;
1095: ctx->ops->destroy = VecScatterDestroy_SGToSG;
1096: ctx->ops->copy = VecScatterCopy_SGToSG;
1097: ctx->ops->view = VecScatterView_SGToSG;
1098: PetscInfo(xin,"Special case: sequential vector general scatter\n");
1099: goto functionend;
1100: } else if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1101: PetscInt nx,ny,to_first,to_step,from_first,from_step;
1102: VecScatter_Seq_Stride *from8 = NULL,*to8 = NULL;
1104: ISGetLocalSize(ix,&nx);
1105: ISGetLocalSize(iy,&ny);
1106: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1107: ISStrideGetInfo(iy,&to_first,&to_step);
1108: ISStrideGetInfo(ix,&from_first,&from_step);
1109: PetscMalloc2(1,&to8,1,&from8);
1110: to8->n = nx;
1111: to8->first = to_first;
1112: to8->step = to_step;
1113: from8->n = nx;
1114: from8->first = from_first;
1115: from8->step = from_step;
1116: to8->type = VEC_SCATTER_SEQ_STRIDE;
1117: from8->type = VEC_SCATTER_SEQ_STRIDE;
1118: ctx->todata = (void*)to8;
1119: ctx->fromdata = (void*)from8;
1120: ctx->ops->begin = VecScatterBegin_SSToSS;
1121: ctx->ops->end = 0;
1122: ctx->ops->destroy = VecScatterDestroy_SSToSS;
1123: ctx->ops->copy = VecScatterCopy_SSToSS;
1124: ctx->ops->view = VecScatterView_SSToSS;
1125: PetscInfo(xin,"Special case: sequential vector stride to stride\n");
1126: goto functionend;
1127: } else if (ix_type == IS_GENERAL_ID && iy_type == IS_STRIDE_ID) {
1128: PetscInt nx,ny,first,step;
1129: const PetscInt *idx;
1130: VecScatter_Seq_General *from9 = NULL;
1131: VecScatter_Seq_Stride *to9 = NULL;
1133: ISGetLocalSize(ix,&nx);
1134: ISGetIndices(ix,&idx);
1135: ISGetLocalSize(iy,&ny);
1136: ISStrideGetInfo(iy,&first,&step);
1137: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1138: PetscMalloc2(1,&to9,1,&from9);
1139: PetscMalloc1(nx,&from9->vslots);
1140: to9->n = nx;
1141: to9->first = first;
1142: to9->step = step;
1143: from9->n = nx;
1144: #if defined(PETSC_USE_DEBUG)
1145: VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1146: #endif
1147: PetscMemcpy(from9->vslots,idx,nx*sizeof(PetscInt));
1148: ctx->todata = (void*)to9; ctx->fromdata = (void*)from9;
1149: if (step == 1) ctx->ops->begin = VecScatterBegin_SGToSS_Stride1;
1150: else ctx->ops->begin = VecScatterBegin_SGToSS;
1151: ctx->ops->destroy = VecScatterDestroy_SGToSS;
1152: ctx->ops->end = 0;
1153: ctx->ops->copy = VecScatterCopy_SGToSS;
1154: ctx->ops->view = VecScatterView_SGToSS;
1155: to9->type = VEC_SCATTER_SEQ_STRIDE;
1156: from9->type = VEC_SCATTER_SEQ_GENERAL;
1157: PetscInfo(xin,"Special case: sequential vector general to stride\n");
1158: goto functionend;
1159: } else if (ix_type == IS_STRIDE_ID && iy_type == IS_GENERAL_ID) {
1160: PetscInt nx,ny,first,step;
1161: const PetscInt *idy;
1162: VecScatter_Seq_General *to10 = NULL;
1163: VecScatter_Seq_Stride *from10 = NULL;
1165: ISGetLocalSize(ix,&nx);
1166: ISGetIndices(iy,&idy);
1167: ISGetLocalSize(iy,&ny);
1168: ISStrideGetInfo(ix,&first,&step);
1169: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1170: PetscMalloc2(1,&to10,1,&from10);
1171: PetscMalloc1(nx,&to10->vslots);
1172: from10->n = nx;
1173: from10->first = first;
1174: from10->step = step;
1175: to10->n = nx;
1176: #if defined(PETSC_USE_DEBUG)
1177: VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1178: #endif
1179: PetscMemcpy(to10->vslots,idy,nx*sizeof(PetscInt));
1180: ctx->todata = (void*)to10;
1181: ctx->fromdata = (void*)from10;
1182: if (step == 1) ctx->ops->begin = VecScatterBegin_SSToSG_Stride1;
1183: else ctx->ops->begin = VecScatterBegin_SSToSG;
1184: ctx->ops->destroy = VecScatterDestroy_SSToSG;
1185: ctx->ops->end = 0;
1186: ctx->ops->copy = 0;
1187: ctx->ops->view = VecScatterView_SSToSG;
1188: to10->type = VEC_SCATTER_SEQ_GENERAL;
1189: from10->type = VEC_SCATTER_SEQ_STRIDE;
1190: PetscInfo(xin,"Special case: sequential vector stride to general\n");
1191: goto functionend;
1192: } else {
1193: PetscInt nx,ny;
1194: const PetscInt *idx,*idy;
1195: VecScatter_Seq_General *to11 = NULL,*from11 = NULL;
1196: PetscBool idnx,idny;
1198: ISGetLocalSize(ix,&nx);
1199: ISGetLocalSize(iy,&ny);
1200: if (nx != ny) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match, in %D out %D",nx,ny);
1202: ISIdentity(ix,&idnx);
1203: ISIdentity(iy,&idny);
1204: if (idnx && idny) {
1205: VecScatter_Seq_Stride *to13 = NULL,*from13 = NULL;
1206: PetscMalloc2(1,&to13,1,&from13);
1207: to13->n = nx;
1208: to13->first = 0;
1209: to13->step = 1;
1210: from13->n = nx;
1211: from13->first = 0;
1212: from13->step = 1;
1213: to13->type = VEC_SCATTER_SEQ_STRIDE;
1214: from13->type = VEC_SCATTER_SEQ_STRIDE;
1215: ctx->todata = (void*)to13;
1216: ctx->fromdata = (void*)from13;
1217: ctx->ops->begin = VecScatterBegin_SSToSS;
1218: ctx->ops->end = 0;
1219: ctx->ops->destroy = VecScatterDestroy_SSToSS;
1220: ctx->ops->copy = VecScatterCopy_SSToSS;
1221: ctx->ops->view = VecScatterView_SSToSS;
1222: PetscInfo(xin,"Special case: sequential copy\n");
1223: goto functionend;
1224: }
1226: ISGetIndices(iy,&idy);
1227: ISGetIndices(ix,&idx);
1228: PetscMalloc2(1,&to11,1,&from11);
1229: PetscMalloc2(nx,&to11->vslots,nx,&from11->vslots);
1230: to11->n = nx;
1231: #if defined(PETSC_USE_DEBUG)
1232: VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1233: #endif
1234: PetscMemcpy(to11->vslots,idy,nx*sizeof(PetscInt));
1235: from11->n = nx;
1236: #if defined(PETSC_USE_DEBUG)
1237: VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1238: #endif
1239: PetscMemcpy(from11->vslots,idx,nx*sizeof(PetscInt));
1240: to11->type = VEC_SCATTER_SEQ_GENERAL;
1241: from11->type = VEC_SCATTER_SEQ_GENERAL;
1242: ctx->todata = (void*)to11;
1243: ctx->fromdata = (void*)from11;
1244: ctx->ops->begin = VecScatterBegin_SGToSG;
1245: ctx->ops->end = 0;
1246: ctx->ops->destroy = VecScatterDestroy_SGToSG;
1247: ctx->ops->copy = VecScatterCopy_SGToSG;
1248: ctx->ops->view = VecScatterView_SGToSG;
1249: ISRestoreIndices(ix,&idx);
1250: ISRestoreIndices(iy,&idy);
1251: PetscInfo(xin,"Sequential vector scatter with block indices\n");
1252: goto functionend;
1253: }
1254: }
1255: /* ---------------------------------------------------------------------------*/
1256: if (xin_type == VEC_MPI_ID && yin_type == VEC_SEQ_ID) {
1258: /* ===========================================================================================================
1259: Check for special cases
1260: ==========================================================================================================*/
1261: islocal = PETSC_FALSE;
1262: /* special case extracting (subset of) local portion */
1263: if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1264: PetscInt nx,ny,to_first,to_step,from_first,from_step;
1265: PetscInt start,end,min,max;
1266: VecScatter_Seq_Stride *from12 = NULL,*to12 = NULL;
1268: VecGetOwnershipRange(xin,&start,&end);
1269: ISGetLocalSize(ix,&nx);
1270: ISStrideGetInfo(ix,&from_first,&from_step);
1271: ISGetLocalSize(iy,&ny);
1272: ISStrideGetInfo(iy,&to_first,&to_step);
1273: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1274: ISGetMinMax(ix,&min,&max);
1275: if (min >= start && max < end) islocal = PETSC_TRUE;
1276: else islocal = PETSC_FALSE;
1277: MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1278: if (cando) {
1279: PetscMalloc2(1,&to12,1,&from12);
1280: to12->n = nx;
1281: to12->first = to_first;
1282: to12->step = to_step;
1283: from12->n = nx;
1284: from12->first = from_first-start;
1285: from12->step = from_step;
1286: to12->type = VEC_SCATTER_SEQ_STRIDE;
1287: from12->type = VEC_SCATTER_SEQ_STRIDE;
1288: ctx->todata = (void*)to12;
1289: ctx->fromdata = (void*)from12;
1290: ctx->ops->begin = VecScatterBegin_SSToSS;
1291: ctx->ops->end = 0;
1292: ctx->ops->destroy = VecScatterDestroy_SSToSS;
1293: ctx->ops->copy = VecScatterCopy_SSToSS;
1294: ctx->ops->view = VecScatterView_SSToSS;
1295: PetscInfo(xin,"Special case: processors only getting local values\n");
1296: goto functionend;
1297: }
1298: } else {
1299: MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1300: }
1302: /* test for special case of all processors getting entire vector */
1303: /* contains check that PetscMPIInt can handle the sizes needed */
1304: totalv = PETSC_FALSE;
1305: if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1306: PetscInt i,nx,ny,to_first,to_step,from_first,from_step,N;
1307: PetscMPIInt *count = NULL,*displx;
1308: VecScatter_MPI_ToAll *sto = NULL;
1310: ISGetLocalSize(ix,&nx);
1311: ISStrideGetInfo(ix,&from_first,&from_step);
1312: ISGetLocalSize(iy,&ny);
1313: ISStrideGetInfo(iy,&to_first,&to_step);
1314: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1315: VecGetSize(xin,&N);
1316: if (nx != N) totalv = PETSC_FALSE;
1317: else if (from_first == 0 && from_step == 1 && from_first == to_first && from_step == to_step) totalv = PETSC_TRUE;
1318: else totalv = PETSC_FALSE;
1319: MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1321: #if defined(PETSC_USE_64BIT_INDICES)
1322: if (cando && (yin->map->N < PETSC_MPI_INT_MAX)) {
1323: #else
1324: if (cando) {
1325: #endif
1326: MPI_Comm_size(PetscObjectComm((PetscObject)ctx),&size);
1327: PetscMalloc3(1,&sto,size,&count,size,&displx);
1328: range = xin->map->range;
1329: for (i=0; i<size; i++) {
1330: PetscMPIIntCast(range[i+1] - range[i],count+i);
1331: PetscMPIIntCast(range[i],displx+i);
1332: }
1333: sto->count = count;
1334: sto->displx = displx;
1335: sto->work1 = 0;
1336: sto->work2 = 0;
1337: sto->type = VEC_SCATTER_MPI_TOALL;
1338: ctx->todata = (void*)sto;
1339: ctx->fromdata = 0;
1340: ctx->ops->begin = VecScatterBegin_MPI_ToAll;
1341: ctx->ops->end = 0;
1342: ctx->ops->destroy = VecScatterDestroy_MPI_ToAll;
1343: ctx->ops->copy = VecScatterCopy_MPI_ToAll;
1344: ctx->ops->view = VecScatterView_MPI_ToAll;
1345: PetscInfo(xin,"Special case: all processors get entire parallel vector\n");
1346: goto functionend;
1347: }
1348: } else {
1349: MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1350: }
1352: /* test for special case of processor 0 getting entire vector */
1353: /* contains check that PetscMPIInt can handle the sizes needed */
1354: totalv = PETSC_FALSE;
1355: if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1356: PetscInt i,nx,ny,to_first,to_step,from_first,from_step,N;
1357: PetscMPIInt rank,*count = NULL,*displx;
1358: VecScatter_MPI_ToAll *sto = NULL;
1360: PetscObjectGetComm((PetscObject)xin,&comm);
1361: MPI_Comm_rank(comm,&rank);
1362: ISGetLocalSize(ix,&nx);
1363: ISStrideGetInfo(ix,&from_first,&from_step);
1364: ISGetLocalSize(iy,&ny);
1365: ISStrideGetInfo(iy,&to_first,&to_step);
1366: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1367: if (!rank) {
1368: VecGetSize(xin,&N);
1369: if (nx != N) totalv = PETSC_FALSE;
1370: else if (from_first == 0 && from_step == 1 &&
1371: from_first == to_first && from_step == to_step) totalv = PETSC_TRUE;
1372: else totalv = PETSC_FALSE;
1373: } else {
1374: if (!nx) totalv = PETSC_TRUE;
1375: else totalv = PETSC_FALSE;
1376: }
1377: MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1379: #if defined(PETSC_USE_64BIT_INDICES)
1380: if (cando && (yin->map->N < PETSC_MPI_INT_MAX)) {
1381: #else
1382: if (cando) {
1383: #endif
1384: MPI_Comm_size(PetscObjectComm((PetscObject)ctx),&size);
1385: PetscMalloc3(1,&sto,size,&count,size,&displx);
1386: range = xin->map->range;
1387: for (i=0; i<size; i++) {
1388: PetscMPIIntCast(range[i+1] - range[i],count+i);
1389: PetscMPIIntCast(range[i],displx+i);
1390: }
1391: sto->count = count;
1392: sto->displx = displx;
1393: sto->work1 = 0;
1394: sto->work2 = 0;
1395: sto->type = VEC_SCATTER_MPI_TOONE;
1396: ctx->todata = (void*)sto;
1397: ctx->fromdata = 0;
1398: ctx->ops->begin = VecScatterBegin_MPI_ToOne;
1399: ctx->ops->end = 0;
1400: ctx->ops->destroy = VecScatterDestroy_MPI_ToAll;
1401: ctx->ops->copy = VecScatterCopy_MPI_ToAll;
1402: ctx->ops->view = VecScatterView_MPI_ToAll;
1403: PetscInfo(xin,"Special case: processor zero gets entire parallel vector, rest get none\n");
1404: goto functionend;
1405: }
1406: } else {
1407: MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1408: }
1410: PetscObjectTypeCompare((PetscObject)ix,ISBLOCK,&ixblock);
1411: PetscObjectTypeCompare((PetscObject)iy,ISBLOCK,&iyblock);
1412: PetscObjectTypeCompare((PetscObject)iy,ISSTRIDE,&iystride);
1413: if (ixblock) {
1414: /* special case block to block */
1415: if (iyblock) {
1416: PetscInt nx,ny,bsx,bsy;
1417: const PetscInt *idx,*idy;
1418: ISGetBlockSize(iy,&bsy);
1419: ISGetBlockSize(ix,&bsx);
1420: if (bsx == bsy && VecScatterOptimizedBS(bsx)) {
1421: ISBlockGetLocalSize(ix,&nx);
1422: ISBlockGetIndices(ix,&idx);
1423: ISBlockGetLocalSize(iy,&ny);
1424: ISBlockGetIndices(iy,&idy);
1425: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1426: VecScatterCreate_PtoS(nx,idx,ny,idy,xin,yin,bsx,ctx);
1427: ISBlockRestoreIndices(ix,&idx);
1428: ISBlockRestoreIndices(iy,&idy);
1429: PetscInfo(xin,"Special case: blocked indices\n");
1430: goto functionend;
1431: }
1432: /* special case block to stride */
1433: } else if (iystride) {
1434: PetscInt ystart,ystride,ysize,bsx;
1435: ISStrideGetInfo(iy,&ystart,&ystride);
1436: ISGetLocalSize(iy,&ysize);
1437: ISGetBlockSize(ix,&bsx);
1438: /* see if stride index set is equivalent to block index set */
1439: if (VecScatterOptimizedBS(bsx) && ((ystart % bsx) == 0) && (ystride == 1) && ((ysize % bsx) == 0)) {
1440: PetscInt nx,il,*idy;
1441: const PetscInt *idx;
1442: ISBlockGetLocalSize(ix,&nx);
1443: ISBlockGetIndices(ix,&idx);
1444: if (ysize != bsx*nx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1445: PetscMalloc1(nx,&idy);
1446: if (nx) {
1447: idy[0] = ystart/bsx;
1448: for (il=1; il<nx; il++) idy[il] = idy[il-1] + 1;
1449: }
1450: VecScatterCreate_PtoS(nx,idx,nx,idy,xin,yin,bsx,ctx);
1451: PetscFree(idy);
1452: ISBlockRestoreIndices(ix,&idx);
1453: PetscInfo(xin,"Special case: blocked indices to stride\n");
1454: goto functionend;
1455: }
1456: }
1457: }
1458: /* left over general case */
1459: {
1460: PetscInt nx,ny;
1461: const PetscInt *idx,*idy;
1462: ISGetLocalSize(ix,&nx);
1463: ISGetIndices(ix,&idx);
1464: ISGetLocalSize(iy,&ny);
1465: ISGetIndices(iy,&idy);
1466: if (nx != ny) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match (%d %d)",nx,ny);
1467: VecScatterCreate_PtoS(nx,idx,ny,idy,xin,yin,1,ctx);
1468: ISRestoreIndices(ix,&idx);
1469: ISRestoreIndices(iy,&idy);
1470: PetscInfo(xin,"General case: MPI to Seq\n");
1471: goto functionend;
1472: }
1473: }
1474: /* ---------------------------------------------------------------------------*/
1475: if (xin_type == VEC_SEQ_ID && yin_type == VEC_MPI_ID) {
1476: /* ===========================================================================================================
1477: Check for special cases
1478: ==========================================================================================================*/
1479: /* special case local copy portion */
1480: islocal = PETSC_FALSE;
1481: if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1482: PetscInt nx,ny,to_first,to_step,from_step,start,end,from_first,min,max;
1483: VecScatter_Seq_Stride *from = NULL,*to = NULL;
1485: VecGetOwnershipRange(yin,&start,&end);
1486: ISGetLocalSize(ix,&nx);
1487: ISStrideGetInfo(ix,&from_first,&from_step);
1488: ISGetLocalSize(iy,&ny);
1489: ISStrideGetInfo(iy,&to_first,&to_step);
1490: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1491: ISGetMinMax(iy,&min,&max);
1492: if (min >= start && max < end) islocal = PETSC_TRUE;
1493: else islocal = PETSC_FALSE;
1494: MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)yin));
1495: if (cando) {
1496: PetscMalloc2(1,&to,1,&from);
1497: to->n = nx;
1498: to->first = to_first-start;
1499: to->step = to_step;
1500: from->n = nx;
1501: from->first = from_first;
1502: from->step = from_step;
1503: to->type = VEC_SCATTER_SEQ_STRIDE;
1504: from->type = VEC_SCATTER_SEQ_STRIDE;
1505: ctx->todata = (void*)to;
1506: ctx->fromdata = (void*)from;
1507: ctx->ops->begin = VecScatterBegin_SSToSS;
1508: ctx->ops->end = 0;
1509: ctx->ops->destroy = VecScatterDestroy_SSToSS;
1510: ctx->ops->copy = VecScatterCopy_SSToSS;
1511: ctx->ops->view = VecScatterView_SSToSS;
1512: PetscInfo(xin,"Special case: sequential stride to MPI stride\n");
1513: goto functionend;
1514: }
1515: } else {
1516: MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)yin));
1517: }
1518: /* special case block to stride */
1519: if (ix_type == IS_BLOCK_ID && iy_type == IS_STRIDE_ID) {
1520: PetscInt ystart,ystride,ysize,bsx;
1521: ISStrideGetInfo(iy,&ystart,&ystride);
1522: ISGetLocalSize(iy,&ysize);
1523: ISGetBlockSize(ix,&bsx);
1524: /* see if stride index set is equivalent to block index set */
1525: if (VecScatterOptimizedBS(bsx) && ((ystart % bsx) == 0) && (ystride == 1) && ((ysize % bsx) == 0)) {
1526: PetscInt nx,il,*idy;
1527: const PetscInt *idx;
1528: ISBlockGetLocalSize(ix,&nx);
1529: ISBlockGetIndices(ix,&idx);
1530: if (ysize != bsx*nx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1531: PetscMalloc1(nx,&idy);
1532: if (nx) {
1533: idy[0] = ystart/bsx;
1534: for (il=1; il<nx; il++) idy[il] = idy[il-1] + 1;
1535: }
1536: VecScatterCreate_StoP(nx,idx,nx,idy,xin,yin,bsx,ctx);
1537: PetscFree(idy);
1538: ISBlockRestoreIndices(ix,&idx);
1539: PetscInfo(xin,"Special case: Blocked indices to stride\n");
1540: goto functionend;
1541: }
1542: }
1544: /* general case */
1545: {
1546: PetscInt nx,ny;
1547: const PetscInt *idx,*idy;
1548: ISGetLocalSize(ix,&nx);
1549: ISGetIndices(ix,&idx);
1550: ISGetLocalSize(iy,&ny);
1551: ISGetIndices(iy,&idy);
1552: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1553: VecScatterCreate_StoP(nx,idx,ny,idy,xin,yin,1,ctx);
1554: ISRestoreIndices(ix,&idx);
1555: ISRestoreIndices(iy,&idy);
1556: PetscInfo(xin,"General case: Seq to MPI\n");
1557: goto functionend;
1558: }
1559: }
1560: /* ---------------------------------------------------------------------------*/
1561: if (xin_type == VEC_MPI_ID && yin_type == VEC_MPI_ID) {
1562: /* no special cases for now */
1563: PetscInt nx,ny;
1564: const PetscInt *idx,*idy;
1565: ISGetLocalSize(ix,&nx);
1566: ISGetIndices(ix,&idx);
1567: ISGetLocalSize(iy,&ny);
1568: ISGetIndices(iy,&idy);
1569: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1570: VecScatterCreate_PtoP(nx,idx,ny,idy,xin,yin,1,ctx);
1571: ISRestoreIndices(ix,&idx);
1572: ISRestoreIndices(iy,&idy);
1573: PetscInfo(xin,"General case: MPI to MPI\n");
1574: goto functionend;
1575: }
1577: functionend:
1578: *newctx = ctx;
1579: ISDestroy(&tix);
1580: ISDestroy(&tiy);
1581: VecScatterViewFromOptions(ctx,NULL,"-vecscatter_view");
1582: return(0);
1583: }
1585: /* ------------------------------------------------------------------*/
1588: /*@
1589: VecScatterGetMerged - Returns true if the scatter is completed in the VecScatterBegin()
1590: and the VecScatterEnd() does nothing
1592: Not Collective
1594: Input Parameter:
1595: . ctx - scatter context created with VecScatterCreate()
1597: Output Parameter:
1598: . flg - PETSC_TRUE if the VecScatterBegin/End() are all done during the VecScatterBegin()
1600: Level: developer
1602: .seealso: VecScatterCreate(), VecScatterEnd(), VecScatterBegin()
1603: @*/
1604: PetscErrorCode VecScatterGetMerged(VecScatter ctx,PetscBool *flg)
1605: {
1608: *flg = ctx->beginandendtogether;
1609: return(0);
1610: }
1614: /*@
1615: VecScatterBegin - Begins a generalized scatter from one vector to
1616: another. Complete the scattering phase with VecScatterEnd().
1618: Neighbor-wise Collective on VecScatter and Vec
1620: Input Parameters:
1621: + inctx - scatter context generated by VecScatterCreate()
1622: . x - the vector from which we scatter
1623: . y - the vector to which we scatter
1624: . addv - either ADD_VALUES or INSERT_VALUES, with INSERT_VALUES mode any location
1625: not scattered to retains its old value; i.e. the vector is NOT first zeroed.
1626: - mode - the scattering mode, usually SCATTER_FORWARD. The available modes are:
1627: SCATTER_FORWARD or SCATTER_REVERSE
1630: Level: intermediate
1632: Options Database: See VecScatterCreate()
1634: Notes:
1635: The vectors x and y need not be the same vectors used in the call
1636: to VecScatterCreate(), but x must have the same parallel data layout
1637: as that passed in as the x to VecScatterCreate(), similarly for the y.
1638: Most likely they have been obtained from VecDuplicate().
1640: You cannot change the values in the input vector between the calls to VecScatterBegin()
1641: and VecScatterEnd().
1643: If you use SCATTER_REVERSE the two arguments x and y should be reversed, from
1644: the SCATTER_FORWARD.
1646: y[iy[i]] = x[ix[i]], for i=0,...,ni-1
1648: This scatter is far more general than the conventional
1649: scatter, since it can be a gather or a scatter or a combination,
1650: depending on the indices ix and iy. If x is a parallel vector and y
1651: is sequential, VecScatterBegin() can serve to gather values to a
1652: single processor. Similarly, if y is parallel and x sequential, the
1653: routine can scatter from one processor to many processors.
1655: Concepts: scatter^between vectors
1656: Concepts: gather^between vectors
1658: .seealso: VecScatterCreate(), VecScatterEnd()
1659: @*/
1660: PetscErrorCode VecScatterBegin(VecScatter inctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
1661: {
1663: #if defined(PETSC_USE_DEBUG)
1664: PetscInt to_n,from_n;
1665: #endif
1670: if (inctx->inuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE," Scatter ctx already in use");
1672: #if defined(PETSC_USE_DEBUG)
1673: /*
1674: Error checking to make sure these vectors match the vectors used
1675: to create the vector scatter context. -1 in the from_n and to_n indicate the
1676: vector lengths are unknown (for example with mapped scatters) and thus
1677: no error checking is performed.
1678: */
1679: if (inctx->from_n >= 0 && inctx->to_n >= 0) {
1680: VecGetLocalSize(x,&from_n);
1681: VecGetLocalSize(y,&to_n);
1682: if (mode & SCATTER_REVERSE) {
1683: if (to_n != inctx->from_n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter reverse and vector to != ctx from size)",to_n,inctx->from_n);
1684: if (from_n != inctx->to_n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter reverse and vector from != ctx to size)",from_n,inctx->to_n);
1685: } else {
1686: if (to_n != inctx->to_n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter forward and vector to != ctx to size)",to_n,inctx->to_n);
1687: if (from_n != inctx->from_n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter forward and vector from != ctx from size)",from_n,inctx->from_n);
1688: }
1689: }
1690: #endif
1692: inctx->inuse = PETSC_TRUE;
1693: PetscLogEventBarrierBegin(VEC_ScatterBarrier,0,0,0,0,PetscObjectComm((PetscObject)inctx));
1694: (*inctx->ops->begin)(inctx,x,y,addv,mode);
1695: if (inctx->beginandendtogether && inctx->ops->end) {
1696: inctx->inuse = PETSC_FALSE;
1697: (*inctx->ops->end)(inctx,x,y,addv,mode);
1698: }
1699: PetscLogEventBarrierEnd(VEC_ScatterBarrier,0,0,0,0,PetscObjectComm((PetscObject)inctx));
1700: return(0);
1701: }
1703: /* --------------------------------------------------------------------*/
1706: /*@
1707: VecScatterEnd - Ends a generalized scatter from one vector to another. Call
1708: after first calling VecScatterBegin().
1710: Neighbor-wise Collective on VecScatter and Vec
1712: Input Parameters:
1713: + ctx - scatter context generated by VecScatterCreate()
1714: . x - the vector from which we scatter
1715: . y - the vector to which we scatter
1716: . addv - either ADD_VALUES or INSERT_VALUES.
1717: - mode - the scattering mode, usually SCATTER_FORWARD. The available modes are:
1718: SCATTER_FORWARD, SCATTER_REVERSE
1720: Level: intermediate
1722: Notes:
1723: If you use SCATTER_REVERSE the arguments x and y should be reversed, from the SCATTER_FORWARD.
1725: y[iy[i]] = x[ix[i]], for i=0,...,ni-1
1727: .seealso: VecScatterBegin(), VecScatterCreate()
1728: @*/
1729: PetscErrorCode VecScatterEnd(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
1730: {
1737: ctx->inuse = PETSC_FALSE;
1738: if (!ctx->ops->end) return(0);
1739: if (!ctx->beginandendtogether) {
1740: PetscLogEventBegin(VEC_ScatterEnd,ctx,x,y,0);
1741: (*(ctx)->ops->end)(ctx,x,y,addv,mode);
1742: PetscLogEventEnd(VEC_ScatterEnd,ctx,x,y,0);
1743: }
1744: return(0);
1745: }
1749: /*@C
1750: VecScatterDestroy - Destroys a scatter context created by
1751: VecScatterCreate().
1753: Collective on VecScatter
1755: Input Parameter:
1756: . ctx - the scatter context
1758: Level: intermediate
1760: .seealso: VecScatterCreate(), VecScatterCopy()
1761: @*/
1762: PetscErrorCode VecScatterDestroy(VecScatter *ctx)
1763: {
1767: if (!*ctx) return(0);
1769: if ((*ctx)->inuse) SETERRQ(((PetscObject)(*ctx))->comm,PETSC_ERR_ARG_WRONGSTATE,"Scatter context is in use");
1770: if (--((PetscObject)(*ctx))->refct > 0) {*ctx = 0; return(0);}
1772: /* if memory was published with SAWs then destroy it */
1773: PetscObjectSAWsViewOff((PetscObject)(*ctx));
1774: (*(*ctx)->ops->destroy)(*ctx);
1775: #if defined(PETSC_HAVE_CUSP)
1776: VecScatterCUSPIndicesDestroy((PetscCUSPIndices*)&((*ctx)->spptr));
1777: #endif
1778: PetscHeaderDestroy(ctx);
1779: return(0);
1780: }
1784: /*@
1785: VecScatterCopy - Makes a copy of a scatter context.
1787: Collective on VecScatter
1789: Input Parameter:
1790: . sctx - the scatter context
1792: Output Parameter:
1793: . ctx - the context copy
1795: Level: advanced
1797: .seealso: VecScatterCreate(), VecScatterDestroy()
1798: @*/
1799: PetscErrorCode VecScatterCopy(VecScatter sctx,VecScatter *ctx)
1800: {
1806: if (!sctx->ops->copy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot copy this type");
1807: PetscHeaderCreate(*ctx,VEC_SCATTER_CLASSID,"VecScatter","VecScatter","Vec",PetscObjectComm((PetscObject)sctx),VecScatterDestroy,VecScatterView);
1808: (*ctx)->to_n = sctx->to_n;
1809: (*ctx)->from_n = sctx->from_n;
1810: (*sctx->ops->copy)(sctx,*ctx);
1811: return(0);
1812: }
1815: /* ------------------------------------------------------------------*/
1818: /*@
1819: VecScatterView - Views a vector scatter context.
1821: Collective on VecScatter
1823: Input Parameters:
1824: + ctx - the scatter context
1825: - viewer - the viewer for displaying the context
1827: Level: intermediate
1829: @*/
1830: PetscErrorCode VecScatterView(VecScatter ctx,PetscViewer viewer)
1831: {
1836: if (!viewer) {
1837: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ctx),&viewer);
1838: }
1840: if (ctx->ops->view) {
1841: (*ctx->ops->view)(ctx,viewer);
1842: }
1843: return(0);
1844: }
1848: /*@C
1849: VecScatterRemap - Remaps the "from" and "to" indices in a
1850: vector scatter context. FOR EXPERTS ONLY!
1852: Collective on VecScatter
1854: Input Parameters:
1855: + scat - vector scatter context
1856: . from - remapping for "from" indices (may be NULL)
1857: - to - remapping for "to" indices (may be NULL)
1859: Level: developer
1861: Notes: In the parallel case the todata is actually the indices
1862: from which the data is TAKEN! The from stuff is where the
1863: data is finally put. This is VERY VERY confusing!
1865: In the sequential case the todata is the indices where the
1866: data is put and the fromdata is where it is taken from.
1867: This is backwards from the paralllel case! CRY! CRY! CRY!
1869: @*/
1870: PetscErrorCode VecScatterRemap(VecScatter scat,PetscInt *rto,PetscInt *rfrom)
1871: {
1872: VecScatter_Seq_General *to,*from;
1873: VecScatter_MPI_General *mto;
1874: PetscInt i;
1881: from = (VecScatter_Seq_General*)scat->fromdata;
1882: mto = (VecScatter_MPI_General*)scat->todata;
1884: if (mto->type == VEC_SCATTER_MPI_TOALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Not for to all scatter");
1886: if (rto) {
1887: if (mto->type == VEC_SCATTER_MPI_GENERAL) {
1888: /* handle off processor parts */
1889: for (i=0; i<mto->starts[mto->n]; i++) mto->indices[i] = rto[mto->indices[i]];
1891: /* handle local part */
1892: to = &mto->local;
1893: for (i=0; i<to->n; i++) to->vslots[i] = rto[to->vslots[i]];
1894: } else if (from->type == VEC_SCATTER_SEQ_GENERAL) {
1895: for (i=0; i<from->n; i++) from->vslots[i] = rto[from->vslots[i]];
1896: } else if (from->type == VEC_SCATTER_SEQ_STRIDE) {
1897: VecScatter_Seq_Stride *sto = (VecScatter_Seq_Stride*)from;
1899: /* if the remapping is the identity and stride is identity then skip remap */
1900: if (sto->step == 1 && sto->first == 0) {
1901: for (i=0; i<sto->n; i++) {
1902: if (rto[i] != i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1903: }
1904: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1905: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1906: }
1908: if (rfrom) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unable to remap the FROM in scatters yet");
1910: /*
1911: Mark then vector lengths as unknown because we do not know the
1912: lengths of the remapped vectors
1913: */
1914: scat->from_n = -1;
1915: scat->to_n = -1;
1916: return(0);
1917: }
1921: /*
1922: VecScatterGetTypes_Private - Returns the scatter types.
1924: scatter - The scatter.
1925: from - Upon exit this contains the type of the from scatter.
1926: to - Upon exit this contains the type of the to scatter.
1927: */
1928: PetscErrorCode VecScatterGetTypes_Private(VecScatter scatter,VecScatterType *from,VecScatterType *to)
1929: {
1930: VecScatter_Common* fromdata = (VecScatter_Common*)scatter->fromdata;
1931: VecScatter_Common* todata = (VecScatter_Common*)scatter->todata;
1934: *from = fromdata->type;
1935: *to = todata->type;
1936: return(0);
1937: }
1942: /*
1943: VecScatterIsSequential_Private - Returns true if the scatter is sequential.
1945: scatter - The scatter.
1946: flag - Upon exit flag is true if the scatter is of type VecScatter_Seq_General
1947: or VecScatter_Seq_Stride; otherwise flag is false.
1948: */
1949: PetscErrorCode VecScatterIsSequential_Private(VecScatter_Common *scatter,PetscBool *flag)
1950: {
1951: VecScatterType scatterType = scatter->type;
1954: if (scatterType == VEC_SCATTER_SEQ_GENERAL || scatterType == VEC_SCATTER_SEQ_STRIDE) {
1955: *flag = PETSC_TRUE;
1956: } else {
1957: *flag = PETSC_FALSE;
1958: }
1959: return(0);
1960: }