Actual source code: vscat.c
petsc-3.5.4 2015-05-23
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: VecGetArray(y,&yv);
56: if (x != y) {VecGetArray(x,&xv);}
58: if (mode & SCATTER_REVERSE) {
59: PetscScalar *xvt,*xvt2;
60: VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
61: PetscInt i;
62: PetscMPIInt *disply = scat->displx;
64: if (addv == INSERT_VALUES) {
65: PetscInt rstart,rend;
66: /*
67: copy the correct part of the local vector into the local storage of
68: the MPI one. Note: this operation only makes sense if all the local
69: vectors have the same values
70: */
71: VecGetOwnershipRange(y,&rstart,&rend);
72: PetscMemcpy(yv,xv+rstart,yy_n*sizeof(PetscScalar));
73: } else {
74: MPI_Comm comm;
75: PetscMPIInt rank;
76: PetscObjectGetComm((PetscObject)y,&comm);
77: MPI_Comm_rank(comm,&rank);
78: if (scat->work1) xvt = scat->work1;
79: else {
80: PetscMalloc1(xx_n,&xvt);
81: scat->work1 = xvt;
82: }
83: if (!rank) { /* I am the zeroth processor, values are accumulated here */
84: if (scat->work2) xvt2 = scat->work2;
85: else {
86: PetscMalloc1(xx_n,&xvt2);
87: scat->work2 = xvt2;
88: }
89: MPI_Gatherv(yv,yy_n,MPIU_SCALAR,xvt2,scat->count,disply,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
90: MPI_Reduce(xv,xvt,xx_n,MPIU_SCALAR,MPIU_SUM,0,PetscObjectComm((PetscObject)ctx));
91: if (addv == ADD_VALUES) {
92: for (i=0; i<xx_n; i++) xvt[i] += xvt2[i];
93: #if !defined(PETSC_USE_COMPLEX)
94: } else if (addv == MAX_VALUES) {
95: for (i=0; i<xx_n; i++) xvt[i] = PetscMax(xvt[i],xvt2[i]);
96: #endif
97: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
98: MPI_Scatterv(xvt,scat->count,disply,MPIU_SCALAR,yv,yy_n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
99: } else {
100: MPI_Gatherv(yv,yy_n,MPIU_SCALAR,0, 0,0,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
101: MPI_Reduce(xv,xvt,xx_n,MPIU_SCALAR,MPIU_SUM,0,PetscObjectComm((PetscObject)ctx));
102: MPI_Scatterv(0,scat->count,disply,MPIU_SCALAR,yv,yy_n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
103: }
104: }
105: } else {
106: PetscScalar *yvt;
107: VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
108: PetscInt i;
109: PetscMPIInt *displx = scat->displx;
111: if (addv == INSERT_VALUES) {
112: MPI_Allgatherv(xv,xx_n,MPIU_SCALAR,yv,scat->count,displx,MPIU_SCALAR,PetscObjectComm((PetscObject)ctx));
113: } else {
114: if (scat->work1) yvt = scat->work1;
115: else {
116: PetscMalloc1(yy_n,&yvt);
117: scat->work1 = yvt;
118: }
119: MPI_Allgatherv(xv,xx_n,MPIU_SCALAR,yvt,scat->count,displx,MPIU_SCALAR,PetscObjectComm((PetscObject)ctx));
120: if (addv == ADD_VALUES) {
121: for (i=0; i<yy_n; i++) yv[i] += yvt[i];
122: #if !defined(PETSC_USE_COMPLEX)
123: } else if (addv == MAX_VALUES) {
124: for (i=0; i<yy_n; i++) yv[i] = PetscMax(yv[i],yvt[i]);
125: #endif
126: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
127: }
128: }
129: VecRestoreArray(y,&yv);
130: if (x != y) {VecRestoreArray(x,&xv);}
131: return(0);
132: }
136: PetscErrorCode VecScatterView_MPI_ToAll(VecScatter in,PetscViewer viewer)
137: {
139: PetscBool isascii;
142: PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
143: if (isascii) {
144: PetscViewerASCIIPrintf(viewer,"Entire parallel vector is copied to each process\n");
145: }
146: return(0);
147: }
149: /*
150: This is special scatter code for when the entire parallel vector is copied to processor 0.
152: */
155: PetscErrorCode VecScatterBegin_MPI_ToOne(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
156: {
158: PetscMPIInt rank;
159: PetscInt yy_n,xx_n;
160: PetscScalar *xv,*yv;
161: MPI_Comm comm;
164: VecGetLocalSize(y,&yy_n);
165: VecGetLocalSize(x,&xx_n);
166: VecGetArray(x,&xv);
167: VecGetArray(y,&yv);
169: PetscObjectGetComm((PetscObject)x,&comm);
170: MPI_Comm_rank(comm,&rank);
172: /* -------- Reverse scatter; spread from processor 0 to other processors */
173: if (mode & SCATTER_REVERSE) {
174: PetscScalar *yvt;
175: VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
176: PetscInt i;
177: PetscMPIInt *disply = scat->displx;
179: if (addv == INSERT_VALUES) {
180: MPI_Scatterv(xv,scat->count,disply,MPIU_SCALAR,yv,yy_n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
181: } else {
182: if (scat->work2) yvt = scat->work2;
183: else {
184: PetscMalloc1(xx_n,&yvt);
185: scat->work2 = yvt;
186: }
187: MPI_Scatterv(xv,scat->count,disply,MPIU_SCALAR,yvt,yy_n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
188: if (addv == ADD_VALUES) {
189: for (i=0; i<yy_n; i++) yv[i] += yvt[i];
190: #if !defined(PETSC_USE_COMPLEX)
191: } else if (addv == MAX_VALUES) {
192: for (i=0; i<yy_n; i++) yv[i] = PetscMax(yv[i],yvt[i]);
193: #endif
194: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
195: }
196: /* --------- Forward scatter; gather all values onto processor 0 */
197: } else {
198: PetscScalar *yvt = 0;
199: VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
200: PetscInt i;
201: PetscMPIInt *displx = scat->displx;
203: if (addv == INSERT_VALUES) {
204: MPI_Gatherv(xv,xx_n,MPIU_SCALAR,yv,scat->count,displx,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
205: } else {
206: if (!rank) {
207: if (scat->work1) yvt = scat->work1;
208: else {
209: PetscMalloc1(yy_n,&yvt);
210: scat->work1 = yvt;
211: }
212: }
213: MPI_Gatherv(xv,xx_n,MPIU_SCALAR,yvt,scat->count,displx,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
214: if (!rank) {
215: if (addv == ADD_VALUES) {
216: for (i=0; i<yy_n; i++) yv[i] += yvt[i];
217: #if !defined(PETSC_USE_COMPLEX)
218: } else if (addv == MAX_VALUES) {
219: for (i=0; i<yy_n; i++) yv[i] = PetscMax(yv[i],yvt[i]);
220: #endif
221: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
222: }
223: }
224: }
225: VecRestoreArray(x,&xv);
226: VecRestoreArray(y,&yv);
227: return(0);
228: }
230: /*
231: The follow to are used for both VecScatterBegin_MPI_ToAll() and VecScatterBegin_MPI_ToOne()
232: */
235: PetscErrorCode VecScatterDestroy_MPI_ToAll(VecScatter ctx)
236: {
237: VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
238: PetscErrorCode ierr;
241: PetscFree(scat->work1);
242: PetscFree(scat->work2);
243: PetscFree3(ctx->todata,scat->count,scat->displx);
244: return(0);
245: }
249: PetscErrorCode VecScatterDestroy_SGToSG(VecScatter ctx)
250: {
254: PetscFree2(((VecScatter_Seq_General*)ctx->todata)->vslots,((VecScatter_Seq_General*)ctx->fromdata)->vslots);
255: PetscFree2(ctx->todata,ctx->fromdata);
256: return(0);
257: }
261: PetscErrorCode VecScatterDestroy_SGToSS(VecScatter ctx)
262: {
266: PetscFree(((VecScatter_Seq_General*)ctx->fromdata)->vslots);
267: PetscFree2(ctx->todata,ctx->fromdata);
268: return(0);
269: }
273: PetscErrorCode VecScatterDestroy_SSToSG(VecScatter ctx)
274: {
278: PetscFree(((VecScatter_Seq_General*)ctx->todata)->vslots);
279: PetscFree2(ctx->todata,ctx->fromdata);
280: return(0);
281: }
285: PetscErrorCode VecScatterDestroy_SSToSS(VecScatter ctx)
286: {
290: PetscFree2(ctx->todata,ctx->fromdata);
291: return(0);
292: }
294: /* -------------------------------------------------------------------------*/
297: PetscErrorCode VecScatterCopy_MPI_ToAll(VecScatter in,VecScatter out)
298: {
299: VecScatter_MPI_ToAll *in_to = (VecScatter_MPI_ToAll*)in->todata,*sto;
300: PetscErrorCode ierr;
301: PetscMPIInt size,*count,*displx;
302: PetscInt i;
305: out->begin = in->begin;
306: out->end = in->end;
307: out->copy = in->copy;
308: out->destroy = in->destroy;
309: out->view = in->view;
311: MPI_Comm_size(PetscObjectComm((PetscObject)out),&size);
312: PetscMalloc3(1,&sto,size,&count,size,&displx);
313: sto->type = in_to->type;
314: sto->count = count;
315: sto->displx = displx;
316: for (i=0; i<size; i++) {
317: sto->count[i] = in_to->count[i];
318: sto->displx[i] = in_to->displx[i];
319: }
320: sto->work1 = 0;
321: sto->work2 = 0;
322: out->todata = (void*)sto;
323: out->fromdata = (void*)0;
324: return(0);
325: }
327: /* --------------------------------------------------------------------------------------*/
328: /*
329: Scatter: sequential general to sequential general
330: */
333: PetscErrorCode VecScatterBegin_SGToSG(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
334: {
335: VecScatter_Seq_General *gen_to = (VecScatter_Seq_General*)ctx->todata;
336: VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
337: PetscErrorCode ierr;
338: PetscInt i,n = gen_from->n,*fslots,*tslots;
339: PetscScalar *xv,*yv;
342: #if defined(PETSC_HAVE_CUSP)
343: if (x->valid_GPU_array == PETSC_CUSP_GPU) {
344: /* create the scatter indices if not done already */
345: if (!ctx->spptr) {
346: PetscInt tofirst = 0,tostep = 0,fromfirst = 0,fromstep = 0;
347: fslots = gen_from->vslots;
348: tslots = gen_to->vslots;
349: VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
350: }
351: /* next do the scatter */
352: VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
353: return(0);
354: }
355: #endif
357: VecGetArray(x,&xv);
358: if (x != y) {
359: VecGetArray(y,&yv);
360: } else yv = xv;
362: if (mode & SCATTER_REVERSE) {
363: gen_to = (VecScatter_Seq_General*)ctx->fromdata;
364: gen_from = (VecScatter_Seq_General*)ctx->todata;
365: }
366: fslots = gen_from->vslots;
367: tslots = gen_to->vslots;
369: if (addv == INSERT_VALUES) {
370: for (i=0; i<n; i++) yv[tslots[i]] = xv[fslots[i]];
371: } else if (addv == ADD_VALUES) {
372: for (i=0; i<n; i++) yv[tslots[i]] += xv[fslots[i]];
373: #if !defined(PETSC_USE_COMPLEX)
374: } else if (addv == MAX_VALUES) {
375: for (i=0; i<n; i++) yv[tslots[i]] = PetscMax(yv[tslots[i]],xv[fslots[i]]);
376: #endif
377: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
378: VecRestoreArray(x,&xv);
379: if (x != y) {VecRestoreArray(y,&yv);}
380: return(0);
381: }
383: /*
384: Scatter: sequential general to sequential stride 1
385: */
388: PetscErrorCode VecScatterBegin_SGToSS_Stride1(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
389: {
390: VecScatter_Seq_Stride *gen_to = (VecScatter_Seq_Stride*)ctx->todata;
391: VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
392: PetscInt i,n = gen_from->n,*fslots = gen_from->vslots;
393: PetscErrorCode ierr;
394: PetscInt first = gen_to->first;
395: PetscScalar *xv,*yv;
398: #if defined(PETSC_HAVE_CUSP)
399: if (x->valid_GPU_array == PETSC_CUSP_GPU) {
400: /* create the scatter indices if not done already */
401: if (!ctx->spptr) {
402: PetscInt tofirst = first,tostep = 1,fromfirst = 0,fromstep = 0;
403: PetscInt *tslots = 0;
404: VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
405: }
406: /* next do the scatter */
407: VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
408: return(0);
409: }
410: #endif
412: VecGetArray(x,&xv);
413: if (x != y) {
414: VecGetArray(y,&yv);
415: } else yv = xv;
416: if (mode & SCATTER_REVERSE) {
417: xv += first;
418: if (addv == INSERT_VALUES) {
419: for (i=0; i<n; i++) yv[fslots[i]] = xv[i];
420: } else if (addv == ADD_VALUES) {
421: for (i=0; i<n; i++) yv[fslots[i]] += xv[i];
422: #if !defined(PETSC_USE_COMPLEX)
423: } else if (addv == MAX_VALUES) {
424: for (i=0; i<n; i++) yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[i]);
425: #endif
426: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
427: } else {
428: yv += first;
429: if (addv == INSERT_VALUES) {
430: for (i=0; i<n; i++) yv[i] = xv[fslots[i]];
431: } else if (addv == ADD_VALUES) {
432: for (i=0; i<n; i++) yv[i] += xv[fslots[i]];
433: #if !defined(PETSC_USE_COMPLEX)
434: } else if (addv == MAX_VALUES) {
435: for (i=0; i<n; i++) yv[i] = PetscMax(yv[i],xv[fslots[i]]);
436: #endif
437: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
438: }
439: VecRestoreArray(x,&xv);
440: if (x != y) {VecRestoreArray(y,&yv);}
441: return(0);
442: }
444: /*
445: Scatter: sequential general to sequential stride
446: */
449: PetscErrorCode VecScatterBegin_SGToSS(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
450: {
451: VecScatter_Seq_Stride *gen_to = (VecScatter_Seq_Stride*)ctx->todata;
452: VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
453: PetscInt i,n = gen_from->n,*fslots = gen_from->vslots;
454: PetscErrorCode ierr;
455: PetscInt first = gen_to->first,step = gen_to->step;
456: PetscScalar *xv,*yv;
459: #if defined(PETSC_HAVE_CUSP)
460: if (x->valid_GPU_array == PETSC_CUSP_GPU) {
461: /* create the scatter indices if not done already */
462: if (!ctx->spptr) {
463: PetscInt tofirst = first,tostep = step,fromfirst = 0,fromstep = 0;
464: PetscInt * tslots = 0;
465: VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
466: }
467: /* next do the scatter */
468: VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
469: return(0);
470: }
471: #endif
473: VecGetArray(x,&xv);
474: if (x != y) {
475: VecGetArray(y,&yv);
476: } else yv = xv;
478: if (mode & SCATTER_REVERSE) {
479: if (addv == INSERT_VALUES) {
480: for (i=0; i<n; i++) yv[fslots[i]] = xv[first + i*step];
481: } else if (addv == ADD_VALUES) {
482: for (i=0; i<n; i++) yv[fslots[i]] += xv[first + i*step];
483: #if !defined(PETSC_USE_COMPLEX)
484: } else if (addv == MAX_VALUES) {
485: for (i=0; i<n; i++) yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[first + i*step]);
486: #endif
487: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
488: } else {
489: if (addv == INSERT_VALUES) {
490: for (i=0; i<n; i++) yv[first + i*step] = xv[fslots[i]];
491: } else if (addv == ADD_VALUES) {
492: for (i=0; i<n; i++) yv[first + i*step] += xv[fslots[i]];
493: #if !defined(PETSC_USE_COMPLEX)
494: } else if (addv == MAX_VALUES) {
495: for (i=0; i<n; i++) yv[first + i*step] = PetscMax(yv[first + i*step],xv[fslots[i]]);
496: #endif
497: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
498: }
499: VecRestoreArray(x,&xv);
500: if (x != y) {VecRestoreArray(y,&yv);}
501: return(0);
502: }
504: /*
505: Scatter: sequential stride 1 to sequential general
506: */
509: PetscErrorCode VecScatterBegin_SSToSG_Stride1(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
510: {
511: VecScatter_Seq_Stride *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
512: VecScatter_Seq_General *gen_to = (VecScatter_Seq_General*)ctx->todata;
513: PetscInt i,n = gen_from->n,*fslots = gen_to->vslots;
514: PetscErrorCode ierr;
515: PetscInt first = gen_from->first;
516: PetscScalar *xv,*yv;
519: #if defined(PETSC_HAVE_CUSP)
520: if (x->valid_GPU_array == PETSC_CUSP_GPU) {
521: /* create the scatter indices if not done already */
522: if (!ctx->spptr) {
523: PetscInt tofirst = 0,tostep = 0,fromfirst = first,fromstep = 1;
524: PetscInt * tslots = 0;
525: VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
526: }
527: /* next do the scatter */
528: VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
529: return(0);
530: }
531: #endif
533: VecGetArray(x,&xv);
534: if (x != y) {
535: VecGetArray(y,&yv);
536: } else yv = xv;
538: if (mode & SCATTER_REVERSE) {
539: yv += first;
540: if (addv == INSERT_VALUES) {
541: for (i=0; i<n; i++) yv[i] = xv[fslots[i]];
542: } else if (addv == ADD_VALUES) {
543: for (i=0; i<n; i++) yv[i] += xv[fslots[i]];
544: #if !defined(PETSC_USE_COMPLEX)
545: } else if (addv == MAX_VALUES) {
546: for (i=0; i<n; i++) yv[i] = PetscMax(yv[i],xv[fslots[i]]);
547: #endif
548: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
549: } else {
550: xv += first;
551: if (addv == INSERT_VALUES) {
552: for (i=0; i<n; i++) yv[fslots[i]] = xv[i];
553: } else if (addv == ADD_VALUES) {
554: for (i=0; i<n; i++) yv[fslots[i]] += xv[i];
555: #if !defined(PETSC_USE_COMPLEX)
556: } else if (addv == MAX_VALUES) {
557: for (i=0; i<n; i++) yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[i]);
558: #endif
559: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
560: }
561: VecRestoreArray(x,&xv);
562: if (x != y) {VecRestoreArray(y,&yv);}
563: return(0);
564: }
568: /*
569: Scatter: sequential stride to sequential general
570: */
571: PetscErrorCode VecScatterBegin_SSToSG(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
572: {
573: VecScatter_Seq_Stride *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
574: VecScatter_Seq_General *gen_to = (VecScatter_Seq_General*)ctx->todata;
575: PetscInt i,n = gen_from->n,*fslots = gen_to->vslots;
576: PetscErrorCode ierr;
577: PetscInt first = gen_from->first,step = gen_from->step;
578: PetscScalar *xv,*yv;
581: #if defined(PETSC_HAVE_CUSP)
582: if (x->valid_GPU_array == PETSC_CUSP_GPU) {
583: /* create the scatter indices if not done already */
584: if (!ctx->spptr) {
585: PetscInt tofirst = 0,tostep = 0,fromfirst = first,fromstep = step;
586: PetscInt * tslots = 0;
587: VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,fslots,tslots,(PetscCUSPIndices*)&(ctx->spptr));
588: }
589: /* next do the scatter */
590: VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
591: return(0);
592: }
593: #endif
595: VecGetArray(x,&xv);
596: if (x != y) {
597: VecGetArray(y,&yv);
598: } else yv = xv;
600: if (mode & SCATTER_REVERSE) {
601: if (addv == INSERT_VALUES) {
602: for (i=0; i<n; i++) yv[first + i*step] = xv[fslots[i]];
603: } else if (addv == ADD_VALUES) {
604: for (i=0; i<n; i++) yv[first + i*step] += xv[fslots[i]];
605: #if !defined(PETSC_USE_COMPLEX)
606: } else if (addv == MAX_VALUES) {
607: for (i=0; i<n; i++) yv[first + i*step] = PetscMax(yv[first + i*step],xv[fslots[i]]);
608: #endif
609: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
610: } else {
611: if (addv == INSERT_VALUES) {
612: for (i=0; i<n; i++) yv[fslots[i]] = xv[first + i*step];
613: } else if (addv == ADD_VALUES) {
614: for (i=0; i<n; i++) yv[fslots[i]] += xv[first + i*step];
615: #if !defined(PETSC_USE_COMPLEX)
616: } else if (addv == MAX_VALUES) {
617: for (i=0; i<n; i++) yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[first + i*step]);
618: #endif
619: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
620: }
621: VecRestoreArray(x,&xv);
622: if (x != y) {VecRestoreArray(y,&yv);}
623: return(0);
624: }
628: PetscErrorCode VecScatterView_SSToSG(VecScatter in,PetscViewer viewer)
629: {
630: PetscErrorCode ierr;
631: VecScatter_Seq_Stride *in_from = (VecScatter_Seq_Stride*)in->todata;
632: VecScatter_Seq_General *in_to = (VecScatter_Seq_General*)in->fromdata;
633: PetscInt i;
634: PetscBool isascii;
637: PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
638: if (isascii) {
639: PetscViewerASCIIPrintf(viewer,"Sequential stride to general scatter\n");
640: for (i=0; i<in_to->n; i++) {
641: PetscViewerASCIIPrintf(viewer,"%D to %D\n",in_from->first + in_from->step*i,in_to->vslots[i]);
642: }
643: }
644: return(0);
645: }
646: /*
647: Scatter: sequential stride to sequential stride
648: */
651: PetscErrorCode VecScatterBegin_SSToSS(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
652: {
653: VecScatter_Seq_Stride *gen_to = (VecScatter_Seq_Stride*)ctx->todata;
654: VecScatter_Seq_Stride *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
655: PetscInt i,n = gen_from->n,to_first = gen_to->first,to_step = gen_to->step;
656: PetscErrorCode ierr;
657: PetscInt from_first = gen_from->first,from_step = gen_from->step;
658: PetscScalar *xv,*yv;
661: #if defined(PETSC_HAVE_CUSP)
662: if (x->valid_GPU_array == PETSC_CUSP_GPU) {
663: /* create the scatter indices if not done already */
664: if (!ctx->spptr) {
665: PetscInt *tslots = 0,*fslots = 0;
666: VecScatterCUSPIndicesCreate_StoS(n,to_first,from_first,to_step,from_step,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
667: }
668: /* next do the scatter */
669: VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
670: return(0);
671: }
672: #endif
674: VecGetArrayRead(x,(const PetscScalar **)&xv);
675: if (x != y) {
676: VecGetArray(y,&yv);
677: } else yv = xv;
679: if (mode & SCATTER_REVERSE) {
680: from_first = gen_to->first;
681: to_first = gen_from->first;
682: from_step = gen_to->step;
683: to_step = gen_from->step;
684: }
686: if (addv == INSERT_VALUES) {
687: if (to_step == 1 && from_step == 1) {
688: PetscMemcpy(yv+to_first,xv+from_first,n*sizeof(PetscScalar));
689: } else {
690: for (i=0; i<n; i++) yv[to_first + i*to_step] = xv[from_first+i*from_step];
691: }
692: } else if (addv == ADD_VALUES) {
693: if (to_step == 1 && from_step == 1) {
694: yv += to_first; xv += from_first;
695: for (i=0; i<n; i++) yv[i] += xv[i];
696: } else {
697: for (i=0; i<n; i++) yv[to_first + i*to_step] += xv[from_first+i*from_step];
698: }
699: #if !defined(PETSC_USE_COMPLEX)
700: } else if (addv == MAX_VALUES) {
701: if (to_step == 1 && from_step == 1) {
702: yv += to_first; xv += from_first;
703: for (i=0; i<n; i++) yv[i] = PetscMax(yv[i],xv[i]);
704: } else {
705: 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]);
706: }
707: #endif
708: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
709: VecRestoreArrayRead(x,(const PetscScalar**)&xv);
710: if (x != y) {VecRestoreArray(y,&yv);}
711: return(0);
712: }
714: /* --------------------------------------------------------------------------*/
719: PetscErrorCode VecScatterCopy_SGToSG(VecScatter in,VecScatter out)
720: {
721: PetscErrorCode ierr;
722: VecScatter_Seq_General *in_to = (VecScatter_Seq_General*)in->todata,*out_to = NULL;
723: VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from = NULL;
726: out->begin = in->begin;
727: out->end = in->end;
728: out->copy = in->copy;
729: out->destroy = in->destroy;
730: out->view = in->view;
732: PetscMalloc2(1,&out_to,1,&out_from);
733: PetscMalloc2(in_to->n,&out_to->vslots,in_from->n,&out_from->vslots);
734: out_to->n = in_to->n;
735: out_to->type = in_to->type;
736: out_to->nonmatching_computed = PETSC_FALSE;
737: out_to->slots_nonmatching = 0;
738: out_to->is_copy = PETSC_FALSE;
739: PetscMemcpy(out_to->vslots,in_to->vslots,(out_to->n)*sizeof(PetscInt));
741: out_from->n = in_from->n;
742: out_from->type = in_from->type;
743: out_from->nonmatching_computed = PETSC_FALSE;
744: out_from->slots_nonmatching = 0;
745: out_from->is_copy = PETSC_FALSE;
746: PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));
748: out->todata = (void*)out_to;
749: out->fromdata = (void*)out_from;
750: return(0);
751: }
755: PetscErrorCode VecScatterView_SGToSG(VecScatter in,PetscViewer viewer)
756: {
757: PetscErrorCode ierr;
758: VecScatter_Seq_General *in_to = (VecScatter_Seq_General*)in->todata;
759: VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata;
760: PetscInt i;
761: PetscBool isascii;
764: PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
765: if (isascii) {
766: PetscViewerASCIIPrintf(viewer,"Sequential general scatter\n");
767: for (i=0; i<in_to->n; i++) {
768: PetscViewerASCIIPrintf(viewer,"%D to %D\n",in_from->vslots[i],in_to->vslots[i]);
769: }
770: }
771: return(0);
772: }
777: PetscErrorCode VecScatterCopy_SGToSS(VecScatter in,VecScatter out)
778: {
779: PetscErrorCode ierr;
780: VecScatter_Seq_Stride *in_to = (VecScatter_Seq_Stride*)in->todata,*out_to = NULL;
781: VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from = NULL;
784: out->begin = in->begin;
785: out->end = in->end;
786: out->copy = in->copy;
787: out->destroy = in->destroy;
788: out->view = in->view;
790: PetscMalloc2(1,&out_to,1,&out_from);
791: PetscMalloc1(in_from->n,&out_from->vslots);
792: out_to->n = in_to->n;
793: out_to->type = in_to->type;
794: out_to->first = in_to->first;
795: out_to->step = in_to->step;
796: out_to->type = in_to->type;
798: out_from->n = in_from->n;
799: out_from->type = in_from->type;
800: out_from->nonmatching_computed = PETSC_FALSE;
801: out_from->slots_nonmatching = 0;
802: out_from->is_copy = PETSC_FALSE;
803: PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));
805: out->todata = (void*)out_to;
806: out->fromdata = (void*)out_from;
807: return(0);
808: }
812: PetscErrorCode VecScatterView_SGToSS(VecScatter in,PetscViewer viewer)
813: {
814: PetscErrorCode ierr;
815: VecScatter_Seq_Stride *in_to = (VecScatter_Seq_Stride*)in->todata;
816: VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata;
817: PetscInt i;
818: PetscBool isascii;
821: PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
822: if (isascii) {
823: PetscViewerASCIIPrintf(viewer,"Sequential general scatter to stride\n");
824: for (i=0; i<in_to->n; i++) {
825: PetscViewerASCIIPrintf(viewer,"%D to %D\n",in_from->vslots[i],in_to->first + in_to->step*i);
826: }
827: }
828: return(0);
829: }
831: /* --------------------------------------------------------------------------*/
832: /*
833: Scatter: parallel to sequential vector, sequential strides for both.
834: */
837: PetscErrorCode VecScatterCopy_SSToSS(VecScatter in,VecScatter out)
838: {
839: VecScatter_Seq_Stride *in_to = (VecScatter_Seq_Stride*)in->todata,*out_to = NULL;
840: VecScatter_Seq_Stride *in_from = (VecScatter_Seq_Stride*)in->fromdata,*out_from = NULL;
841: PetscErrorCode ierr;
844: out->begin = in->begin;
845: out->end = in->end;
846: out->copy = in->copy;
847: out->destroy = in->destroy;
848: out->view = in->view;
850: PetscMalloc2(1,&out_to,1,&out_from);
851: out_to->n = in_to->n;
852: out_to->type = in_to->type;
853: out_to->first = in_to->first;
854: out_to->step = in_to->step;
855: out_to->type = in_to->type;
856: out_from->n = in_from->n;
857: out_from->type = in_from->type;
858: out_from->first = in_from->first;
859: out_from->step = in_from->step;
860: out_from->type = in_from->type;
861: out->todata = (void*)out_to;
862: out->fromdata = (void*)out_from;
863: return(0);
864: }
868: PetscErrorCode VecScatterView_SSToSS(VecScatter in,PetscViewer viewer)
869: {
870: VecScatter_Seq_Stride *in_to = (VecScatter_Seq_Stride*)in->todata;
871: VecScatter_Seq_Stride *in_from = (VecScatter_Seq_Stride*)in->fromdata;
872: PetscErrorCode ierr;
873: PetscBool isascii;
876: PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
877: if (isascii) {
878: 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);
879: }
880: return(0);
881: }
884: extern PetscErrorCode VecScatterCreate_PtoS(PetscInt,const PetscInt*,PetscInt,const PetscInt*,Vec,Vec,PetscInt,VecScatter);
885: extern PetscErrorCode VecScatterCreate_PtoP(PetscInt,const PetscInt*,PetscInt,const PetscInt*,Vec,Vec,PetscInt,VecScatter);
886: extern PetscErrorCode VecScatterCreate_StoP(PetscInt,const PetscInt*,PetscInt,const PetscInt*,Vec,Vec,PetscInt,VecScatter);
888: /* =======================================================================*/
889: #define VEC_SEQ_ID 0
890: #define VEC_MPI_ID 1
891: #define IS_GENERAL_ID 0
892: #define IS_STRIDE_ID 1
893: #define IS_BLOCK_ID 2
895: /*
896: Blocksizes we have optimized scatters for
897: */
899: #define VecScatterOptimizedBS(mbs) (2 <= mbs)
901: PetscErrorCode VecScatterCreateEmpty(MPI_Comm comm,VecScatter *newctx)
902: {
903: VecScatter ctx;
907: PetscHeaderCreate(ctx,_p_VecScatter,int,VEC_SCATTER_CLASSID,"VecScatter","VecScatter","Vec",comm,VecScatterDestroy,VecScatterView);
908: ctx->inuse = PETSC_FALSE;
909: ctx->beginandendtogether = PETSC_FALSE;
910: PetscOptionsGetBool(NULL,"-vecscatter_merge",&ctx->beginandendtogether,NULL);
911: if (ctx->beginandendtogether) {
912: PetscInfo(ctx,"Using combined (merged) vector scatter begin and end\n");
913: }
914: ctx->packtogether = PETSC_FALSE;
915: PetscOptionsGetBool(NULL,"-vecscatter_packtogether",&ctx->packtogether,NULL);
916: if (ctx->packtogether) {
917: PetscInfo(ctx,"Pack all messages before sending\n");
918: }
919: *newctx = ctx;
920: return(0);
921: }
925: /*@C
926: VecScatterCreate - Creates a vector scatter context.
928: Collective on Vec
930: Input Parameters:
931: + xin - a vector that defines the shape (parallel data layout of the vector)
932: of vectors from which we scatter
933: . yin - a vector that defines the shape (parallel data layout of the vector)
934: of vectors to which we scatter
935: . ix - the indices of xin to scatter (if NULL scatters all values)
936: - iy - the indices of yin to hold results (if NULL fills entire vector yin)
938: Output Parameter:
939: . newctx - location to store the new scatter context
941: Options Database Keys: (uses regular MPI_Sends by default)
942: + -vecscatter_view - Prints detail of communications
943: . -vecscatter_view ::ascii_info - Print less details about communication
944: . -vecscatter_ssend - Uses MPI_Ssend_init() instead of MPI_Send_init()
945: . -vecscatter_rsend - use ready receiver mode for MPI sends
946: . -vecscatter_merge - VecScatterBegin() handles all of the communication, VecScatterEnd() is a nop
947: eliminates the chance for overlap of computation and communication
948: . -vecscatter_sendfirst - Posts sends before receives
949: . -vecscatter_packtogether - Pack all messages before sending, receive all messages before unpacking
950: . -vecscatter_alltoall - Uses MPI all to all communication for scatter
951: . -vecscatter_window - Use MPI 2 window operations to move data
952: . -vecscatter_nopack - Avoid packing to work vector when possible (if used with -vecscatter_alltoall then will use MPI_Alltoallw()
953: - -vecscatter_reproduce - insure that the order of the communications are done the same for each scatter, this under certain circumstances
954: will make the results of scatters deterministic when otherwise they are not (it may be slower also).
956: $
957: $ --When packing is used--
958: $ MPI Datatypes (no packing) sendfirst merge packtogether persistent*
959: $ _nopack _sendfirst _merge _packtogether -vecscatter_
960: $ ----------------------------------------------------------------------------------------------------------------------------
961: $ Message passing Send p X X X always
962: $ Ssend p X X X always _ssend
963: $ Rsend p nonsense X X always _rsend
964: $ AlltoAll v or w X nonsense always X nonsense _alltoall
965: $ MPI_Win p nonsense p p nonsense _window
966: $
967: $ Since persistent sends and receives require a constant memory address they can only be used when data is packed into the work vector
968: $ because the in and out array may be different for each call to VecScatterBegin/End().
969: $
970: $ p indicates possible, but not implemented. X indicates implemented
971: $
973: Level: intermediate
975: Notes:
976: In calls to VecScatter() you can use different vectors than the xin and
977: yin you used above; BUT they must have the same parallel data layout, for example,
978: they could be obtained from VecDuplicate().
979: A VecScatter context CANNOT be used in two or more simultaneous scatters;
980: that is you cannot call a second VecScatterBegin() with the same scatter
981: context until the VecScatterEnd() has been called on the first VecScatterBegin().
982: In this case a separate VecScatter is needed for each concurrent scatter.
984: Currently the MPI_Send(), MPI_Ssend() and MPI_Rsend() all use PERSISTENT versions.
985: (this unfortunately requires that the same in and out arrays be used for each use, this
986: is why when not using MPI_alltoallw() we always need to pack the input into the work array before sending
987: and unpack upon receeving instead of using MPI datatypes to avoid the packing/unpacking).
989: Both ix and iy cannot be NULL at the same time.
991: Concepts: scatter^between vectors
992: Concepts: gather^between vectors
994: .seealso: VecScatterDestroy(), VecScatterCreateToAll(), VecScatterCreateToZero()
995: @*/
996: PetscErrorCode VecScatterCreate(Vec xin,IS ix,Vec yin,IS iy,VecScatter *newctx)
997: {
998: VecScatter ctx;
999: PetscErrorCode ierr;
1000: PetscMPIInt size;
1001: PetscInt xin_type = VEC_SEQ_ID,yin_type = VEC_SEQ_ID,*range;
1002: PetscInt ix_type = IS_GENERAL_ID,iy_type = IS_GENERAL_ID;
1003: MPI_Comm comm,ycomm;
1004: PetscBool totalv,ixblock,iyblock,iystride,islocal,cando,flag;
1005: IS tix = 0,tiy = 0;
1008: if (!ix && !iy) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_SUP,"Cannot pass default in for both input and output indices");
1010: /*
1011: Determine if the vectors are "parallel", ie. it shares a comm with other processors, or
1012: sequential (it does not share a comm). The difference is that parallel vectors treat the
1013: index set as providing indices in the global parallel numbering of the vector, with
1014: sequential vectors treat the index set as providing indices in the local sequential
1015: numbering
1016: */
1017: PetscObjectGetComm((PetscObject)xin,&comm);
1018: MPI_Comm_size(comm,&size);
1019: if (size > 1) xin_type = VEC_MPI_ID;
1021: PetscObjectGetComm((PetscObject)yin,&ycomm);
1022: MPI_Comm_size(ycomm,&size);
1023: if (size > 1) {comm = ycomm; yin_type = VEC_MPI_ID;}
1025: /* generate the Scatter context */
1026: PetscHeaderCreate(ctx,_p_VecScatter,int,VEC_SCATTER_CLASSID,"VecScatter","VecScatter","Vec",comm,VecScatterDestroy,VecScatterView);
1027: ctx->inuse = PETSC_FALSE;
1029: ctx->beginandendtogether = PETSC_FALSE;
1030: PetscOptionsGetBool(NULL,"-vecscatter_merge",&ctx->beginandendtogether,NULL);
1031: if (ctx->beginandendtogether) {
1032: PetscInfo(ctx,"Using combined (merged) vector scatter begin and end\n");
1033: }
1034: ctx->packtogether = PETSC_FALSE;
1035: PetscOptionsGetBool(NULL,"-vecscatter_packtogether",&ctx->packtogether,NULL);
1036: if (ctx->packtogether) {
1037: PetscInfo(ctx,"Pack all messages before sending\n");
1038: }
1040: VecGetLocalSize(xin,&ctx->from_n);
1041: VecGetLocalSize(yin,&ctx->to_n);
1043: /*
1044: if ix or iy is not included; assume just grabbing entire vector
1045: */
1046: if (!ix && xin_type == VEC_SEQ_ID) {
1047: ISCreateStride(comm,ctx->from_n,0,1,&ix);
1048: tix = ix;
1049: } else if (!ix && xin_type == VEC_MPI_ID) {
1050: if (yin_type == VEC_MPI_ID) {
1051: PetscInt ntmp, low;
1052: VecGetLocalSize(xin,&ntmp);
1053: VecGetOwnershipRange(xin,&low,NULL);
1054: ISCreateStride(comm,ntmp,low,1,&ix);
1055: } else {
1056: PetscInt Ntmp;
1057: VecGetSize(xin,&Ntmp);
1058: ISCreateStride(comm,Ntmp,0,1,&ix);
1059: }
1060: tix = ix;
1061: } else if (!ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ix not given, but not Seq or MPI vector");
1063: if (!iy && yin_type == VEC_SEQ_ID) {
1064: ISCreateStride(comm,ctx->to_n,0,1,&iy);
1065: tiy = iy;
1066: } else if (!iy && yin_type == VEC_MPI_ID) {
1067: if (xin_type == VEC_MPI_ID) {
1068: PetscInt ntmp, low;
1069: VecGetLocalSize(yin,&ntmp);
1070: VecGetOwnershipRange(yin,&low,NULL);
1071: ISCreateStride(comm,ntmp,low,1,&iy);
1072: } else {
1073: PetscInt Ntmp;
1074: VecGetSize(yin,&Ntmp);
1075: ISCreateStride(comm,Ntmp,0,1,&iy);
1076: }
1077: tiy = iy;
1078: } else if (!iy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"iy not given, but not Seq or MPI vector");
1080: /*
1081: Determine types of index sets
1082: */
1083: PetscObjectTypeCompare((PetscObject)ix,ISBLOCK,&flag);
1084: if (flag) ix_type = IS_BLOCK_ID;
1085: PetscObjectTypeCompare((PetscObject)iy,ISBLOCK,&flag);
1086: if (flag) iy_type = IS_BLOCK_ID;
1087: PetscObjectTypeCompare((PetscObject)ix,ISSTRIDE,&flag);
1088: if (flag) ix_type = IS_STRIDE_ID;
1089: PetscObjectTypeCompare((PetscObject)iy,ISSTRIDE,&flag);
1090: if (flag) iy_type = IS_STRIDE_ID;
1092: /* ===========================================================================================================
1093: Check for special cases
1094: ==========================================================================================================*/
1095: /* ---------------------------------------------------------------------------*/
1096: if (xin_type == VEC_SEQ_ID && yin_type == VEC_SEQ_ID) {
1097: if (ix_type == IS_GENERAL_ID && iy_type == IS_GENERAL_ID) {
1098: PetscInt nx,ny;
1099: const PetscInt *idx,*idy;
1100: VecScatter_Seq_General *to = NULL,*from = NULL;
1102: ISGetLocalSize(ix,&nx);
1103: ISGetLocalSize(iy,&ny);
1104: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1105: ISGetIndices(ix,&idx);
1106: ISGetIndices(iy,&idy);
1107: PetscMalloc2(1,&to,1,&from);
1108: PetscMalloc2(nx,&to->vslots,nx,&from->vslots);
1109: to->n = nx;
1110: #if defined(PETSC_USE_DEBUG)
1111: VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1112: #endif
1113: PetscMemcpy(to->vslots,idy,nx*sizeof(PetscInt));
1114: from->n = nx;
1115: #if defined(PETSC_USE_DEBUG)
1116: VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1117: #endif
1118: PetscMemcpy(from->vslots,idx,nx*sizeof(PetscInt));
1119: to->type = VEC_SCATTER_SEQ_GENERAL;
1120: from->type = VEC_SCATTER_SEQ_GENERAL;
1121: ctx->todata = (void*)to;
1122: ctx->fromdata = (void*)from;
1123: ctx->begin = VecScatterBegin_SGToSG;
1124: ctx->end = 0;
1125: ctx->destroy = VecScatterDestroy_SGToSG;
1126: ctx->copy = VecScatterCopy_SGToSG;
1127: ctx->view = VecScatterView_SGToSG;
1128: PetscInfo(xin,"Special case: sequential vector general scatter\n");
1129: goto functionend;
1130: } else if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1131: PetscInt nx,ny,to_first,to_step,from_first,from_step;
1132: VecScatter_Seq_Stride *from8 = NULL,*to8 = NULL;
1134: ISGetLocalSize(ix,&nx);
1135: ISGetLocalSize(iy,&ny);
1136: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1137: ISStrideGetInfo(iy,&to_first,&to_step);
1138: ISStrideGetInfo(ix,&from_first,&from_step);
1139: PetscMalloc2(1,&to8,1,&from8);
1140: to8->n = nx;
1141: to8->first = to_first;
1142: to8->step = to_step;
1143: from8->n = nx;
1144: from8->first = from_first;
1145: from8->step = from_step;
1146: to8->type = VEC_SCATTER_SEQ_STRIDE;
1147: from8->type = VEC_SCATTER_SEQ_STRIDE;
1148: ctx->todata = (void*)to8;
1149: ctx->fromdata = (void*)from8;
1150: ctx->begin = VecScatterBegin_SSToSS;
1151: ctx->end = 0;
1152: ctx->destroy = VecScatterDestroy_SSToSS;
1153: ctx->copy = VecScatterCopy_SSToSS;
1154: ctx->view = VecScatterView_SSToSS;
1155: PetscInfo(xin,"Special case: sequential vector stride to stride\n");
1156: goto functionend;
1157: } else if (ix_type == IS_GENERAL_ID && iy_type == IS_STRIDE_ID) {
1158: PetscInt nx,ny,first,step;
1159: const PetscInt *idx;
1160: VecScatter_Seq_General *from9 = NULL;
1161: VecScatter_Seq_Stride *to9 = NULL;
1163: ISGetLocalSize(ix,&nx);
1164: ISGetIndices(ix,&idx);
1165: ISGetLocalSize(iy,&ny);
1166: ISStrideGetInfo(iy,&first,&step);
1167: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1168: PetscMalloc2(1,&to9,1,&from9);
1169: PetscMalloc1(nx,&from9->vslots);
1170: to9->n = nx;
1171: to9->first = first;
1172: to9->step = step;
1173: from9->n = nx;
1174: #if defined(PETSC_USE_DEBUG)
1175: VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1176: #endif
1177: PetscMemcpy(from9->vslots,idx,nx*sizeof(PetscInt));
1178: ctx->todata = (void*)to9; ctx->fromdata = (void*)from9;
1179: if (step == 1) ctx->begin = VecScatterBegin_SGToSS_Stride1;
1180: else ctx->begin = VecScatterBegin_SGToSS;
1181: ctx->destroy = VecScatterDestroy_SGToSS;
1182: ctx->end = 0;
1183: ctx->copy = VecScatterCopy_SGToSS;
1184: ctx->view = VecScatterView_SGToSS;
1185: to9->type = VEC_SCATTER_SEQ_STRIDE;
1186: from9->type = VEC_SCATTER_SEQ_GENERAL;
1187: PetscInfo(xin,"Special case: sequential vector general to stride\n");
1188: goto functionend;
1189: } else if (ix_type == IS_STRIDE_ID && iy_type == IS_GENERAL_ID) {
1190: PetscInt nx,ny,first,step;
1191: const PetscInt *idy;
1192: VecScatter_Seq_General *to10 = NULL;
1193: VecScatter_Seq_Stride *from10 = NULL;
1195: ISGetLocalSize(ix,&nx);
1196: ISGetIndices(iy,&idy);
1197: ISGetLocalSize(iy,&ny);
1198: ISStrideGetInfo(ix,&first,&step);
1199: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1200: PetscMalloc2(1,&to10,1,&from10);
1201: PetscMalloc1(nx,&to10->vslots);
1202: from10->n = nx;
1203: from10->first = first;
1204: from10->step = step;
1205: to10->n = nx;
1206: #if defined(PETSC_USE_DEBUG)
1207: VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1208: #endif
1209: PetscMemcpy(to10->vslots,idy,nx*sizeof(PetscInt));
1210: ctx->todata = (void*)to10;
1211: ctx->fromdata = (void*)from10;
1212: if (step == 1) ctx->begin = VecScatterBegin_SSToSG_Stride1;
1213: else ctx->begin = VecScatterBegin_SSToSG;
1214: ctx->destroy = VecScatterDestroy_SSToSG;
1215: ctx->end = 0;
1216: ctx->copy = 0;
1217: ctx->view = VecScatterView_SSToSG;
1218: to10->type = VEC_SCATTER_SEQ_GENERAL;
1219: from10->type = VEC_SCATTER_SEQ_STRIDE;
1220: PetscInfo(xin,"Special case: sequential vector stride to general\n");
1221: goto functionend;
1222: } else {
1223: PetscInt nx,ny;
1224: const PetscInt *idx,*idy;
1225: VecScatter_Seq_General *to11 = NULL,*from11 = NULL;
1226: PetscBool idnx,idny;
1228: ISGetLocalSize(ix,&nx);
1229: ISGetLocalSize(iy,&ny);
1230: if (nx != ny) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match, in %D out %D",nx,ny);
1232: ISIdentity(ix,&idnx);
1233: ISIdentity(iy,&idny);
1234: if (idnx && idny) {
1235: VecScatter_Seq_Stride *to13 = NULL,*from13 = NULL;
1236: PetscMalloc2(1,&to13,1,&from13);
1237: to13->n = nx;
1238: to13->first = 0;
1239: to13->step = 1;
1240: from13->n = nx;
1241: from13->first = 0;
1242: from13->step = 1;
1243: to13->type = VEC_SCATTER_SEQ_STRIDE;
1244: from13->type = VEC_SCATTER_SEQ_STRIDE;
1245: ctx->todata = (void*)to13;
1246: ctx->fromdata = (void*)from13;
1247: ctx->begin = VecScatterBegin_SSToSS;
1248: ctx->end = 0;
1249: ctx->destroy = VecScatterDestroy_SSToSS;
1250: ctx->copy = VecScatterCopy_SSToSS;
1251: ctx->view = VecScatterView_SSToSS;
1252: PetscInfo(xin,"Special case: sequential copy\n");
1253: goto functionend;
1254: }
1256: ISGetIndices(iy,&idy);
1257: ISGetIndices(ix,&idx);
1258: PetscMalloc2(1,&to11,1,&from11);
1259: PetscMalloc2(nx,&to11->vslots,nx,&from11->vslots);
1260: to11->n = nx;
1261: #if defined(PETSC_USE_DEBUG)
1262: VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1263: #endif
1264: PetscMemcpy(to11->vslots,idy,nx*sizeof(PetscInt));
1265: from11->n = nx;
1266: #if defined(PETSC_USE_DEBUG)
1267: VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1268: #endif
1269: PetscMemcpy(from11->vslots,idx,nx*sizeof(PetscInt));
1270: to11->type = VEC_SCATTER_SEQ_GENERAL;
1271: from11->type = VEC_SCATTER_SEQ_GENERAL;
1272: ctx->todata = (void*)to11;
1273: ctx->fromdata = (void*)from11;
1274: ctx->begin = VecScatterBegin_SGToSG;
1275: ctx->end = 0;
1276: ctx->destroy = VecScatterDestroy_SGToSG;
1277: ctx->copy = VecScatterCopy_SGToSG;
1278: ctx->view = VecScatterView_SGToSG;
1279: ISRestoreIndices(ix,&idx);
1280: ISRestoreIndices(iy,&idy);
1281: PetscInfo(xin,"Sequential vector scatter with block indices\n");
1282: goto functionend;
1283: }
1284: }
1285: /* ---------------------------------------------------------------------------*/
1286: if (xin_type == VEC_MPI_ID && yin_type == VEC_SEQ_ID) {
1288: /* ===========================================================================================================
1289: Check for special cases
1290: ==========================================================================================================*/
1291: islocal = PETSC_FALSE;
1292: /* special case extracting (subset of) local portion */
1293: if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1294: PetscInt nx,ny,to_first,to_step,from_first,from_step;
1295: PetscInt start,end,min,max;
1296: VecScatter_Seq_Stride *from12 = NULL,*to12 = NULL;
1298: VecGetOwnershipRange(xin,&start,&end);
1299: ISGetLocalSize(ix,&nx);
1300: ISStrideGetInfo(ix,&from_first,&from_step);
1301: ISGetLocalSize(iy,&ny);
1302: ISStrideGetInfo(iy,&to_first,&to_step);
1303: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1304: ISGetMinMax(ix,&min,&max);
1305: if (min >= start && max < end) islocal = PETSC_TRUE;
1306: else islocal = PETSC_FALSE;
1307: MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1308: if (cando) {
1309: PetscMalloc2(1,&to12,1,&from12);
1310: to12->n = nx;
1311: to12->first = to_first;
1312: to12->step = to_step;
1313: from12->n = nx;
1314: from12->first = from_first-start;
1315: from12->step = from_step;
1316: to12->type = VEC_SCATTER_SEQ_STRIDE;
1317: from12->type = VEC_SCATTER_SEQ_STRIDE;
1318: ctx->todata = (void*)to12;
1319: ctx->fromdata = (void*)from12;
1320: ctx->begin = VecScatterBegin_SSToSS;
1321: ctx->end = 0;
1322: ctx->destroy = VecScatterDestroy_SSToSS;
1323: ctx->copy = VecScatterCopy_SSToSS;
1324: ctx->view = VecScatterView_SSToSS;
1325: PetscInfo(xin,"Special case: processors only getting local values\n");
1326: goto functionend;
1327: }
1328: } else {
1329: MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1330: }
1332: /* test for special case of all processors getting entire vector */
1333: /* contains check that PetscMPIInt can handle the sizes needed */
1334: totalv = PETSC_FALSE;
1335: if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1336: PetscInt i,nx,ny,to_first,to_step,from_first,from_step,N;
1337: PetscMPIInt *count = NULL,*displx;
1338: VecScatter_MPI_ToAll *sto = NULL;
1340: ISGetLocalSize(ix,&nx);
1341: ISStrideGetInfo(ix,&from_first,&from_step);
1342: ISGetLocalSize(iy,&ny);
1343: ISStrideGetInfo(iy,&to_first,&to_step);
1344: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1345: VecGetSize(xin,&N);
1346: if (nx != N) totalv = PETSC_FALSE;
1347: else if (from_first == 0 && from_step == 1 && from_first == to_first && from_step == to_step) totalv = PETSC_TRUE;
1348: else totalv = PETSC_FALSE;
1349: MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1351: #if defined(PETSC_USE_64BIT_INDICES)
1352: if (cando && (yin->map->N < PETSC_MPI_INT_MAX)) {
1353: #else
1354: if (cando) {
1355: #endif
1356: MPI_Comm_size(PetscObjectComm((PetscObject)ctx),&size);
1357: PetscMalloc3(1,&sto,size,&count,size,&displx);
1358: range = xin->map->range;
1359: for (i=0; i<size; i++) {
1360: PetscMPIIntCast(range[i+1] - range[i],count+i);
1361: PetscMPIIntCast(range[i],displx+i);
1362: }
1363: sto->count = count;
1364: sto->displx = displx;
1365: sto->work1 = 0;
1366: sto->work2 = 0;
1367: sto->type = VEC_SCATTER_MPI_TOALL;
1368: ctx->todata = (void*)sto;
1369: ctx->fromdata = 0;
1370: ctx->begin = VecScatterBegin_MPI_ToAll;
1371: ctx->end = 0;
1372: ctx->destroy = VecScatterDestroy_MPI_ToAll;
1373: ctx->copy = VecScatterCopy_MPI_ToAll;
1374: ctx->view = VecScatterView_MPI_ToAll;
1375: PetscInfo(xin,"Special case: all processors get entire parallel vector\n");
1376: goto functionend;
1377: }
1378: } else {
1379: MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1380: }
1382: /* test for special case of processor 0 getting entire vector */
1383: /* contains check that PetscMPIInt can handle the sizes needed */
1384: totalv = PETSC_FALSE;
1385: if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1386: PetscInt i,nx,ny,to_first,to_step,from_first,from_step,N;
1387: PetscMPIInt rank,*count = NULL,*displx;
1388: VecScatter_MPI_ToAll *sto = NULL;
1390: PetscObjectGetComm((PetscObject)xin,&comm);
1391: MPI_Comm_rank(comm,&rank);
1392: ISGetLocalSize(ix,&nx);
1393: ISStrideGetInfo(ix,&from_first,&from_step);
1394: ISGetLocalSize(iy,&ny);
1395: ISStrideGetInfo(iy,&to_first,&to_step);
1396: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1397: if (!rank) {
1398: VecGetSize(xin,&N);
1399: if (nx != N) totalv = PETSC_FALSE;
1400: else if (from_first == 0 && from_step == 1 &&
1401: from_first == to_first && from_step == to_step) totalv = PETSC_TRUE;
1402: else totalv = PETSC_FALSE;
1403: } else {
1404: if (!nx) totalv = PETSC_TRUE;
1405: else totalv = PETSC_FALSE;
1406: }
1407: MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1409: #if defined(PETSC_USE_64BIT_INDICES)
1410: if (cando && (yin->map->N < PETSC_MPI_INT_MAX)) {
1411: #else
1412: if (cando) {
1413: #endif
1414: MPI_Comm_size(PetscObjectComm((PetscObject)ctx),&size);
1415: PetscMalloc3(1,&sto,size,&count,size,&displx);
1416: range = xin->map->range;
1417: for (i=0; i<size; i++) {
1418: PetscMPIIntCast(range[i+1] - range[i],count+i);
1419: PetscMPIIntCast(range[i],displx+i);
1420: }
1421: sto->count = count;
1422: sto->displx = displx;
1423: sto->work1 = 0;
1424: sto->work2 = 0;
1425: sto->type = VEC_SCATTER_MPI_TOONE;
1426: ctx->todata = (void*)sto;
1427: ctx->fromdata = 0;
1428: ctx->begin = VecScatterBegin_MPI_ToOne;
1429: ctx->end = 0;
1430: ctx->destroy = VecScatterDestroy_MPI_ToAll;
1431: ctx->copy = VecScatterCopy_MPI_ToAll;
1432: ctx->view = VecScatterView_MPI_ToAll;
1433: PetscInfo(xin,"Special case: processor zero gets entire parallel vector, rest get none\n");
1434: goto functionend;
1435: }
1436: } else {
1437: MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1438: }
1440: PetscObjectTypeCompare((PetscObject)ix,ISBLOCK,&ixblock);
1441: PetscObjectTypeCompare((PetscObject)iy,ISBLOCK,&iyblock);
1442: PetscObjectTypeCompare((PetscObject)iy,ISSTRIDE,&iystride);
1443: if (ixblock) {
1444: /* special case block to block */
1445: if (iyblock) {
1446: PetscInt nx,ny,bsx,bsy;
1447: const PetscInt *idx,*idy;
1448: ISGetBlockSize(iy,&bsy);
1449: ISGetBlockSize(ix,&bsx);
1450: if (bsx == bsy && VecScatterOptimizedBS(bsx)) {
1451: ISBlockGetLocalSize(ix,&nx);
1452: ISBlockGetIndices(ix,&idx);
1453: ISBlockGetLocalSize(iy,&ny);
1454: ISBlockGetIndices(iy,&idy);
1455: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1456: VecScatterCreate_PtoS(nx,idx,ny,idy,xin,yin,bsx,ctx);
1457: ISBlockRestoreIndices(ix,&idx);
1458: ISBlockRestoreIndices(iy,&idy);
1459: PetscInfo(xin,"Special case: blocked indices\n");
1460: goto functionend;
1461: }
1462: /* special case block to stride */
1463: } else if (iystride) {
1464: PetscInt ystart,ystride,ysize,bsx;
1465: ISStrideGetInfo(iy,&ystart,&ystride);
1466: ISGetLocalSize(iy,&ysize);
1467: ISGetBlockSize(ix,&bsx);
1468: /* see if stride index set is equivalent to block index set */
1469: if (VecScatterOptimizedBS(bsx) && ((ystart % bsx) == 0) && (ystride == 1) && ((ysize % bsx) == 0)) {
1470: PetscInt nx,il,*idy;
1471: const PetscInt *idx;
1472: ISBlockGetLocalSize(ix,&nx);
1473: ISBlockGetIndices(ix,&idx);
1474: if (ysize != bsx*nx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1475: PetscMalloc1(nx,&idy);
1476: if (nx) {
1477: idy[0] = ystart/bsx;
1478: for (il=1; il<nx; il++) idy[il] = idy[il-1] + 1;
1479: }
1480: VecScatterCreate_PtoS(nx,idx,nx,idy,xin,yin,bsx,ctx);
1481: PetscFree(idy);
1482: ISBlockRestoreIndices(ix,&idx);
1483: PetscInfo(xin,"Special case: blocked indices to stride\n");
1484: goto functionend;
1485: }
1486: }
1487: }
1488: /* left over general case */
1489: {
1490: PetscInt nx,ny;
1491: const PetscInt *idx,*idy;
1492: ISGetLocalSize(ix,&nx);
1493: ISGetIndices(ix,&idx);
1494: ISGetLocalSize(iy,&ny);
1495: ISGetIndices(iy,&idy);
1496: if (nx != ny) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match (%d %d)",nx,ny);
1497: VecScatterCreate_PtoS(nx,idx,ny,idy,xin,yin,1,ctx);
1498: ISRestoreIndices(ix,&idx);
1499: ISRestoreIndices(iy,&idy);
1500: PetscInfo(xin,"General case: MPI to Seq\n");
1501: goto functionend;
1502: }
1503: }
1504: /* ---------------------------------------------------------------------------*/
1505: if (xin_type == VEC_SEQ_ID && yin_type == VEC_MPI_ID) {
1506: /* ===========================================================================================================
1507: Check for special cases
1508: ==========================================================================================================*/
1509: /* special case local copy portion */
1510: islocal = PETSC_FALSE;
1511: if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1512: PetscInt nx,ny,to_first,to_step,from_step,start,end,from_first,min,max;
1513: VecScatter_Seq_Stride *from = NULL,*to = NULL;
1515: VecGetOwnershipRange(yin,&start,&end);
1516: ISGetLocalSize(ix,&nx);
1517: ISStrideGetInfo(ix,&from_first,&from_step);
1518: ISGetLocalSize(iy,&ny);
1519: ISStrideGetInfo(iy,&to_first,&to_step);
1520: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1521: ISGetMinMax(iy,&min,&max);
1522: if (min >= start && max < end) islocal = PETSC_TRUE;
1523: else islocal = PETSC_FALSE;
1524: MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)yin));
1525: if (cando) {
1526: PetscMalloc2(1,&to,1,&from);
1527: to->n = nx;
1528: to->first = to_first-start;
1529: to->step = to_step;
1530: from->n = nx;
1531: from->first = from_first;
1532: from->step = from_step;
1533: to->type = VEC_SCATTER_SEQ_STRIDE;
1534: from->type = VEC_SCATTER_SEQ_STRIDE;
1535: ctx->todata = (void*)to;
1536: ctx->fromdata = (void*)from;
1537: ctx->begin = VecScatterBegin_SSToSS;
1538: ctx->end = 0;
1539: ctx->destroy = VecScatterDestroy_SSToSS;
1540: ctx->copy = VecScatterCopy_SSToSS;
1541: ctx->view = VecScatterView_SSToSS;
1542: PetscInfo(xin,"Special case: sequential stride to MPI stride\n");
1543: goto functionend;
1544: }
1545: } else {
1546: MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)yin));
1547: }
1548: /* special case block to stride */
1549: if (ix_type == IS_BLOCK_ID && iy_type == IS_STRIDE_ID) {
1550: PetscInt ystart,ystride,ysize,bsx;
1551: ISStrideGetInfo(iy,&ystart,&ystride);
1552: ISGetLocalSize(iy,&ysize);
1553: ISGetBlockSize(ix,&bsx);
1554: /* see if stride index set is equivalent to block index set */
1555: if (VecScatterOptimizedBS(bsx) && ((ystart % bsx) == 0) && (ystride == 1) && ((ysize % bsx) == 0)) {
1556: PetscInt nx,il,*idy;
1557: const PetscInt *idx;
1558: ISBlockGetLocalSize(ix,&nx);
1559: ISBlockGetIndices(ix,&idx);
1560: if (ysize != bsx*nx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1561: PetscMalloc1(nx,&idy);
1562: if (nx) {
1563: idy[0] = ystart/bsx;
1564: for (il=1; il<nx; il++) idy[il] = idy[il-1] + 1;
1565: }
1566: VecScatterCreate_StoP(nx,idx,nx,idy,xin,yin,bsx,ctx);
1567: PetscFree(idy);
1568: ISBlockRestoreIndices(ix,&idx);
1569: PetscInfo(xin,"Special case: Blocked indices to stride\n");
1570: goto functionend;
1571: }
1572: }
1574: /* general case */
1575: {
1576: PetscInt nx,ny;
1577: const PetscInt *idx,*idy;
1578: ISGetLocalSize(ix,&nx);
1579: ISGetIndices(ix,&idx);
1580: ISGetLocalSize(iy,&ny);
1581: ISGetIndices(iy,&idy);
1582: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1583: VecScatterCreate_StoP(nx,idx,ny,idy,xin,yin,1,ctx);
1584: ISRestoreIndices(ix,&idx);
1585: ISRestoreIndices(iy,&idy);
1586: PetscInfo(xin,"General case: Seq to MPI\n");
1587: goto functionend;
1588: }
1589: }
1590: /* ---------------------------------------------------------------------------*/
1591: if (xin_type == VEC_MPI_ID && yin_type == VEC_MPI_ID) {
1592: /* no special cases for now */
1593: PetscInt nx,ny;
1594: const PetscInt *idx,*idy;
1595: ISGetLocalSize(ix,&nx);
1596: ISGetIndices(ix,&idx);
1597: ISGetLocalSize(iy,&ny);
1598: ISGetIndices(iy,&idy);
1599: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1600: VecScatterCreate_PtoP(nx,idx,ny,idy,xin,yin,1,ctx);
1601: ISRestoreIndices(ix,&idx);
1602: ISRestoreIndices(iy,&idy);
1603: PetscInfo(xin,"General case: MPI to MPI\n");
1604: goto functionend;
1605: }
1607: functionend:
1608: *newctx = ctx;
1609: ISDestroy(&tix);
1610: ISDestroy(&tiy);
1611: VecScatterViewFromOptions(ctx,NULL,"-vecscatter_view");
1612: return(0);
1613: }
1615: /* ------------------------------------------------------------------*/
1618: /*@
1619: VecScatterGetMerged - Returns true if the scatter is completed in the VecScatterBegin()
1620: and the VecScatterEnd() does nothing
1622: Not Collective
1624: Input Parameter:
1625: . ctx - scatter context created with VecScatterCreate()
1627: Output Parameter:
1628: . flg - PETSC_TRUE if the VecScatterBegin/End() are all done during the VecScatterBegin()
1630: Level: developer
1632: .seealso: VecScatterCreate(), VecScatterEnd(), VecScatterBegin()
1633: @*/
1634: PetscErrorCode VecScatterGetMerged(VecScatter ctx,PetscBool *flg)
1635: {
1638: *flg = ctx->beginandendtogether;
1639: return(0);
1640: }
1644: /*@
1645: VecScatterBegin - Begins a generalized scatter from one vector to
1646: another. Complete the scattering phase with VecScatterEnd().
1648: Neighbor-wise Collective on VecScatter and Vec
1650: Input Parameters:
1651: + inctx - scatter context generated by VecScatterCreate()
1652: . x - the vector from which we scatter
1653: . y - the vector to which we scatter
1654: . addv - either ADD_VALUES or INSERT_VALUES, with INSERT_VALUES mode any location
1655: not scattered to retains its old value; i.e. the vector is NOT first zeroed.
1656: - mode - the scattering mode, usually SCATTER_FORWARD. The available modes are:
1657: SCATTER_FORWARD or SCATTER_REVERSE
1660: Level: intermediate
1662: Options Database: See VecScatterCreate()
1664: Notes:
1665: The vectors x and y need not be the same vectors used in the call
1666: to VecScatterCreate(), but x must have the same parallel data layout
1667: as that passed in as the x to VecScatterCreate(), similarly for the y.
1668: Most likely they have been obtained from VecDuplicate().
1670: You cannot change the values in the input vector between the calls to VecScatterBegin()
1671: and VecScatterEnd().
1673: If you use SCATTER_REVERSE the two arguments x and y should be reversed, from
1674: the SCATTER_FORWARD.
1676: y[iy[i]] = x[ix[i]], for i=0,...,ni-1
1678: This scatter is far more general than the conventional
1679: scatter, since it can be a gather or a scatter or a combination,
1680: depending on the indices ix and iy. If x is a parallel vector and y
1681: is sequential, VecScatterBegin() can serve to gather values to a
1682: single processor. Similarly, if y is parallel and x sequential, the
1683: routine can scatter from one processor to many processors.
1685: Concepts: scatter^between vectors
1686: Concepts: gather^between vectors
1688: .seealso: VecScatterCreate(), VecScatterEnd()
1689: @*/
1690: PetscErrorCode VecScatterBegin(VecScatter inctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
1691: {
1693: #if defined(PETSC_USE_DEBUG)
1694: PetscInt to_n,from_n;
1695: #endif
1700: if (inctx->inuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE," Scatter ctx already in use");
1702: #if defined(PETSC_USE_DEBUG)
1703: /*
1704: Error checking to make sure these vectors match the vectors used
1705: to create the vector scatter context. -1 in the from_n and to_n indicate the
1706: vector lengths are unknown (for example with mapped scatters) and thus
1707: no error checking is performed.
1708: */
1709: if (inctx->from_n >= 0 && inctx->to_n >= 0) {
1710: VecGetLocalSize(x,&from_n);
1711: VecGetLocalSize(y,&to_n);
1712: if (mode & SCATTER_REVERSE) {
1713: 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);
1714: 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);
1715: } else {
1716: 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);
1717: 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);
1718: }
1719: }
1720: #endif
1722: inctx->inuse = PETSC_TRUE;
1723: PetscLogEventBarrierBegin(VEC_ScatterBarrier,0,0,0,0,PetscObjectComm((PetscObject)inctx));
1724: (*inctx->begin)(inctx,x,y,addv,mode);
1725: if (inctx->beginandendtogether && inctx->end) {
1726: inctx->inuse = PETSC_FALSE;
1727: (*inctx->end)(inctx,x,y,addv,mode);
1728: }
1729: PetscLogEventBarrierEnd(VEC_ScatterBarrier,0,0,0,0,PetscObjectComm((PetscObject)inctx));
1730: return(0);
1731: }
1733: /* --------------------------------------------------------------------*/
1736: /*@
1737: VecScatterEnd - Ends a generalized scatter from one vector to another. Call
1738: after first calling VecScatterBegin().
1740: Neighbor-wise Collective on VecScatter and Vec
1742: Input Parameters:
1743: + ctx - scatter context generated by VecScatterCreate()
1744: . x - the vector from which we scatter
1745: . y - the vector to which we scatter
1746: . addv - either ADD_VALUES or INSERT_VALUES.
1747: - mode - the scattering mode, usually SCATTER_FORWARD. The available modes are:
1748: SCATTER_FORWARD, SCATTER_REVERSE
1750: Level: intermediate
1752: Notes:
1753: If you use SCATTER_REVERSE the arguments x and y should be reversed, from the SCATTER_FORWARD.
1755: y[iy[i]] = x[ix[i]], for i=0,...,ni-1
1757: .seealso: VecScatterBegin(), VecScatterCreate()
1758: @*/
1759: PetscErrorCode VecScatterEnd(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
1760: {
1767: ctx->inuse = PETSC_FALSE;
1768: if (!ctx->end) return(0);
1769: if (!ctx->beginandendtogether) {
1770: PetscLogEventBegin(VEC_ScatterEnd,ctx,x,y,0);
1771: (*(ctx)->end)(ctx,x,y,addv,mode);
1772: PetscLogEventEnd(VEC_ScatterEnd,ctx,x,y,0);
1773: }
1774: return(0);
1775: }
1779: /*@C
1780: VecScatterDestroy - Destroys a scatter context created by
1781: VecScatterCreate().
1783: Collective on VecScatter
1785: Input Parameter:
1786: . ctx - the scatter context
1788: Level: intermediate
1790: .seealso: VecScatterCreate(), VecScatterCopy()
1791: @*/
1792: PetscErrorCode VecScatterDestroy(VecScatter *ctx)
1793: {
1797: if (!*ctx) return(0);
1799: if ((*ctx)->inuse) SETERRQ(((PetscObject)(*ctx))->comm,PETSC_ERR_ARG_WRONGSTATE,"Scatter context is in use");
1800: if (--((PetscObject)(*ctx))->refct > 0) {*ctx = 0; return(0);}
1802: /* if memory was published with SAWs then destroy it */
1803: PetscObjectSAWsViewOff((PetscObject)(*ctx));
1804: (*(*ctx)->destroy)(*ctx);
1805: #if defined(PETSC_HAVE_CUSP)
1806: VecScatterCUSPIndicesDestroy((PetscCUSPIndices*)&((*ctx)->spptr));
1807: #endif
1808: PetscHeaderDestroy(ctx);
1809: return(0);
1810: }
1814: /*@
1815: VecScatterCopy - Makes a copy of a scatter context.
1817: Collective on VecScatter
1819: Input Parameter:
1820: . sctx - the scatter context
1822: Output Parameter:
1823: . ctx - the context copy
1825: Level: advanced
1827: .seealso: VecScatterCreate(), VecScatterDestroy()
1828: @*/
1829: PetscErrorCode VecScatterCopy(VecScatter sctx,VecScatter *ctx)
1830: {
1836: if (!sctx->copy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot copy this type");
1837: PetscHeaderCreate(*ctx,_p_VecScatter,int,VEC_SCATTER_CLASSID,"VecScatter","VecScatter","Vec",PetscObjectComm((PetscObject)sctx),VecScatterDestroy,VecScatterView);
1838: (*ctx)->to_n = sctx->to_n;
1839: (*ctx)->from_n = sctx->from_n;
1840: (*sctx->copy)(sctx,*ctx);
1841: return(0);
1842: }
1845: /* ------------------------------------------------------------------*/
1848: /*@
1849: VecScatterView - Views a vector scatter context.
1851: Collective on VecScatter
1853: Input Parameters:
1854: + ctx - the scatter context
1855: - viewer - the viewer for displaying the context
1857: Level: intermediate
1859: @*/
1860: PetscErrorCode VecScatterView(VecScatter ctx,PetscViewer viewer)
1861: {
1866: if (!viewer) {
1867: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ctx),&viewer);
1868: }
1870: if (ctx->view) {
1871: (*ctx->view)(ctx,viewer);
1872: }
1873: return(0);
1874: }
1878: /*@C
1879: VecScatterRemap - Remaps the "from" and "to" indices in a
1880: vector scatter context. FOR EXPERTS ONLY!
1882: Collective on VecScatter
1884: Input Parameters:
1885: + scat - vector scatter context
1886: . from - remapping for "from" indices (may be NULL)
1887: - to - remapping for "to" indices (may be NULL)
1889: Level: developer
1891: Notes: In the parallel case the todata is actually the indices
1892: from which the data is TAKEN! The from stuff is where the
1893: data is finally put. This is VERY VERY confusing!
1895: In the sequential case the todata is the indices where the
1896: data is put and the fromdata is where it is taken from.
1897: This is backwards from the paralllel case! CRY! CRY! CRY!
1899: @*/
1900: PetscErrorCode VecScatterRemap(VecScatter scat,PetscInt *rto,PetscInt *rfrom)
1901: {
1902: VecScatter_Seq_General *to,*from;
1903: VecScatter_MPI_General *mto;
1904: PetscInt i;
1911: from = (VecScatter_Seq_General*)scat->fromdata;
1912: mto = (VecScatter_MPI_General*)scat->todata;
1914: if (mto->type == VEC_SCATTER_MPI_TOALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Not for to all scatter");
1916: if (rto) {
1917: if (mto->type == VEC_SCATTER_MPI_GENERAL) {
1918: /* handle off processor parts */
1919: for (i=0; i<mto->starts[mto->n]; i++) mto->indices[i] = rto[mto->indices[i]];
1921: /* handle local part */
1922: to = &mto->local;
1923: for (i=0; i<to->n; i++) to->vslots[i] = rto[to->vslots[i]];
1924: } else if (from->type == VEC_SCATTER_SEQ_GENERAL) {
1925: for (i=0; i<from->n; i++) from->vslots[i] = rto[from->vslots[i]];
1926: } else if (from->type == VEC_SCATTER_SEQ_STRIDE) {
1927: VecScatter_Seq_Stride *sto = (VecScatter_Seq_Stride*)from;
1929: /* if the remapping is the identity and stride is identity then skip remap */
1930: if (sto->step == 1 && sto->first == 0) {
1931: for (i=0; i<sto->n; i++) {
1932: if (rto[i] != i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1933: }
1934: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1935: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1936: }
1938: if (rfrom) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unable to remap the FROM in scatters yet");
1940: /*
1941: Mark then vector lengths as unknown because we do not know the
1942: lengths of the remapped vectors
1943: */
1944: scat->from_n = -1;
1945: scat->to_n = -1;
1946: return(0);
1947: }