Actual source code: vscat.c
petsc-3.3-p7 2013-05-11
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/isimpl.h> /*I "petscis.h" I*/
9: #include <petsc-private/vecimpl.h> /*I "petscvec.h" I*/
11: /* Logging support */
12: PetscClassId VEC_SCATTER_CLASSID;
14: #if defined(PETSC_USE_DEBUG)
15: /*
16: Checks if any indices are less than zero and generates an error
17: */
20: static PetscErrorCode VecScatterCheckIndices_Private(PetscInt nmax,PetscInt n,const PetscInt *idx)
21: {
22: PetscInt i;
25: for (i=0; i<n; i++) {
26: if (idx[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative index %D at %D location",idx[i],i);
27: 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);
28: }
29: return(0);
30: }
31: #endif
33: /*
34: This is special scatter code for when the entire parallel vector is copied to each processor.
36: This code was written by Cameron Cooper, Occidental College, Fall 1995,
37: will working at ANL as a SERS student.
38: */
41: PetscErrorCode VecScatterBegin_MPI_ToAll(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
42: {
44: PetscInt yy_n,xx_n;
45: PetscScalar *xv,*yv;
48: VecGetLocalSize(y,&yy_n);
49: VecGetLocalSize(x,&xx_n);
50: VecGetArray(y,&yv);
51: if (x != y) {VecGetArray(x,&xv);}
53: if (mode & SCATTER_REVERSE) {
54: PetscScalar *xvt,*xvt2;
55: VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
56: PetscInt i;
57: PetscMPIInt *disply = scat->displx;
59: if (addv == INSERT_VALUES) {
60: PetscInt rstart,rend;
61: /*
62: copy the correct part of the local vector into the local storage of
63: the MPI one. Note: this operation only makes sense if all the local
64: vectors have the same values
65: */
66: VecGetOwnershipRange(y,&rstart,&rend);
67: PetscMemcpy(yv,xv+rstart,yy_n*sizeof(PetscScalar));
68: } else {
69: MPI_Comm comm;
70: PetscMPIInt rank;
71: PetscObjectGetComm((PetscObject)y,&comm);
72: MPI_Comm_rank(comm,&rank);
73: if (scat->work1) xvt = scat->work1;
74: else {
75: PetscMalloc(xx_n*sizeof(PetscScalar),&xvt);
76: scat->work1 = xvt;
77: }
78: if (!rank) { /* I am the zeroth processor, values are accumulated here */
79: if (scat->work2) xvt2 = scat->work2;
80: else {
81: PetscMalloc(xx_n*sizeof(PetscScalar),& xvt2);
82: scat->work2 = xvt2;
83: }
84: MPI_Gatherv(yv,yy_n,MPIU_SCALAR,xvt2,scat->count,disply,MPIU_SCALAR,0,((PetscObject)ctx)->comm);
85: #if defined(PETSC_USE_COMPLEX)
86: MPI_Reduce(xv,xvt,2*xx_n,MPIU_REAL,MPIU_SUM,0,((PetscObject)ctx)->comm);
87: #else
88: MPI_Reduce(xv,xvt,xx_n,MPIU_SCALAR,MPI_SUM,0,((PetscObject)ctx)->comm);
89: #endif
90: if (addv == ADD_VALUES) {
91: for (i=0; i<xx_n; i++) {
92: xvt[i] += xvt2[i];
93: }
94: #if !defined(PETSC_USE_COMPLEX)
95: } else if (addv == MAX_VALUES) {
96: for (i=0; i<xx_n; i++) {
97: xvt[i] = PetscMax(xvt[i],xvt2[i]);
98: }
99: #endif
100: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
101: MPI_Scatterv(xvt,scat->count,disply,MPIU_SCALAR,yv,yy_n,MPIU_SCALAR,0,((PetscObject)ctx)->comm);
102: } else {
103: MPI_Gatherv(yv,yy_n,MPIU_SCALAR,0, 0,0,MPIU_SCALAR,0,((PetscObject)ctx)->comm);
104: #if defined(PETSC_USE_COMPLEX)
105: MPI_Reduce(xv,xvt,2*xx_n,MPIU_REAL,MPIU_SUM,0,((PetscObject)ctx)->comm);
106: #else
107: MPI_Reduce(xv,xvt,xx_n,MPIU_SCALAR,MPI_SUM,0,((PetscObject)ctx)->comm);
108: #endif
109: MPI_Scatterv(0,scat->count,disply,MPIU_SCALAR,yv,yy_n,MPIU_SCALAR,0,((PetscObject)ctx)->comm);
110: }
111: }
112: } else {
113: PetscScalar *yvt;
114: VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
115: PetscInt i;
116: PetscMPIInt *displx = scat->displx;
118: if (addv == INSERT_VALUES) {
119: MPI_Allgatherv(xv,xx_n,MPIU_SCALAR,yv,scat->count,displx,MPIU_SCALAR,((PetscObject)ctx)->comm);
120: } else {
121: if (scat->work1) yvt = scat->work1;
122: else {
123: PetscMalloc(yy_n*sizeof(PetscScalar),&yvt);
124: scat->work1 = yvt;
125: }
126: MPI_Allgatherv(xv,xx_n,MPIU_SCALAR,yvt,scat->count,displx,MPIU_SCALAR,((PetscObject)ctx)->comm);
127: if (addv == ADD_VALUES){
128: for (i=0; i<yy_n; i++) {
129: yv[i] += yvt[i];
130: }
131: #if !defined(PETSC_USE_COMPLEX)
132: } else if (addv == MAX_VALUES) {
133: for (i=0; i<yy_n; i++) {
134: yv[i] = PetscMax(yv[i],yvt[i]);
135: }
136: #endif
137: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
138: }
139: }
140: VecRestoreArray(y,&yv);
141: if (x != y) {VecRestoreArray(x,&xv);}
142: return(0);
143: }
145: /*
146: This is special scatter code for when the entire parallel vector is copied to processor 0.
148: */
151: PetscErrorCode VecScatterBegin_MPI_ToOne(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
152: {
154: PetscMPIInt rank;
155: PetscInt yy_n,xx_n;
156: PetscScalar *xv,*yv;
157: MPI_Comm comm;
160: VecGetLocalSize(y,&yy_n);
161: VecGetLocalSize(x,&xx_n);
162: VecGetArray(x,&xv);
163: VecGetArray(y,&yv);
165: PetscObjectGetComm((PetscObject)x,&comm);
166: MPI_Comm_rank(comm,&rank);
168: /* -------- Reverse scatter; spread from processor 0 to other processors */
169: if (mode & SCATTER_REVERSE) {
170: PetscScalar *yvt;
171: VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
172: PetscInt i;
173: PetscMPIInt *disply = scat->displx;
175: if (addv == INSERT_VALUES) {
176: MPI_Scatterv(xv,scat->count,disply,MPIU_SCALAR,yv,yy_n,MPIU_SCALAR,0,((PetscObject)ctx)->comm);
177: } else {
178: if (scat->work2) yvt = scat->work2;
179: else {
180: PetscMalloc(xx_n*sizeof(PetscScalar),&yvt);
181: scat->work2 = yvt;
182: }
183: MPI_Scatterv(xv,scat->count,disply,MPIU_SCALAR,yvt,yy_n,MPIU_SCALAR,0,((PetscObject)ctx)->comm);
184: if (addv == ADD_VALUES) {
185: for (i=0; i<yy_n; i++) {
186: yv[i] += yvt[i];
187: }
188: #if !defined(PETSC_USE_COMPLEX)
189: } else if (addv == MAX_VALUES) {
190: for (i=0; i<yy_n; i++) {
191: yv[i] = PetscMax(yv[i],yvt[i]);
192: }
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,((PetscObject)ctx)->comm);
205: } else {
206: if (!rank) {
207: if (scat->work1) yvt = scat->work1;
208: else {
209: PetscMalloc(yy_n*sizeof(PetscScalar),&yvt);
210: scat->work1 = yvt;
211: }
212: }
213: MPI_Gatherv(xv,xx_n,MPIU_SCALAR,yvt,scat->count,displx,MPIU_SCALAR,0,((PetscObject)ctx)->comm);
214: if (!rank) {
215: if (addv == ADD_VALUES) {
216: for (i=0; i<yy_n; i++) {
217: yv[i] += yvt[i];
218: }
219: #if !defined(PETSC_USE_COMPLEX)
220: } else if (addv == MAX_VALUES) {
221: for (i=0; i<yy_n; i++) {
222: yv[i] = PetscMax(yv[i],yvt[i]);
223: }
224: #endif
225: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
226: }
227: }
228: }
229: VecRestoreArray(x,&xv);
230: VecRestoreArray(y,&yv);
231: return(0);
232: }
234: /*
235: The follow to are used for both VecScatterBegin_MPI_ToAll() and VecScatterBegin_MPI_ToOne()
236: */
239: PetscErrorCode VecScatterDestroy_MPI_ToAll(VecScatter ctx)
240: {
241: VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
242: PetscErrorCode ierr;
245: PetscFree(scat->work1);
246: PetscFree(scat->work2);
247: PetscFree3(ctx->todata,scat->count,scat->displx);
248: return(0);
249: }
253: PetscErrorCode VecScatterDestroy_SGtoSG(VecScatter ctx)
254: {
255: PetscErrorCode ierr;
258: PetscFree2(((VecScatter_Seq_General*)ctx->todata)->vslots,((VecScatter_Seq_General*)ctx->fromdata)->vslots);
259: PetscFree2(ctx->todata,ctx->fromdata);
260: return(0);
261: }
265: PetscErrorCode VecScatterDestroy_SGtoSS(VecScatter ctx)
266: {
267: PetscErrorCode ierr;
270: PetscFree(((VecScatter_Seq_General*)ctx->fromdata)->vslots);
271: PetscFree2(ctx->todata,ctx->fromdata);
272: return(0);
273: }
277: PetscErrorCode VecScatterDestroy_SStoSG(VecScatter ctx)
278: {
279: PetscErrorCode ierr;
282: PetscFree(((VecScatter_Seq_General*)ctx->todata)->vslots);
283: PetscFree2(ctx->todata,ctx->fromdata);
284: return(0);
285: }
289: PetscErrorCode VecScatterDestroy_SStoSS(VecScatter ctx)
290: {
294: PetscFree2(ctx->todata,ctx->fromdata);
295: return(0);
296: }
298: /* -------------------------------------------------------------------------*/
301: PetscErrorCode VecScatterCopy_MPI_ToAll(VecScatter in,VecScatter out)
302: {
303: VecScatter_MPI_ToAll *in_to = (VecScatter_MPI_ToAll*)in->todata,*sto;
304: PetscErrorCode ierr;
305: PetscMPIInt size,*count,*displx;
306: PetscInt i;
309: out->begin = in->begin;
310: out->end = in->end;
311: out->copy = in->copy;
312: out->destroy = in->destroy;
313: out->view = in->view;
315: MPI_Comm_size(((PetscObject)out)->comm,&size);
316: PetscMalloc3(1,VecScatter_MPI_ToAll,&sto,size,PetscMPIInt,&count,size,PetscMPIInt,&displx);
317: sto->type = in_to->type;
318: sto->count = count;
319: sto->displx = displx;
320: for (i=0; i<size; i++) {
321: sto->count[i] = in_to->count[i];
322: sto->displx[i] = in_to->displx[i];
323: }
324: sto->work1 = 0;
325: sto->work2 = 0;
326: out->todata = (void*)sto;
327: out->fromdata = (void*)0;
328: return(0);
329: }
331: /* --------------------------------------------------------------------------------------*/
332: /*
333: Scatter: sequential general to sequential general
334: */
337: PetscErrorCode VecScatterBegin_SGtoSG(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
338: {
339: VecScatter_Seq_General *gen_to = (VecScatter_Seq_General*)ctx->todata;
340: VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
341: PetscErrorCode ierr;
342: PetscInt i,n = gen_from->n,*fslots,*tslots;
343: PetscScalar *xv,*yv;
344:
346: VecGetArray(x,&xv);
347: if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}
348: if (mode & SCATTER_REVERSE){
349: gen_to = (VecScatter_Seq_General*)ctx->fromdata;
350: gen_from = (VecScatter_Seq_General*)ctx->todata;
351: }
352: fslots = gen_from->vslots;
353: tslots = gen_to->vslots;
355: if (addv == INSERT_VALUES) {
356: for (i=0; i<n; i++) {yv[tslots[i]] = xv[fslots[i]];}
357: } else if (addv == ADD_VALUES) {
358: for (i=0; i<n; i++) {yv[tslots[i]] += xv[fslots[i]];}
359: #if !defined(PETSC_USE_COMPLEX)
360: } else if (addv == MAX_VALUES) {
361: for (i=0; i<n; i++) {yv[tslots[i]] = PetscMax(yv[tslots[i]],xv[fslots[i]]);}
362: #endif
363: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
364: VecRestoreArray(x,&xv);
365: if (x != y) {VecRestoreArray(y,&yv);}
366: return(0);
367: }
369: /*
370: Scatter: sequential general to sequential stride 1
371: */
374: PetscErrorCode VecScatterBegin_SGtoSS_Stride1(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
375: {
376: VecScatter_Seq_Stride *gen_to = (VecScatter_Seq_Stride*)ctx->todata;
377: VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
378: PetscInt i,n = gen_from->n,*fslots = gen_from->vslots;
379: PetscErrorCode ierr;
380: PetscInt first = gen_to->first;
381: PetscScalar *xv,*yv;
382:
384: VecGetArray(x,&xv);
385: if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}
386: if (mode & SCATTER_REVERSE){
387: xv += first;
388: if (addv == INSERT_VALUES) {
389: for (i=0; i<n; i++) {yv[fslots[i]] = xv[i];}
390: } else if (addv == ADD_VALUES) {
391: for (i=0; i<n; i++) {yv[fslots[i]] += xv[i];}
392: #if !defined(PETSC_USE_COMPLEX)
393: } else if (addv == MAX_VALUES) {
394: for (i=0; i<n; i++) {yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[i]);}
395: #endif
396: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
397: } else {
398: yv += first;
399: if (addv == INSERT_VALUES) {
400: for (i=0; i<n; i++) {yv[i] = xv[fslots[i]];}
401: } else if (addv == ADD_VALUES) {
402: for (i=0; i<n; i++) {yv[i] += xv[fslots[i]];}
403: #if !defined(PETSC_USE_COMPLEX)
404: } else if (addv == MAX_VALUES) {
405: for (i=0; i<n; i++) {yv[i] = PetscMax(yv[i],xv[fslots[i]]);}
406: #endif
407: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
408: }
409: VecRestoreArray(x,&xv);
410: if (x != y) {VecRestoreArray(y,&yv);}
411: return(0);
412: }
414: /*
415: Scatter: sequential general to sequential stride
416: */
419: PetscErrorCode VecScatterBegin_SGtoSS(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
420: {
421: VecScatter_Seq_Stride *gen_to = (VecScatter_Seq_Stride*)ctx->todata;
422: VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
423: PetscInt i,n = gen_from->n,*fslots = gen_from->vslots;
424: PetscErrorCode ierr;
425: PetscInt first = gen_to->first,step = gen_to->step;
426: PetscScalar *xv,*yv;
427:
429: VecGetArray(x,&xv);
430: if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}
432: if (mode & SCATTER_REVERSE){
433: if (addv == INSERT_VALUES) {
434: for (i=0; i<n; i++) {yv[fslots[i]] = xv[first + i*step];}
435: } else if (addv == ADD_VALUES) {
436: for (i=0; i<n; i++) {yv[fslots[i]] += xv[first + i*step];}
437: #if !defined(PETSC_USE_COMPLEX)
438: } else if (addv == MAX_VALUES) {
439: for (i=0; i<n; i++) {yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[first + i*step]);}
440: #endif
441: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
442: } else {
443: if (addv == INSERT_VALUES) {
444: for (i=0; i<n; i++) {yv[first + i*step] = xv[fslots[i]];}
445: } else if (addv == ADD_VALUES) {
446: for (i=0; i<n; i++) {yv[first + i*step] += xv[fslots[i]];}
447: #if !defined(PETSC_USE_COMPLEX)
448: } else if (addv == MAX_VALUES) {
449: for (i=0; i<n; i++) {yv[first + i*step] = PetscMax(yv[first + i*step],xv[fslots[i]]);}
450: #endif
451: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
452: }
453: VecRestoreArray(x,&xv);
454: if (x != y) {VecRestoreArray(y,&yv);}
455: return(0);
456: }
458: /*
459: Scatter: sequential stride 1 to sequential general
460: */
463: PetscErrorCode VecScatterBegin_SStoSG_Stride1(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
464: {
465: VecScatter_Seq_Stride *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
466: VecScatter_Seq_General *gen_to = (VecScatter_Seq_General*)ctx->todata;
467: PetscInt i,n = gen_from->n,*fslots = gen_to->vslots;
468: PetscErrorCode ierr;
469: PetscInt first = gen_from->first;
470: PetscScalar *xv,*yv;
471:
473: VecGetArray(x,&xv);
474: if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}
476: if (mode & SCATTER_REVERSE){
477: yv += first;
478: if (addv == INSERT_VALUES) {
479: for (i=0; i<n; i++) {yv[i] = xv[fslots[i]];}
480: } else if (addv == ADD_VALUES) {
481: for (i=0; i<n; i++) {yv[i] += xv[fslots[i]];}
482: #if !defined(PETSC_USE_COMPLEX)
483: } else if (addv == MAX_VALUES) {
484: for (i=0; i<n; i++) {yv[i] = PetscMax(yv[i],xv[fslots[i]]);}
485: #endif
486: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
487: } else {
488: xv += first;
489: if (addv == INSERT_VALUES) {
490: for (i=0; i<n; i++) {yv[fslots[i]] = xv[i];}
491: } else if (addv == ADD_VALUES) {
492: for (i=0; i<n; i++) {yv[fslots[i]] += xv[i];}
493: #if !defined(PETSC_USE_COMPLEX)
494: } else if (addv == MAX_VALUES) {
495: for (i=0; i<n; i++) {yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[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: }
506: /*
507: Scatter: sequential stride to sequential general
508: */
509: PetscErrorCode VecScatterBegin_SStoSG(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,step = gen_from->step;
516: PetscScalar *xv,*yv;
517:
519: VecGetArray(x,&xv);
520: if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}
522: if (mode & SCATTER_REVERSE){
523: if (addv == INSERT_VALUES) {
524: for (i=0; i<n; i++) {yv[first + i*step] = xv[fslots[i]];}
525: } else if (addv == ADD_VALUES) {
526: for (i=0; i<n; i++) {yv[first + i*step] += xv[fslots[i]];}
527: #if !defined(PETSC_USE_COMPLEX)
528: } else if (addv == MAX_VALUES) {
529: for (i=0; i<n; i++) {yv[first + i*step] = PetscMax(yv[first + i*step],xv[fslots[i]]);}
530: #endif
531: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
532: } else {
533: if (addv == INSERT_VALUES) {
534: for (i=0; i<n; i++) {yv[fslots[i]] = xv[first + i*step];}
535: } else if (addv == ADD_VALUES) {
536: for (i=0; i<n; i++) {yv[fslots[i]] += xv[first + i*step];}
537: #if !defined(PETSC_USE_COMPLEX)
538: } else if (addv == MAX_VALUES) {
539: for (i=0; i<n; i++) {yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[first + i*step]);}
540: #endif
541: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
542: }
543: VecRestoreArray(x,&xv);
544: if (x != y) {VecRestoreArray(y,&yv);}
545: return(0);
546: }
548: /*
549: Scatter: sequential stride to sequential stride
550: */
553: PetscErrorCode VecScatterBegin_SStoSS(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
554: {
555: VecScatter_Seq_Stride *gen_to = (VecScatter_Seq_Stride*)ctx->todata;
556: VecScatter_Seq_Stride *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
557: PetscInt i,n = gen_from->n,to_first = gen_to->first,to_step = gen_to->step;
558: PetscErrorCode ierr;
559: PetscInt from_first = gen_from->first,from_step = gen_from->step;
560: PetscScalar *xv,*yv;
561:
563: VecGetArrayRead(x,(const PetscScalar **)&xv);
564: if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}
566: if (mode & SCATTER_REVERSE){
567: from_first = gen_to->first;
568: to_first = gen_from->first;
569: from_step = gen_to->step;
570: to_step = gen_from->step;
571: }
573: if (addv == INSERT_VALUES) {
574: if (to_step == 1 && from_step == 1) {
575: PetscMemcpy(yv+to_first,xv+from_first,n*sizeof(PetscScalar));
576: } else {
577: for (i=0; i<n; i++) {
578: yv[to_first + i*to_step] = xv[from_first+i*from_step];
579: }
580: }
581: } else if (addv == ADD_VALUES) {
582: if (to_step == 1 && from_step == 1) {
583: yv += to_first; xv += from_first;
584: for (i=0; i<n; i++) {
585: yv[i] += xv[i];
586: }
587: } else {
588: for (i=0; i<n; i++) {
589: yv[to_first + i*to_step] += xv[from_first+i*from_step];
590: }
591: }
592: #if !defined(PETSC_USE_COMPLEX)
593: } else if (addv == MAX_VALUES) {
594: if (to_step == 1 && from_step == 1) {
595: yv += to_first; xv += from_first;
596: for (i=0; i<n; i++) {
597: yv[i] = PetscMax(yv[i],xv[i]);
598: }
599: } else {
600: for (i=0; i<n; i++) {
601: yv[to_first + i*to_step] = PetscMax(yv[to_first + i*to_step],xv[from_first+i*from_step]);
602: }
603: }
604: #endif
605: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
606: VecRestoreArrayRead(x,(const PetscScalar **)&xv);
607: if (x != y) {VecRestoreArray(y,&yv);}
608: return(0);
609: }
611: /* --------------------------------------------------------------------------*/
616: PetscErrorCode VecScatterCopy_SGToSG(VecScatter in,VecScatter out)
617: {
618: PetscErrorCode ierr;
619: VecScatter_Seq_General *in_to = (VecScatter_Seq_General*)in->todata,*out_to = PETSC_NULL;
620: VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from = PETSC_NULL;
621:
623: out->begin = in->begin;
624: out->end = in->end;
625: out->copy = in->copy;
626: out->destroy = in->destroy;
627: out->view = in->view;
629: PetscMalloc2(1,VecScatter_Seq_General,&out_to,1,VecScatter_Seq_General,&out_from);
630: PetscMalloc2(in_to->n,PetscInt,&out_to->vslots,in_from->n,PetscInt,&out_from->vslots);
631: out_to->n = in_to->n;
632: out_to->type = in_to->type;
633: out_to->nonmatching_computed = PETSC_FALSE;
634: out_to->slots_nonmatching = 0;
635: out_to->is_copy = PETSC_FALSE;
636: PetscMemcpy(out_to->vslots,in_to->vslots,(out_to->n)*sizeof(PetscInt));
638: out_from->n = in_from->n;
639: out_from->type = in_from->type;
640: out_from->nonmatching_computed = PETSC_FALSE;
641: out_from->slots_nonmatching = 0;
642: out_from->is_copy = PETSC_FALSE;
643: PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));
645: out->todata = (void*)out_to;
646: out->fromdata = (void*)out_from;
647: return(0);
648: }
653: PetscErrorCode VecScatterCopy_SGToSS(VecScatter in,VecScatter out)
654: {
655: PetscErrorCode ierr;
656: VecScatter_Seq_Stride *in_to = (VecScatter_Seq_Stride*)in->todata,*out_to = PETSC_NULL;
657: VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from = PETSC_NULL;
658:
660: out->begin = in->begin;
661: out->end = in->end;
662: out->copy = in->copy;
663: out->destroy = in->destroy;
664: out->view = in->view;
666: PetscMalloc2(1,VecScatter_Seq_Stride,&out_to,1,VecScatter_Seq_General,&out_from);
667: PetscMalloc(in_from->n*sizeof(PetscInt),&out_from->vslots);
668: out_to->n = in_to->n;
669: out_to->type = in_to->type;
670: out_to->first = in_to->first;
671: out_to->step = in_to->step;
672: out_to->type = in_to->type;
674: out_from->n = in_from->n;
675: out_from->type = in_from->type;
676: out_from->nonmatching_computed = PETSC_FALSE;
677: out_from->slots_nonmatching = 0;
678: out_from->is_copy = PETSC_FALSE;
679: PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));
681: out->todata = (void*)out_to;
682: out->fromdata = (void*)out_from;
683: return(0);
684: }
686: /* --------------------------------------------------------------------------*/
687: /*
688: Scatter: parallel to sequential vector, sequential strides for both.
689: */
692: PetscErrorCode VecScatterCopy_SStoSS(VecScatter in,VecScatter out)
693: {
694: VecScatter_Seq_Stride *in_to = (VecScatter_Seq_Stride*)in->todata,*out_to = PETSC_NULL;
695: VecScatter_Seq_Stride *in_from = (VecScatter_Seq_Stride*)in->fromdata,*out_from = PETSC_NULL;
696: PetscErrorCode ierr;
699: out->begin = in->begin;
700: out->end = in->end;
701: out->copy = in->copy;
702: out->destroy = in->destroy;
703: out->view = in->view;
705: PetscMalloc2(1,VecScatter_Seq_Stride,&out_to,1,VecScatter_Seq_Stride,&out_from);
706: out_to->n = in_to->n;
707: out_to->type = in_to->type;
708: out_to->first = in_to->first;
709: out_to->step = in_to->step;
710: out_to->type = in_to->type;
711: out_from->n = in_from->n;
712: out_from->type = in_from->type;
713: out_from->first = in_from->first;
714: out_from->step = in_from->step;
715: out_from->type = in_from->type;
716: out->todata = (void*)out_to;
717: out->fromdata = (void*)out_from;
718: return(0);
719: }
721: extern PetscErrorCode VecScatterCreate_PtoS(PetscInt,const PetscInt *,PetscInt,const PetscInt *,Vec,Vec,PetscInt,VecScatter);
722: extern PetscErrorCode VecScatterCreate_PtoP(PetscInt,const PetscInt *,PetscInt,const PetscInt *,Vec,Vec,PetscInt,VecScatter);
723: extern PetscErrorCode VecScatterCreate_StoP(PetscInt,const PetscInt *,PetscInt,const PetscInt *,Vec,Vec,PetscInt,VecScatter);
725: /* =======================================================================*/
726: #define VEC_SEQ_ID 0
727: #define VEC_MPI_ID 1
728: #define IS_GENERAL_ID 0
729: #define IS_STRIDE_ID 1
730: #define IS_BLOCK_ID 2
732: /*
733: Blocksizes we have optimized scatters for
734: */
736: #define VecScatterOptimizedBS(mbs) ((2 <= mbs && mbs <= 8) || mbs == 12)
738: PetscErrorCode VecScatterCreateEmpty(MPI_Comm comm,VecScatter *newctx)
739: {
740: VecScatter ctx;
744: PetscHeaderCreate(ctx,_p_VecScatter,int,VEC_SCATTER_CLASSID,0,"VecScatter","VecScatter","Vec",comm,VecScatterDestroy,VecScatterView);
745: ctx->inuse = PETSC_FALSE;
746: ctx->beginandendtogether = PETSC_FALSE;
747: PetscOptionsGetBool(PETSC_NULL,"-vecscatter_merge",&ctx->beginandendtogether,PETSC_NULL);
748: if (ctx->beginandendtogether) {
749: PetscInfo(ctx,"Using combined (merged) vector scatter begin and end\n");
750: }
751: ctx->packtogether = PETSC_FALSE;
752: PetscOptionsGetBool(PETSC_NULL,"-vecscatter_packtogether",&ctx->packtogether,PETSC_NULL);
753: if (ctx->packtogether) {
754: PetscInfo(ctx,"Pack all messages before sending\n");
755: }
756: *newctx = ctx;
757: return(0);
758: }
762: /*@C
763: VecScatterCreate - Creates a vector scatter context.
765: Collective on Vec
767: Input Parameters:
768: + xin - a vector that defines the shape (parallel data layout of the vector)
769: of vectors from which we scatter
770: . yin - a vector that defines the shape (parallel data layout of the vector)
771: of vectors to which we scatter
772: . ix - the indices of xin to scatter (if PETSC_NULL scatters all values)
773: - iy - the indices of yin to hold results (if PETSC_NULL fills entire vector yin)
775: Output Parameter:
776: . newctx - location to store the new scatter context
778: Options Database Keys: (uses regular MPI_Sends by default)
779: + -vecscatter_view - Prints detail of communications
780: . -vecscatter_view_info - Print less details about communication
781: . -vecscatter_ssend - Uses MPI_Ssend_init() instead of MPI_Send_init()
782: . -vecscatter_rsend - use ready receiver mode for MPI sends
783: . -vecscatter_merge - VecScatterBegin() handles all of the communication, VecScatterEnd() is a nop
784: eliminates the chance for overlap of computation and communication
785: . -vecscatter_sendfirst - Posts sends before receives
786: . -vecscatter_packtogether - Pack all messages before sending, receive all messages before unpacking
787: . -vecscatter_alltoall - Uses MPI all to all communication for scatter
788: . -vecscatter_window - Use MPI 2 window operations to move data
789: . -vecscatter_nopack - Avoid packing to work vector when possible (if used with -vecscatter_alltoall then will use MPI_Alltoallw()
790: - -vecscatter_reproduce - insure that the order of the communications are done the same for each scatter, this under certain circumstances
791: will make the results of scatters deterministic when otherwise they are not (it may be slower also).
793: $
794: $ --When packing is used--
795: $ MPI Datatypes (no packing) sendfirst merge packtogether persistent*
796: $ _nopack _sendfirst _merge _packtogether -vecscatter_
797: $ ----------------------------------------------------------------------------------------------------------------------------
798: $ Message passing Send p X X X always
799: $ Ssend p X X X always _ssend
800: $ Rsend p nonsense X X always _rsend
801: $ AlltoAll v or w X nonsense always X nonsense _alltoall
802: $ MPI_Win p nonsense p p nonsense _window
803: $
804: $ Since persistent sends and receives require a constant memory address they can only be used when data is packed into the work vector
805: $ because the in and out array may be different for each call to VecScatterBegin/End().
806: $
807: $ p indicates possible, but not implemented. X indicates implemented
808: $
810: Level: intermediate
812: Notes:
813: In calls to VecScatter() you can use different vectors than the xin and
814: yin you used above; BUT they must have the same parallel data layout, for example,
815: they could be obtained from VecDuplicate().
816: A VecScatter context CANNOT be used in two or more simultaneous scatters;
817: that is you cannot call a second VecScatterBegin() with the same scatter
818: context until the VecScatterEnd() has been called on the first VecScatterBegin().
819: In this case a separate VecScatter is needed for each concurrent scatter.
821: Currently the MPI_Send(), MPI_Ssend() and MPI_Rsend() all use PERSISTENT versions.
822: (this unfortunately requires that the same in and out arrays be used for each use, this
823: is why when not using MPI_alltoallw() we always need to pack the input into the work array before sending
824: and unpack upon receeving instead of using MPI datatypes to avoid the packing/unpacking).
826: Both ix and iy cannot be PETSC_NULL at the same time.
828: Concepts: scatter^between vectors
829: Concepts: gather^between vectors
831: .seealso: VecScatterDestroy(), VecScatterCreateToAll(), VecScatterCreateToZero()
832: @*/
833: PetscErrorCode VecScatterCreate(Vec xin,IS ix,Vec yin,IS iy,VecScatter *newctx)
834: {
835: VecScatter ctx;
837: PetscMPIInt size;
838: PetscInt totalv,xin_type = VEC_SEQ_ID,yin_type = VEC_SEQ_ID,*range;
839: PetscInt ix_type = IS_GENERAL_ID,iy_type = IS_GENERAL_ID;
840: MPI_Comm comm,ycomm;
841: PetscBool ixblock,iyblock,iystride,islocal,cando,flag;
842: IS tix = 0,tiy = 0;
845: if (!ix && !iy) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_SUP,"Cannot pass default in for both input and output indices");
847: /*
848: Determine if the vectors are "parallel", ie. it shares a comm with other processors, or
849: sequential (it does not share a comm). The difference is that parallel vectors treat the
850: index set as providing indices in the global parallel numbering of the vector, with
851: sequential vectors treat the index set as providing indices in the local sequential
852: numbering
853: */
854: PetscObjectGetComm((PetscObject)xin,&comm);
855: MPI_Comm_size(comm,&size);
856: if (size > 1) {xin_type = VEC_MPI_ID;}
858: PetscObjectGetComm((PetscObject)yin,&ycomm);
859: MPI_Comm_size(ycomm,&size);
860: if (size > 1) {comm = ycomm; yin_type = VEC_MPI_ID;}
862: /* generate the Scatter context */
863: PetscHeaderCreate(ctx,_p_VecScatter,int,VEC_SCATTER_CLASSID,0,"VecScatter","VecScatter","Vec",comm,VecScatterDestroy,VecScatterView);
864: ctx->inuse = PETSC_FALSE;
866: ctx->beginandendtogether = PETSC_FALSE;
867: PetscOptionsGetBool(PETSC_NULL,"-vecscatter_merge",&ctx->beginandendtogether,PETSC_NULL);
868: if (ctx->beginandendtogether) {
869: PetscInfo(ctx,"Using combined (merged) vector scatter begin and end\n");
870: }
871: ctx->packtogether = PETSC_FALSE;
872: PetscOptionsGetBool(PETSC_NULL,"-vecscatter_packtogether",&ctx->packtogether,PETSC_NULL);
873: if (ctx->packtogether) {
874: PetscInfo(ctx,"Pack all messages before sending\n");
875: }
877: VecGetLocalSize(xin,&ctx->from_n);
878: VecGetLocalSize(yin,&ctx->to_n);
880: /*
881: if ix or iy is not included; assume just grabbing entire vector
882: */
883: if (!ix && xin_type == VEC_SEQ_ID) {
884: ISCreateStride(comm,ctx->from_n,0,1,&ix);
885: tix = ix;
886: } else if (!ix && xin_type == VEC_MPI_ID) {
887: if (yin_type == VEC_MPI_ID) {
888: PetscInt ntmp, low;
889: VecGetLocalSize(xin,&ntmp);
890: VecGetOwnershipRange(xin,&low,PETSC_NULL);
891: ISCreateStride(comm,ntmp,low,1,&ix);
892: } else{
893: PetscInt Ntmp;
894: VecGetSize(xin,&Ntmp);
895: ISCreateStride(comm,Ntmp,0,1,&ix);
896: }
897: tix = ix;
898: } else if (!ix) {
899: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ix not given, but not Seq or MPI vector");
900: }
902: if (!iy && yin_type == VEC_SEQ_ID) {
903: ISCreateStride(comm,ctx->to_n,0,1,&iy);
904: tiy = iy;
905: } else if (!iy && yin_type == VEC_MPI_ID) {
906: if (xin_type == VEC_MPI_ID) {
907: PetscInt ntmp, low;
908: VecGetLocalSize(yin,&ntmp);
909: VecGetOwnershipRange(yin,&low,PETSC_NULL);
910: ISCreateStride(comm,ntmp,low,1,&iy);
911: } else{
912: PetscInt Ntmp;
913: VecGetSize(yin,&Ntmp);
914: ISCreateStride(comm,Ntmp,0,1,&iy);
915: }
916: tiy = iy;
917: } else if (!iy) {
918: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"iy not given, but not Seq or MPI vector");
919: }
921: /*
922: Determine types of index sets
923: */
924: PetscObjectTypeCompare((PetscObject)ix,ISBLOCK,&flag);
925: if (flag) ix_type = IS_BLOCK_ID;
926: PetscObjectTypeCompare((PetscObject)iy,ISBLOCK,&flag);
927: if (flag) iy_type = IS_BLOCK_ID;
928: PetscObjectTypeCompare((PetscObject)ix,ISSTRIDE,&flag);
929: if (flag) ix_type = IS_STRIDE_ID;
930: PetscObjectTypeCompare((PetscObject)iy,ISSTRIDE,&flag);
931: if (flag) iy_type = IS_STRIDE_ID;
933: /* ===========================================================================================================
934: Check for special cases
935: ==========================================================================================================*/
936: /* ---------------------------------------------------------------------------*/
937: if (xin_type == VEC_SEQ_ID && yin_type == VEC_SEQ_ID) {
938: if (ix_type == IS_GENERAL_ID && iy_type == IS_GENERAL_ID){
939: PetscInt nx,ny;
940: const PetscInt *idx,*idy;
941: VecScatter_Seq_General *to = PETSC_NULL,*from = PETSC_NULL;
943: ISGetLocalSize(ix,&nx);
944: ISGetLocalSize(iy,&ny);
945: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
946: ISGetIndices(ix,&idx);
947: ISGetIndices(iy,&idy);
948: PetscMalloc2(1,VecScatter_Seq_General,&to,1,VecScatter_Seq_General,&from);
949: PetscMalloc2(nx,PetscInt,&to->vslots,nx,PetscInt,&from->vslots);
950: to->n = nx;
951: #if defined(PETSC_USE_DEBUG)
952: VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
953: #endif
954: PetscMemcpy(to->vslots,idy,nx*sizeof(PetscInt));
955: from->n = nx;
956: #if defined(PETSC_USE_DEBUG)
957: VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
958: #endif
959: PetscMemcpy(from->vslots,idx,nx*sizeof(PetscInt));
960: to->type = VEC_SCATTER_SEQ_GENERAL;
961: from->type = VEC_SCATTER_SEQ_GENERAL;
962: ctx->todata = (void*)to;
963: ctx->fromdata = (void*)from;
964: ctx->begin = VecScatterBegin_SGtoSG;
965: ctx->end = 0;
966: ctx->destroy = VecScatterDestroy_SGtoSG;
967: ctx->copy = VecScatterCopy_SGToSG;
968: PetscInfo(xin,"Special case: sequential vector general scatter\n");
969: goto functionend;
970: } else if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID){
971: PetscInt nx,ny,to_first,to_step,from_first,from_step;
972: VecScatter_Seq_Stride *from8 = PETSC_NULL,*to8 = PETSC_NULL;
974: ISGetLocalSize(ix,&nx);
975: ISGetLocalSize(iy,&ny);
976: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
977: ISStrideGetInfo(iy,&to_first,&to_step);
978: ISStrideGetInfo(ix,&from_first,&from_step);
979: PetscMalloc2(1,VecScatter_Seq_Stride,&to8,1,VecScatter_Seq_Stride,&from8);
980: to8->n = nx;
981: to8->first = to_first;
982: to8->step = to_step;
983: from8->n = nx;
984: from8->first = from_first;
985: from8->step = from_step;
986: to8->type = VEC_SCATTER_SEQ_STRIDE;
987: from8->type = VEC_SCATTER_SEQ_STRIDE;
988: ctx->todata = (void*)to8;
989: ctx->fromdata = (void*)from8;
990: ctx->begin = VecScatterBegin_SStoSS;
991: ctx->end = 0;
992: ctx->destroy = VecScatterDestroy_SStoSS;
993: ctx->copy = VecScatterCopy_SStoSS;
994: PetscInfo(xin,"Special case: sequential vector stride to stride\n");
995: goto functionend;
996: } else if (ix_type == IS_GENERAL_ID && iy_type == IS_STRIDE_ID){
997: PetscInt nx,ny,first,step;
998: const PetscInt *idx;
999: VecScatter_Seq_General *from9 = PETSC_NULL;
1000: VecScatter_Seq_Stride *to9 = PETSC_NULL;
1002: ISGetLocalSize(ix,&nx);
1003: ISGetIndices(ix,&idx);
1004: ISGetLocalSize(iy,&ny);
1005: ISStrideGetInfo(iy,&first,&step);
1006: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1007: PetscMalloc2(1,VecScatter_Seq_Stride,&to9,1,VecScatter_Seq_General,&from9);
1008: PetscMalloc(nx*sizeof(PetscInt),&from9->vslots);
1009: to9->n = nx;
1010: to9->first = first;
1011: to9->step = step;
1012: from9->n = nx;
1013: #if defined(PETSC_USE_DEBUG)
1014: VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1015: #endif
1016: PetscMemcpy(from9->vslots,idx,nx*sizeof(PetscInt));
1017: ctx->todata = (void*)to9; ctx->fromdata = (void*)from9;
1018: if (step == 1) ctx->begin = VecScatterBegin_SGtoSS_Stride1;
1019: else ctx->begin = VecScatterBegin_SGtoSS;
1020: ctx->destroy = VecScatterDestroy_SGtoSS;
1021: ctx->end = 0;
1022: ctx->copy = VecScatterCopy_SGToSS;
1023: to9->type = VEC_SCATTER_SEQ_STRIDE;
1024: from9->type = VEC_SCATTER_SEQ_GENERAL;
1025: PetscInfo(xin,"Special case: sequential vector general to stride\n");
1026: goto functionend;
1027: } else if (ix_type == IS_STRIDE_ID && iy_type == IS_GENERAL_ID){
1028: PetscInt nx,ny,first,step;
1029: const PetscInt *idy;
1030: VecScatter_Seq_General *to10 = PETSC_NULL;
1031: VecScatter_Seq_Stride *from10 = PETSC_NULL;
1033: ISGetLocalSize(ix,&nx);
1034: ISGetIndices(iy,&idy);
1035: ISGetLocalSize(iy,&ny);
1036: ISStrideGetInfo(ix,&first,&step);
1037: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1038: PetscMalloc2(1,VecScatter_Seq_General,&to10,1,VecScatter_Seq_Stride,&from10);
1039: PetscMalloc(nx*sizeof(PetscInt),&to10->vslots);
1040: from10->n = nx;
1041: from10->first = first;
1042: from10->step = step;
1043: to10->n = nx;
1044: #if defined(PETSC_USE_DEBUG)
1045: VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1046: #endif
1047: PetscMemcpy(to10->vslots,idy,nx*sizeof(PetscInt));
1048: ctx->todata = (void*)to10;
1049: ctx->fromdata = (void*)from10;
1050: if (step == 1) ctx->begin = VecScatterBegin_SStoSG_Stride1;
1051: else ctx->begin = VecScatterBegin_SStoSG;
1052: ctx->destroy = VecScatterDestroy_SStoSG;
1053: ctx->end = 0;
1054: ctx->copy = 0;
1055: to10->type = VEC_SCATTER_SEQ_GENERAL;
1056: from10->type = VEC_SCATTER_SEQ_STRIDE;
1057: PetscInfo(xin,"Special case: sequential vector stride to general\n");
1058: goto functionend;
1059: } else {
1060: PetscInt nx,ny;
1061: const PetscInt *idx,*idy;
1062: VecScatter_Seq_General *to11 = PETSC_NULL,*from11 = PETSC_NULL;
1063: PetscBool idnx,idny;
1065: ISGetLocalSize(ix,&nx);
1066: ISGetLocalSize(iy,&ny);
1067: if (nx != ny) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match, in %D out %D",nx,ny);
1069: ISIdentity(ix,&idnx);
1070: ISIdentity(iy,&idny);
1071: if (idnx && idny) {
1072: VecScatter_Seq_Stride *to13 = PETSC_NULL,*from13 = PETSC_NULL;
1073: PetscMalloc2(1,VecScatter_Seq_Stride,&to13,1,VecScatter_Seq_Stride,&from13);
1074: to13->n = nx;
1075: to13->first = 0;
1076: to13->step = 1;
1077: from13->n = nx;
1078: from13->first = 0;
1079: from13->step = 1;
1080: to13->type = VEC_SCATTER_SEQ_STRIDE;
1081: from13->type = VEC_SCATTER_SEQ_STRIDE;
1082: ctx->todata = (void*)to13;
1083: ctx->fromdata = (void*)from13;
1084: ctx->begin = VecScatterBegin_SStoSS;
1085: ctx->end = 0;
1086: ctx->destroy = VecScatterDestroy_SStoSS;
1087: ctx->copy = VecScatterCopy_SStoSS;
1088: PetscInfo(xin,"Special case: sequential copy\n");
1089: goto functionend;
1090: }
1092: ISGetIndices(iy,&idy);
1093: ISGetIndices(ix,&idx);
1094: PetscMalloc2(1,VecScatter_Seq_General,&to11,1,VecScatter_Seq_General,&from11);
1095: PetscMalloc2(nx,PetscInt,&to11->vslots,nx,PetscInt,&from11->vslots);
1096: to11->n = nx;
1097: #if defined(PETSC_USE_DEBUG)
1098: VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1099: #endif
1100: PetscMemcpy(to11->vslots,idy,nx*sizeof(PetscInt));
1101: from11->n = nx;
1102: #if defined(PETSC_USE_DEBUG)
1103: VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1104: #endif
1105: PetscMemcpy(from11->vslots,idx,nx*sizeof(PetscInt));
1106: to11->type = VEC_SCATTER_SEQ_GENERAL;
1107: from11->type = VEC_SCATTER_SEQ_GENERAL;
1108: ctx->todata = (void*)to11;
1109: ctx->fromdata = (void*)from11;
1110: ctx->begin = VecScatterBegin_SGtoSG;
1111: ctx->end = 0;
1112: ctx->destroy = VecScatterDestroy_SGtoSG;
1113: ctx->copy = VecScatterCopy_SGToSG;
1114: ISRestoreIndices(ix,&idx);
1115: ISRestoreIndices(iy,&idy);
1116: PetscInfo(xin,"Sequential vector scatter with block indices\n");
1117: goto functionend;
1118: }
1119: }
1120: /* ---------------------------------------------------------------------------*/
1121: if (xin_type == VEC_MPI_ID && yin_type == VEC_SEQ_ID) {
1123: /* ===========================================================================================================
1124: Check for special cases
1125: ==========================================================================================================*/
1126: islocal = PETSC_FALSE;
1127: /* special case extracting (subset of) local portion */
1128: if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID){
1129: PetscInt nx,ny,to_first,to_step,from_first,from_step;
1130: PetscInt start,end;
1131: VecScatter_Seq_Stride *from12 = PETSC_NULL,*to12 = PETSC_NULL;
1133: VecGetOwnershipRange(xin,&start,&end);
1134: ISGetLocalSize(ix,&nx);
1135: ISStrideGetInfo(ix,&from_first,&from_step);
1136: ISGetLocalSize(iy,&ny);
1137: ISStrideGetInfo(iy,&to_first,&to_step);
1138: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1139: if (ix->min >= start && ix->max < end) islocal = PETSC_TRUE; else islocal = PETSC_FALSE;
1140: MPI_Allreduce(&islocal,&cando,1,MPI_INT,MPI_LAND,((PetscObject)xin)->comm);
1141: if (cando) {
1142: PetscMalloc2(1,VecScatter_Seq_Stride,&to12,1,VecScatter_Seq_Stride,&from12);
1143: to12->n = nx;
1144: to12->first = to_first;
1145: to12->step = to_step;
1146: from12->n = nx;
1147: from12->first = from_first-start;
1148: from12->step = from_step;
1149: to12->type = VEC_SCATTER_SEQ_STRIDE;
1150: from12->type = VEC_SCATTER_SEQ_STRIDE;
1151: ctx->todata = (void*)to12;
1152: ctx->fromdata = (void*)from12;
1153: ctx->begin = VecScatterBegin_SStoSS;
1154: ctx->end = 0;
1155: ctx->destroy = VecScatterDestroy_SStoSS;
1156: ctx->copy = VecScatterCopy_SStoSS;
1157: PetscInfo(xin,"Special case: processors only getting local values\n");
1158: goto functionend;
1159: }
1160: } else {
1161: MPI_Allreduce(&islocal,&cando,1,MPI_INT,MPI_LAND,((PetscObject)xin)->comm);
1162: }
1164: /* test for special case of all processors getting entire vector */
1165: /* contains check that PetscMPIInt can handle the sizes needed */
1166: totalv = 0;
1167: if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID){
1168: PetscInt i,nx,ny,to_first,to_step,from_first,from_step,N;
1169: PetscMPIInt *count = PETSC_NULL,*displx;
1170: VecScatter_MPI_ToAll *sto = PETSC_NULL;
1172: ISGetLocalSize(ix,&nx);
1173: ISStrideGetInfo(ix,&from_first,&from_step);
1174: ISGetLocalSize(iy,&ny);
1175: ISStrideGetInfo(iy,&to_first,&to_step);
1176: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1177: VecGetSize(xin,&N);
1178: if (nx != N) {
1179: totalv = 0;
1180: } else if (from_first == 0 && from_step == 1 && from_first == to_first && from_step == to_step){
1181: totalv = 1;
1182: } else totalv = 0;
1183: MPI_Allreduce(&totalv,&cando,1,MPI_INT,MPI_LAND,((PetscObject)xin)->comm);
1185: #if defined(PETSC_USE_64BIT_INDICES)
1186: if (cando && (yin->map->N < PETSC_MPI_INT_MAX)) {
1187: #else
1188: if (cando) {
1189: #endif
1190: MPI_Comm_size(((PetscObject)ctx)->comm,&size);
1191: PetscMalloc3(1,VecScatter_MPI_ToAll,&sto,size,PetscMPIInt,&count,size,PetscMPIInt,&displx);
1192: range = xin->map->range;
1193: for (i=0; i<size; i++) {
1194: count[i] = PetscMPIIntCast(range[i+1] - range[i]);
1195: displx[i] = PetscMPIIntCast(range[i]);
1196: }
1197: sto->count = count;
1198: sto->displx = displx;
1199: sto->work1 = 0;
1200: sto->work2 = 0;
1201: sto->type = VEC_SCATTER_MPI_TOALL;
1202: ctx->todata = (void*)sto;
1203: ctx->fromdata = 0;
1204: ctx->begin = VecScatterBegin_MPI_ToAll;
1205: ctx->end = 0;
1206: ctx->destroy = VecScatterDestroy_MPI_ToAll;
1207: ctx->copy = VecScatterCopy_MPI_ToAll;
1208: PetscInfo(xin,"Special case: all processors get entire parallel vector\n");
1209: goto functionend;
1210: }
1211: } else {
1212: MPI_Allreduce(&totalv,&cando,1,MPI_INT,MPI_LAND,((PetscObject)xin)->comm);
1213: }
1215: /* test for special case of processor 0 getting entire vector */
1216: /* contains check that PetscMPIInt can handle the sizes needed */
1217: totalv = 0;
1218: if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID){
1219: PetscInt i,nx,ny,to_first,to_step,from_first,from_step,N;
1220: PetscMPIInt rank,*count = PETSC_NULL,*displx;
1221: VecScatter_MPI_ToAll *sto = PETSC_NULL;
1223: PetscObjectGetComm((PetscObject)xin,&comm);
1224: MPI_Comm_rank(comm,&rank);
1225: ISGetLocalSize(ix,&nx);
1226: ISStrideGetInfo(ix,&from_first,&from_step);
1227: ISGetLocalSize(iy,&ny);
1228: ISStrideGetInfo(iy,&to_first,&to_step);
1229: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1230: if (!rank) {
1231: VecGetSize(xin,&N);
1232: if (nx != N) {
1233: totalv = 0;
1234: } else if (from_first == 0 && from_step == 1 &&
1235: from_first == to_first && from_step == to_step){
1236: totalv = 1;
1237: } else totalv = 0;
1238: } else {
1239: if (!nx) totalv = 1;
1240: else totalv = 0;
1241: }
1242: MPI_Allreduce(&totalv,&cando,1,MPI_INT,MPI_LAND,((PetscObject)xin)->comm);
1244: #if defined(PETSC_USE_64BIT_INDICES)
1245: if (cando && (yin->map->N < PETSC_MPI_INT_MAX)) {
1246: #else
1247: if (cando) {
1248: #endif
1249: MPI_Comm_size(((PetscObject)ctx)->comm,&size);
1250: PetscMalloc3(1,VecScatter_MPI_ToAll,&sto,size,PetscMPIInt,&count,size,PetscMPIInt,&displx);
1251: range = xin->map->range;
1252: for (i=0; i<size; i++) {
1253: count[i] = PetscMPIIntCast(range[i+1] - range[i]);
1254: displx[i] = PetscMPIIntCast(range[i]);
1255: }
1256: sto->count = count;
1257: sto->displx = displx;
1258: sto->work1 = 0;
1259: sto->work2 = 0;
1260: sto->type = VEC_SCATTER_MPI_TOONE;
1261: ctx->todata = (void*)sto;
1262: ctx->fromdata = 0;
1263: ctx->begin = VecScatterBegin_MPI_ToOne;
1264: ctx->end = 0;
1265: ctx->destroy = VecScatterDestroy_MPI_ToAll;
1266: ctx->copy = VecScatterCopy_MPI_ToAll;
1267: PetscInfo(xin,"Special case: processor zero gets entire parallel vector, rest get none\n");
1268: goto functionend;
1269: }
1270: } else {
1271: MPI_Allreduce(&totalv,&cando,1,MPI_INT,MPI_LAND,((PetscObject)xin)->comm);
1272: }
1274: PetscObjectTypeCompare((PetscObject)ix,ISBLOCK,&ixblock);
1275: PetscObjectTypeCompare((PetscObject)iy,ISBLOCK,&iyblock);
1276: PetscObjectTypeCompare((PetscObject)iy,ISSTRIDE,&iystride);
1277: if (ixblock) {
1278: /* special case block to block */
1279: if (iyblock) {
1280: PetscInt nx,ny,bsx,bsy;
1281: const PetscInt *idx,*idy;
1282: ISGetBlockSize(iy,&bsy);
1283: ISGetBlockSize(ix,&bsx);
1284: if (bsx == bsy && VecScatterOptimizedBS(bsx)) {
1285: ISBlockGetLocalSize(ix,&nx);
1286: ISBlockGetIndices(ix,&idx);
1287: ISBlockGetLocalSize(iy,&ny);
1288: ISBlockGetIndices(iy,&idy);
1289: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1290: VecScatterCreate_PtoS(nx,idx,ny,idy,xin,yin,bsx,ctx);
1291: ISBlockRestoreIndices(ix,&idx);
1292: ISBlockRestoreIndices(iy,&idy);
1293: PetscInfo(xin,"Special case: blocked indices\n");
1294: goto functionend;
1295: }
1296: /* special case block to stride */
1297: } else if (iystride) {
1298: PetscInt ystart,ystride,ysize,bsx;
1299: ISStrideGetInfo(iy,&ystart,&ystride);
1300: ISGetLocalSize(iy,&ysize);
1301: ISGetBlockSize(ix,&bsx);
1302: /* see if stride index set is equivalent to block index set */
1303: if (VecScatterOptimizedBS(bsx) && ((ystart % bsx) == 0) && (ystride == 1) && ((ysize % bsx) == 0)) {
1304: PetscInt nx,il,*idy;
1305: const PetscInt *idx;
1306: ISBlockGetLocalSize(ix,&nx);
1307: ISBlockGetIndices(ix,&idx);
1308: if (ysize != bsx*nx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1309: PetscMalloc(nx*sizeof(PetscInt),&idy);
1310: if (nx) {
1311: idy[0] = ystart/bsx;
1312: for (il=1; il<nx; il++) idy[il] = idy[il-1] + 1;
1313: }
1314: VecScatterCreate_PtoS(nx,idx,nx,idy,xin,yin,bsx,ctx);
1315: PetscFree(idy);
1316: ISBlockRestoreIndices(ix,&idx);
1317: PetscInfo(xin,"Special case: blocked indices to stride\n");
1318: goto functionend;
1319: }
1320: }
1321: }
1322: /* left over general case */
1323: {
1324: PetscInt nx,ny;
1325: const PetscInt *idx,*idy;
1326: ISGetLocalSize(ix,&nx);
1327: ISGetIndices(ix,&idx);
1328: ISGetLocalSize(iy,&ny);
1329: ISGetIndices(iy,&idy);
1330: if (nx != ny) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match (%d %d)",nx,ny);
1331: VecScatterCreate_PtoS(nx,idx,ny,idy,xin,yin,1,ctx);
1332: ISRestoreIndices(ix,&idx);
1333: ISRestoreIndices(iy,&idy);
1334: PetscInfo(xin,"General case: MPI to Seq\n");
1335: goto functionend;
1336: }
1337: }
1338: /* ---------------------------------------------------------------------------*/
1339: if (xin_type == VEC_SEQ_ID && yin_type == VEC_MPI_ID) {
1340: /* ===========================================================================================================
1341: Check for special cases
1342: ==========================================================================================================*/
1343: /* special case local copy portion */
1344: islocal = PETSC_FALSE;
1345: if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID){
1346: PetscInt nx,ny,to_first,to_step,from_step,start,end,from_first;
1347: VecScatter_Seq_Stride *from = PETSC_NULL,*to = PETSC_NULL;
1349: VecGetOwnershipRange(yin,&start,&end);
1350: ISGetLocalSize(ix,&nx);
1351: ISStrideGetInfo(ix,&from_first,&from_step);
1352: ISGetLocalSize(iy,&ny);
1353: ISStrideGetInfo(iy,&to_first,&to_step);
1354: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1355: if (iy->min >= start && iy->max < end) islocal = PETSC_TRUE; else islocal = PETSC_FALSE;
1356: MPI_Allreduce(&islocal,&cando,1,MPI_INT,MPI_LAND,((PetscObject)yin)->comm);
1357: if (cando) {
1358: PetscMalloc2(1,VecScatter_Seq_Stride,&to,1,VecScatter_Seq_Stride,&from);
1359: to->n = nx;
1360: to->first = to_first-start;
1361: to->step = to_step;
1362: from->n = nx;
1363: from->first = from_first;
1364: from->step = from_step;
1365: to->type = VEC_SCATTER_SEQ_STRIDE;
1366: from->type = VEC_SCATTER_SEQ_STRIDE;
1367: ctx->todata = (void*)to;
1368: ctx->fromdata = (void*)from;
1369: ctx->begin = VecScatterBegin_SStoSS;
1370: ctx->end = 0;
1371: ctx->destroy = VecScatterDestroy_SStoSS;
1372: ctx->copy = VecScatterCopy_SStoSS;
1373: PetscInfo(xin,"Special case: sequential stride to MPI stride\n");
1374: goto functionend;
1375: }
1376: } else {
1377: MPI_Allreduce(&islocal,&cando,1,MPI_INT,MPI_LAND,((PetscObject)yin)->comm);
1378: }
1379: /* special case block to stride */
1380: if (ix_type == IS_BLOCK_ID && iy_type == IS_STRIDE_ID){
1381: PetscInt ystart,ystride,ysize,bsx;
1382: ISStrideGetInfo(iy,&ystart,&ystride);
1383: ISGetLocalSize(iy,&ysize);
1384: ISGetBlockSize(ix,&bsx);
1385: /* see if stride index set is equivalent to block index set */
1386: if (VecScatterOptimizedBS(bsx) && ((ystart % bsx) == 0) && (ystride == 1) && ((ysize % bsx) == 0)) {
1387: PetscInt nx,il,*idy;
1388: const PetscInt *idx;
1389: ISBlockGetLocalSize(ix,&nx);
1390: ISBlockGetIndices(ix,&idx);
1391: if (ysize != bsx*nx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1392: PetscMalloc(nx*sizeof(PetscInt),&idy);
1393: if (nx) {
1394: idy[0] = ystart/bsx;
1395: for (il=1; il<nx; il++) idy[il] = idy[il-1] + 1;
1396: }
1397: VecScatterCreate_StoP(nx,idx,nx,idy,xin,yin,bsx,ctx);
1398: PetscFree(idy);
1399: ISBlockRestoreIndices(ix,&idx);
1400: PetscInfo(xin,"Special case: Blocked indices to stride\n");
1401: goto functionend;
1402: }
1403: }
1405: /* general case */
1406: {
1407: PetscInt nx,ny;
1408: const PetscInt *idx,*idy;
1409: ISGetLocalSize(ix,&nx);
1410: ISGetIndices(ix,&idx);
1411: ISGetLocalSize(iy,&ny);
1412: ISGetIndices(iy,&idy);
1413: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1414: VecScatterCreate_StoP(nx,idx,ny,idy,xin,yin,1,ctx);
1415: ISRestoreIndices(ix,&idx);
1416: ISRestoreIndices(iy,&idy);
1417: PetscInfo(xin,"General case: Seq to MPI\n");
1418: goto functionend;
1419: }
1420: }
1421: /* ---------------------------------------------------------------------------*/
1422: if (xin_type == VEC_MPI_ID && yin_type == VEC_MPI_ID) {
1423: /* no special cases for now */
1424: PetscInt nx,ny;
1425: const PetscInt *idx,*idy;
1426: ISGetLocalSize(ix,&nx);
1427: ISGetIndices(ix,&idx);
1428: ISGetLocalSize(iy,&ny);
1429: ISGetIndices(iy,&idy);
1430: if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1431: VecScatterCreate_PtoP(nx,idx,ny,idy,xin,yin,1,ctx);
1432: ISRestoreIndices(ix,&idx);
1433: ISRestoreIndices(iy,&idy);
1434: PetscInfo(xin,"General case: MPI to MPI\n");
1435: goto functionend;
1436: }
1438: functionend:
1439: *newctx = ctx;
1440: ISDestroy(&tix);
1441: ISDestroy(&tiy);
1442: flag = PETSC_FALSE;
1443: PetscOptionsGetBool(PETSC_NULL,"-vecscatter_view_info",&flag,PETSC_NULL);
1444: if (flag) {
1445: PetscViewer viewer;
1446: PetscViewerASCIIGetStdout(comm,&viewer);
1447: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1448: VecScatterView(ctx,viewer);
1449: PetscViewerPopFormat(viewer);
1450: }
1451: flag = PETSC_FALSE;
1452: PetscOptionsGetBool(PETSC_NULL,"-vecscatter_view",&flag,PETSC_NULL);
1453: if (flag) {
1454: PetscViewer viewer;
1455: PetscViewerASCIIGetStdout(comm,&viewer);
1456: VecScatterView(ctx,viewer);
1457: }
1458: return(0);
1459: }
1461: /* ------------------------------------------------------------------*/
1464: /*@
1465: VecScatterGetMerged - Returns true if the scatter is completed in the VecScatterBegin()
1466: and the VecScatterEnd() does nothing
1468: Not Collective
1470: Input Parameter:
1471: . ctx - scatter context created with VecScatterCreate()
1473: Output Parameter:
1474: . flg - PETSC_TRUE if the VecScatterBegin/End() are all done during the VecScatterBegin()
1476: Level: developer
1478: .seealso: VecScatterCreate(), VecScatterEnd(), VecScatterBegin()
1479: @*/
1480: PetscErrorCode VecScatterGetMerged(VecScatter ctx,PetscBool *flg)
1481: {
1484: *flg = ctx->beginandendtogether;
1485: return(0);
1486: }
1490: /*@
1491: VecScatterBegin - Begins a generalized scatter from one vector to
1492: another. Complete the scattering phase with VecScatterEnd().
1494: Neighbor-wise Collective on VecScatter and Vec
1496: Input Parameters:
1497: + inctx - scatter context generated by VecScatterCreate()
1498: . x - the vector from which we scatter
1499: . y - the vector to which we scatter
1500: . addv - either ADD_VALUES or INSERT_VALUES, with INSERT_VALUES mode any location
1501: not scattered to retains its old value; i.e. the vector is NOT first zeroed.
1502: - mode - the scattering mode, usually SCATTER_FORWARD. The available modes are:
1503: SCATTER_FORWARD or SCATTER_REVERSE
1506: Level: intermediate
1508: Options Database: See VecScatterCreate()
1510: Notes:
1511: The vectors x and y need not be the same vectors used in the call
1512: to VecScatterCreate(), but x must have the same parallel data layout
1513: as that passed in as the x to VecScatterCreate(), similarly for the y.
1514: Most likely they have been obtained from VecDuplicate().
1516: You cannot change the values in the input vector between the calls to VecScatterBegin()
1517: and VecScatterEnd().
1519: If you use SCATTER_REVERSE the two arguments x and y should be reversed, from
1520: the SCATTER_FORWARD.
1521:
1522: y[iy[i]] = x[ix[i]], for i=0,...,ni-1
1524: This scatter is far more general than the conventional
1525: scatter, since it can be a gather or a scatter or a combination,
1526: depending on the indices ix and iy. If x is a parallel vector and y
1527: is sequential, VecScatterBegin() can serve to gather values to a
1528: single processor. Similarly, if y is parallel and x sequential, the
1529: routine can scatter from one processor to many processors.
1531: Concepts: scatter^between vectors
1532: Concepts: gather^between vectors
1534: .seealso: VecScatterCreate(), VecScatterEnd()
1535: @*/
1536: PetscErrorCode VecScatterBegin(VecScatter inctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
1537: {
1539: #if defined(PETSC_USE_DEBUG)
1540: PetscInt to_n,from_n;
1541: #endif
1546: if (inctx->inuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE," Scatter ctx already in use");
1547: CHKMEMQ;
1549: #if defined(PETSC_USE_DEBUG)
1550: /*
1551: Error checking to make sure these vectors match the vectors used
1552: to create the vector scatter context. -1 in the from_n and to_n indicate the
1553: vector lengths are unknown (for example with mapped scatters) and thus
1554: no error checking is performed.
1555: */
1556: if (inctx->from_n >= 0 && inctx->to_n >= 0) {
1557: VecGetLocalSize(x,&from_n);
1558: VecGetLocalSize(y,&to_n);
1559: if (mode & SCATTER_REVERSE) {
1560: 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);
1561: 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);
1562: } else {
1563: 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);
1564: 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);
1565: }
1566: }
1567: #endif
1569: inctx->inuse = PETSC_TRUE;
1570: PetscLogEventBarrierBegin(VEC_ScatterBarrier,0,0,0,0,((PetscObject)inctx)->comm);
1571: (*inctx->begin)(inctx,x,y,addv,mode);
1572: if (inctx->beginandendtogether && inctx->end) {
1573: inctx->inuse = PETSC_FALSE;
1574: (*inctx->end)(inctx,x,y,addv,mode);
1575: }
1576: PetscLogEventBarrierEnd(VEC_ScatterBarrier,0,0,0,0,((PetscObject)inctx)->comm);
1577: CHKMEMQ;
1578: return(0);
1579: }
1581: /* --------------------------------------------------------------------*/
1584: /*@
1585: VecScatterEnd - Ends a generalized scatter from one vector to another. Call
1586: after first calling VecScatterBegin().
1588: Neighbor-wise Collective on VecScatter and Vec
1590: Input Parameters:
1591: + ctx - scatter context generated by VecScatterCreate()
1592: . x - the vector from which we scatter
1593: . y - the vector to which we scatter
1594: . addv - either ADD_VALUES or INSERT_VALUES.
1595: - mode - the scattering mode, usually SCATTER_FORWARD. The available modes are:
1596: SCATTER_FORWARD, SCATTER_REVERSE
1598: Level: intermediate
1600: Notes:
1601: If you use SCATTER_REVERSE the arguments x and y should be reversed, from the SCATTER_FORWARD.
1603: y[iy[i]] = x[ix[i]], for i=0,...,ni-1
1605: .seealso: VecScatterBegin(), VecScatterCreate()
1606: @*/
1607: PetscErrorCode VecScatterEnd(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
1608: {
1615: ctx->inuse = PETSC_FALSE;
1616: if (!ctx->end) return(0);
1617: CHKMEMQ;
1618: if (!ctx->beginandendtogether) {
1619: PetscLogEventBegin(VEC_ScatterEnd,ctx,x,y,0);
1620: (*(ctx)->end)(ctx,x,y,addv,mode);
1621: PetscLogEventEnd(VEC_ScatterEnd,ctx,x,y,0);
1622: }
1623: CHKMEMQ;
1624: return(0);
1625: }
1629: /*@C
1630: VecScatterDestroy - Destroys a scatter context created by
1631: VecScatterCreate().
1633: Collective on VecScatter
1635: Input Parameter:
1636: . ctx - the scatter context
1638: Level: intermediate
1640: .seealso: VecScatterCreate(), VecScatterCopy()
1641: @*/
1642: PetscErrorCode VecScatterDestroy(VecScatter *ctx)
1643: {
1647: if (!*ctx) return(0);
1649: if ((*ctx)->inuse) SETERRQ(((PetscObject)(*ctx))->comm,PETSC_ERR_ARG_WRONGSTATE,"Scatter context is in use");
1650: if (--((PetscObject)(*ctx))->refct > 0) {*ctx = 0; return(0);}
1652: /* if memory was published with AMS then destroy it */
1653: PetscObjectDepublish((*ctx));
1655: (*(*ctx)->destroy)(*ctx);
1656: PetscHeaderDestroy(ctx);
1657: return(0);
1658: }
1662: /*@
1663: VecScatterCopy - Makes a copy of a scatter context.
1665: Collective on VecScatter
1667: Input Parameter:
1668: . sctx - the scatter context
1670: Output Parameter:
1671: . ctx - the context copy
1673: Level: advanced
1675: .seealso: VecScatterCreate(), VecScatterDestroy()
1676: @*/
1677: PetscErrorCode VecScatterCopy(VecScatter sctx,VecScatter *ctx)
1678: {
1684: if (!sctx->copy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot copy this type");
1685: PetscHeaderCreate(*ctx,_p_VecScatter,int,VEC_SCATTER_CLASSID,0,"VecScatter","VecScatter","Vec",((PetscObject)sctx)->comm,VecScatterDestroy,VecScatterView);
1686: (*ctx)->to_n = sctx->to_n;
1687: (*ctx)->from_n = sctx->from_n;
1688: (*sctx->copy)(sctx,*ctx);
1689: return(0);
1690: }
1693: /* ------------------------------------------------------------------*/
1696: /*@
1697: VecScatterView - Views a vector scatter context.
1699: Collective on VecScatter
1701: Input Parameters:
1702: + ctx - the scatter context
1703: - viewer - the viewer for displaying the context
1705: Level: intermediate
1707: @*/
1708: PetscErrorCode VecScatterView(VecScatter ctx,PetscViewer viewer)
1709: {
1714: if (!viewer) {
1715: PetscViewerASCIIGetStdout(((PetscObject)ctx)->comm,&viewer);
1716: }
1718: if (ctx->view) {
1719: (*ctx->view)(ctx,viewer);
1720: }
1721: return(0);
1722: }
1726: /*@C
1727: VecScatterRemap - Remaps the "from" and "to" indices in a
1728: vector scatter context. FOR EXPERTS ONLY!
1730: Collective on VecScatter
1732: Input Parameters:
1733: + scat - vector scatter context
1734: . from - remapping for "from" indices (may be PETSC_NULL)
1735: - to - remapping for "to" indices (may be PETSC_NULL)
1737: Level: developer
1739: Notes: In the parallel case the todata is actually the indices
1740: from which the data is TAKEN! The from stuff is where the
1741: data is finally put. This is VERY VERY confusing!
1743: In the sequential case the todata is the indices where the
1744: data is put and the fromdata is where it is taken from.
1745: This is backwards from the paralllel case! CRY! CRY! CRY!
1747: @*/
1748: PetscErrorCode VecScatterRemap(VecScatter scat,PetscInt *rto,PetscInt *rfrom)
1749: {
1750: VecScatter_Seq_General *to,*from;
1751: VecScatter_MPI_General *mto;
1752: PetscInt i;
1759: from = (VecScatter_Seq_General *)scat->fromdata;
1760: mto = (VecScatter_MPI_General *)scat->todata;
1762: if (mto->type == VEC_SCATTER_MPI_TOALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Not for to all scatter");
1764: if (rto) {
1765: if (mto->type == VEC_SCATTER_MPI_GENERAL) {
1766: /* handle off processor parts */
1767: for (i=0; i<mto->starts[mto->n]; i++) {
1768: mto->indices[i] = rto[mto->indices[i]];
1769: }
1770: /* handle local part */
1771: to = &mto->local;
1772: for (i=0; i<to->n; i++) {
1773: to->vslots[i] = rto[to->vslots[i]];
1774: }
1775: } else if (from->type == VEC_SCATTER_SEQ_GENERAL) {
1776: for (i=0; i<from->n; i++) {
1777: from->vslots[i] = rto[from->vslots[i]];
1778: }
1779: } else if (from->type == VEC_SCATTER_SEQ_STRIDE) {
1780: VecScatter_Seq_Stride *sto = (VecScatter_Seq_Stride*)from;
1781:
1782: /* if the remapping is the identity and stride is identity then skip remap */
1783: if (sto->step == 1 && sto->first == 0) {
1784: for (i=0; i<sto->n; i++) {
1785: if (rto[i] != i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1786: }
1787: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1788: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1789: }
1791: if (rfrom) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unable to remap the FROM in scatters yet");
1793: /*
1794: Mark then vector lengths as unknown because we do not know the
1795: lengths of the remapped vectors
1796: */
1797: scat->from_n = -1;
1798: scat->to_n = -1;
1800: return(0);
1801: }