Actual source code: vscat.c

  1: #define PETSCVEC_DLL
  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 src/vec/is/isimpl.h
 9:  #include private/vecimpl.h

 11: /* Logging support */
 12: PetscCookie  VEC_SCATTER_COOKIE = 0;

 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,PetscInt *idx)
 21: {
 22:   PetscInt i;

 25:   for (i=0; i<n; i++) {
 26:     if (idx[i] < 0)     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative index %D at %D location",idx[i],i);
 27:     if (idx[i] >= nmax) SETERRQ3(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 
 35:    copied to each processor.

 37:    This code was written by Cameron Cooper, Occidental College, Fall 1995,
 38:    will working at ANL as a SERS student.
 39: */
 42: PetscErrorCode VecScatterBegin_MPI_ToAll(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
 43: {
 44: #if defined(PETSC_USE_64BIT_INDICES)
 46:   SETERRQ(PETSC_ERR_SUP,"Not implemented due to limited MPI support for extremely long gathers");
 47: #else 
 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;

 63:     if (addv == INSERT_VALUES) {
 64:       PetscInt rstart,rend;
 65:       /* 
 66:          copy the correct part of the local vector into the local storage of 
 67:          the MPI one.  Note: this operation only makes sense if all the local 
 68:          vectors have the same values
 69:       */
 70:       VecGetOwnershipRange(y,&rstart,&rend);
 71:       PetscMemcpy(yv,xv+rstart,yy_n*sizeof(PetscScalar));
 72:     } else {
 73:       MPI_Comm    comm;
 74:       PetscMPIInt rank;
 75:       PetscObjectGetComm((PetscObject)y,&comm);
 76:       MPI_Comm_rank(comm,&rank);
 77:       if (scat->work1) xvt = scat->work1;
 78:       else {
 79:         PetscMalloc(xx_n*sizeof(PetscScalar),&xvt);
 80:         scat->work1 = xvt;
 81:       }
 82:       if (!rank) { /* I am the zeroth processor, values are accumulated here */
 83:         if   (scat->work2) xvt2 = scat->work2;
 84:         else {
 85:           PetscMalloc(xx_n*sizeof(PetscScalar),& xvt2);
 86:           scat->work2 = xvt2;
 87:         }
 88:         MPI_Gatherv(yv,yy_n,MPIU_SCALAR,xvt2,scat->count,y->map.range,MPIU_SCALAR,0,ctx->comm);
 89: #if defined(PETSC_USE_COMPLEX)
 90:         MPI_Reduce(xv,xvt,2*xx_n,MPIU_REAL,MPI_SUM,0,ctx->comm);
 91: #else
 92:         MPI_Reduce(xv,xvt,xx_n,MPIU_SCALAR,MPI_SUM,0,ctx->comm);
 93: #endif
 94:         if (addv == ADD_VALUES) {
 95:           for (i=0; i<xx_n; i++) {
 96:             xvt[i] += xvt2[i];
 97:           }
 98: #if !defined(PETSC_USE_COMPLEX)
 99:         } else if (addv == MAX_VALUES) {
100:           for (i=0; i<xx_n; i++) {
101:             xvt[i] = PetscMax(xvt[i],xvt2[i]);
102:           }
103: #endif
104:         } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
105:         MPI_Scatterv(xvt,scat->count,y->map.range,MPIU_SCALAR,yv,yy_n,MPIU_SCALAR,0,ctx->comm);
106:       } else {
107:         MPI_Gatherv(yv,yy_n,MPIU_SCALAR,0, 0,0,MPIU_SCALAR,0,ctx->comm);
108: #if defined(PETSC_USE_COMPLEX)
109:         MPI_Reduce(xv,xvt,2*xx_n,MPIU_REAL,MPI_SUM,0,ctx->comm);
110: #else
111:         MPI_Reduce(xv,xvt,xx_n,MPIU_SCALAR,MPI_SUM,0,ctx->comm);
112: #endif
113:         MPI_Scatterv(0,scat->count,y->map.range,MPIU_SCALAR,yv,yy_n,MPIU_SCALAR,0,ctx->comm);
114:       }
115:     }
116:   } else {
117:     PetscScalar          *yvt;
118:     VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
119:     PetscInt             i;

121:     if (addv == INSERT_VALUES) {
122:       MPI_Allgatherv(xv,xx_n,MPIU_SCALAR,yv,scat->count,x->map.range,MPIU_SCALAR,ctx->comm);
123:     } else {
124:       if (scat->work1) yvt = scat->work1;
125:       else {
126:         PetscMalloc(yy_n*sizeof(PetscScalar),&yvt);
127:         scat->work1 = yvt;
128:       }
129:       MPI_Allgatherv(xv,xx_n,MPIU_SCALAR,yvt,scat->count,x->map.range,MPIU_SCALAR,ctx->comm);
130:       if (addv == ADD_VALUES){
131:         for (i=0; i<yy_n; i++) {
132:           yv[i] += yvt[i];
133:         }
134: #if !defined(PETSC_USE_COMPLEX)
135:       } else if (addv == MAX_VALUES) {
136:         for (i=0; i<yy_n; i++) {
137:           yv[i] = PetscMax(yv[i],yvt[i]);
138:         }
139: #endif
140:       } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
141:     }
142:   }
143:   VecRestoreArray(y,&yv);
144:   if (x != y) {VecRestoreArray(x,&xv);}
145: #endif
146:   return(0);
147: }

149: /*
150:       This is special scatter code for when the entire parallel vector is 
151:    copied to processor 0.

153: */
156: PetscErrorCode VecScatterBegin_MPI_ToOne(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
157: {
158: #if defined(PETSC_USE_64BIT_INDICES)
160:   SETERRQ(PETSC_ERR_SUP,"Not implemented due to limited MPI support for extremely long gathers");
161: #else 
163:   PetscMPIInt    rank;
164:   PetscInt       yy_n,xx_n;
165:   PetscScalar    *xv,*yv;
166:   MPI_Comm       comm;

169:   VecGetLocalSize(y,&yy_n);
170:   VecGetLocalSize(x,&xx_n);
171:   VecGetArray(x,&xv);
172:   VecGetArray(y,&yv);

174:   PetscObjectGetComm((PetscObject)x,&comm);
175:   MPI_Comm_rank(comm,&rank);

177:   /* --------  Reverse scatter; spread from processor 0 to other processors */
178:   if (mode & SCATTER_REVERSE) {
179:     PetscScalar          *yvt;
180:     VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
181:     PetscInt             i;

183:     if (addv == INSERT_VALUES) {
184:       MPI_Scatterv(xv,scat->count,y->map.range,MPIU_SCALAR,yv,yy_n,MPIU_SCALAR,0,ctx->comm);
185:     } else {
186:       if (scat->work2) yvt = scat->work2;
187:       else {
188:         PetscMalloc(xx_n*sizeof(PetscScalar),&yvt);
189:         scat->work2 = yvt;
190:       }
191:       MPI_Scatterv(xv,scat->count,y->map.range,MPIU_SCALAR,yvt,yy_n,MPIU_SCALAR,0,ctx->comm);
192:       if (addv == ADD_VALUES) {
193:         for (i=0; i<yy_n; i++) {
194:           yv[i] += yvt[i];
195:         }
196: #if !defined(PETSC_USE_COMPLEX)
197:       } else if (addv == MAX_VALUES) {
198:         for (i=0; i<yy_n; i++) {
199:           yv[i] = PetscMax(yv[i],yvt[i]);
200:         }
201: #endif
202:       } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
203:     }
204:   /* ---------  Forward scatter; gather all values onto processor 0 */
205:   } else {
206:     PetscScalar          *yvt  = 0;
207:     VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
208:     PetscInt             i;

210:     if (addv == INSERT_VALUES) {
211:       MPI_Gatherv(xv,xx_n,MPIU_SCALAR,yv,scat->count,x->map.range,MPIU_SCALAR,0,ctx->comm);
212:     } else {
213:       if (!rank) {
214:         if (scat->work1) yvt = scat->work1;
215:         else {
216:           PetscMalloc(yy_n*sizeof(PetscScalar),&yvt);
217:           scat->work1 = yvt;
218:         }
219:       }
220:       MPI_Gatherv(xv,xx_n,MPIU_SCALAR,yvt,scat->count,x->map.range,MPIU_SCALAR,0,ctx->comm);
221:       if (!rank) {
222:         if (addv == ADD_VALUES) {
223:           for (i=0; i<yy_n; i++) {
224:             yv[i] += yvt[i];
225:           }
226: #if !defined(PETSC_USE_COMPLEX)
227:         } else if (addv == MAX_VALUES) {
228:           for (i=0; i<yy_n; i++) {
229:             yv[i] = PetscMax(yv[i],yvt[i]);
230:           }
231: #endif
232:         }  else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
233:       }
234:     }
235:   }
236:   VecRestoreArray(x,&xv);
237:   VecRestoreArray(y,&yv);
238: #endif
239:   return(0);
240: }

242: /*
243:        The follow to are used for both VecScatterBegin_MPI_ToAll() and VecScatterBegin_MPI_ToOne()
244: */
247: PetscErrorCode VecScatterDestroy_MPI_ToAll(VecScatter ctx)
248: {
249:   VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
250:   PetscErrorCode       ierr;

253:   PetscFree(scat->work1);
254:   PetscFree(scat->work2);
255:   PetscFree2(ctx->todata,scat->count);
256:   PetscHeaderDestroy(ctx);
257:   return(0);
258: }

262: PetscErrorCode VecScatterDestroy_SGtoSG(VecScatter ctx)
263: {
264:   PetscErrorCode         ierr;

267:   PetscFree4(ctx->todata,((VecScatter_Seq_General*)ctx->todata)->vslots,ctx->fromdata,((VecScatter_Seq_General*)ctx->fromdata)->vslots);
268:   PetscHeaderDestroy(ctx);
269:   return(0);
270: }

274: PetscErrorCode VecScatterDestroy_SGtoSS(VecScatter ctx)
275: {
276:   PetscErrorCode         ierr;

279:   PetscFree3(ctx->todata,ctx->fromdata,((VecScatter_Seq_General*)ctx->fromdata)->vslots);
280:   PetscHeaderDestroy(ctx);
281:   return(0);
282: }

286: PetscErrorCode VecScatterDestroy_SStoSG(VecScatter ctx)
287: {
288:   PetscErrorCode         ierr;

291:   PetscFree3(ctx->todata,((VecScatter_Seq_General*)ctx->todata)->vslots,ctx->fromdata);
292:   PetscHeaderDestroy(ctx);
293:   return(0);
294: }

298: PetscErrorCode VecScatterDestroy_SStoSS(VecScatter ctx)
299: {

303:   PetscFree2(ctx->todata,ctx->fromdata);
304:   PetscHeaderDestroy(ctx);
305:   return(0);
306: }

308: /* -------------------------------------------------------------------------*/
311: PetscErrorCode VecScatterCopy_MPI_ToAll(VecScatter in,VecScatter out)
312: {
313:   VecScatter_MPI_ToAll *in_to = (VecScatter_MPI_ToAll*)in->todata,*sto;
314:   PetscErrorCode       ierr;
315:   PetscMPIInt          size;
316:   PetscInt             i;

319:   out->begin          = in->begin;
320:   out->end            = in->end;
321:   out->copy           = in->copy;
322:   out->destroy        = in->destroy;
323:   out->view           = in->view;

325:   MPI_Comm_size(out->comm,&size);
326:   PetscMalloc2(1,VecScatter_MPI_ToAll,&sto,size,PetscMPIInt,&sto->count);
327:   sto->type           = in_to->type;
328:   for (i=0; i<size; i++) {
329:     sto->count[i] = in_to->count[i];
330:   }
331:   sto->work1         = 0;
332:   sto->work2         = 0;
333:   out->todata        = (void*)sto;
334:   out->fromdata      = (void*)0;
335:   return(0);
336: }

338: /* --------------------------------------------------------------------------------------*/
339: /* 
340:         Scatter: sequential general to sequential general 
341: */
344: PetscErrorCode VecScatterBegin_SGtoSG(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
345: {
346:   VecScatter_Seq_General *gen_to = (VecScatter_Seq_General*)ctx->todata;
347:   VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
348:   PetscErrorCode         ierr;
349:   PetscInt               i,n = gen_from->n,*fslots,*tslots;
350:   PetscScalar            *xv,*yv;
351: 
353:   VecGetArray(x,&xv);
354:   if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}
355:   if (mode & SCATTER_REVERSE){
356:     gen_to   = (VecScatter_Seq_General*)ctx->fromdata;
357:     gen_from = (VecScatter_Seq_General*)ctx->todata;
358:   }
359:   fslots = gen_from->vslots;
360:   tslots = gen_to->vslots;

362:   if (addv == INSERT_VALUES) {
363:     for (i=0; i<n; i++) {yv[tslots[i]] = xv[fslots[i]];}
364:   } else if (addv == ADD_VALUES) {
365:     for (i=0; i<n; i++) {yv[tslots[i]] += xv[fslots[i]];}
366: #if !defined(PETSC_USE_COMPLEX)
367:   } else  if (addv == MAX_VALUES) {
368:     for (i=0; i<n; i++) {yv[tslots[i]] = PetscMax(yv[tslots[i]],xv[fslots[i]]);}
369: #endif
370:   } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
371:   VecRestoreArray(x,&xv);
372:   if (x != y) {VecRestoreArray(y,&yv);}
373:   return(0);
374: }

376: /* 
377:     Scatter: sequential general to sequential stride 1 
378: */
381: PetscErrorCode VecScatterBegin_SGtoSS_Stride1(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
382: {
383:   VecScatter_Seq_Stride  *gen_to   = (VecScatter_Seq_Stride*)ctx->todata;
384:   VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
385:   PetscInt               i,n = gen_from->n,*fslots = gen_from->vslots;
386:   PetscErrorCode         ierr;
387:   PetscInt               first = gen_to->first;
388:   PetscScalar            *xv,*yv;
389: 
391:   VecGetArray(x,&xv);
392:   if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}
393:   if (mode & SCATTER_REVERSE){
394:     xv += first;
395:     if (addv == INSERT_VALUES) {
396:       for (i=0; i<n; i++) {yv[fslots[i]] = xv[i];}
397:     } else  if (addv == ADD_VALUES) {
398:       for (i=0; i<n; i++) {yv[fslots[i]] += xv[i];}
399: #if !defined(PETSC_USE_COMPLEX)
400:     } else  if (addv == MAX_VALUES) {
401:       for (i=0; i<n; i++) {yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[i]);}
402: #endif
403:     } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
404:   } else {
405:     yv += first;
406:     if (addv == INSERT_VALUES) {
407:       for (i=0; i<n; i++) {yv[i] = xv[fslots[i]];}
408:     } else  if (addv == ADD_VALUES) {
409:       for (i=0; i<n; i++) {yv[i] += xv[fslots[i]];}
410: #if !defined(PETSC_USE_COMPLEX)
411:     } else if (addv == MAX_VALUES) {
412:       for (i=0; i<n; i++) {yv[i] = PetscMax(yv[i],xv[fslots[i]]);}
413: #endif
414:     } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
415:   }
416:   VecRestoreArray(x,&xv);
417:   if (x != y) {VecRestoreArray(y,&yv);}
418:   return(0);
419: }

421: /* 
422:    Scatter: sequential general to sequential stride 
423: */
426: PetscErrorCode VecScatterBegin_SGtoSS(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
427: {
428:   VecScatter_Seq_Stride  *gen_to   = (VecScatter_Seq_Stride*)ctx->todata;
429:   VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
430:   PetscInt               i,n = gen_from->n,*fslots = gen_from->vslots;
431:   PetscErrorCode         ierr;
432:   PetscInt               first = gen_to->first,step = gen_to->step;
433:   PetscScalar            *xv,*yv;
434: 
436:   VecGetArray(x,&xv);
437:   if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}

439:   if (mode & SCATTER_REVERSE){
440:     if (addv == INSERT_VALUES) {
441:       for (i=0; i<n; i++) {yv[fslots[i]] = xv[first + i*step];}
442:     } else if (addv == ADD_VALUES) {
443:       for (i=0; i<n; i++) {yv[fslots[i]] += xv[first + i*step];}
444: #if !defined(PETSC_USE_COMPLEX)
445:     } else if (addv == MAX_VALUES) {
446:       for (i=0; i<n; i++) {yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[first + i*step]);}
447: #endif
448:     } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
449:   } else {
450:     if (addv == INSERT_VALUES) {
451:       for (i=0; i<n; i++) {yv[first + i*step] = xv[fslots[i]];}
452:     } else if (addv == ADD_VALUES) {
453:       for (i=0; i<n; i++) {yv[first + i*step] += xv[fslots[i]];}
454: #if !defined(PETSC_USE_COMPLEX)
455:     } else if (addv == MAX_VALUES) {
456:       for (i=0; i<n; i++) {yv[first + i*step] = PetscMax(yv[first + i*step],xv[fslots[i]]);}
457: #endif
458:     } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
459:   }
460:   VecRestoreArray(x,&xv);
461:   if (x != y) {VecRestoreArray(y,&yv);}
462:   return(0);
463: }

465: /* 
466:     Scatter: sequential stride 1 to sequential general 
467: */
470: PetscErrorCode VecScatterBegin_SStoSG_Stride1(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
471: {
472:   VecScatter_Seq_Stride  *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
473:   VecScatter_Seq_General *gen_to   = (VecScatter_Seq_General*)ctx->todata;
474:   PetscInt               i,n = gen_from->n,*fslots = gen_to->vslots;
475:   PetscErrorCode         ierr;
476:   PetscInt               first = gen_from->first;
477:   PetscScalar            *xv,*yv;
478: 
480:   VecGetArray(x,&xv);
481:   if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}

483:   if (mode & SCATTER_REVERSE){
484:     yv += first;
485:     if (addv == INSERT_VALUES) {
486:       for (i=0; i<n; i++) {yv[i] = xv[fslots[i]];}
487:     } else  if (addv == ADD_VALUES) {
488:       for (i=0; i<n; i++) {yv[i] += xv[fslots[i]];}
489: #if !defined(PETSC_USE_COMPLEX)
490:     } else  if (addv == MAX_VALUES) {
491:       for (i=0; i<n; i++) {yv[i] = PetscMax(yv[i],xv[fslots[i]]);}
492: #endif
493:     } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
494:   } else {
495:     xv += first;
496:     if (addv == INSERT_VALUES) {
497:       for (i=0; i<n; i++) {yv[fslots[i]] = xv[i];}
498:     } else  if (addv == ADD_VALUES) {
499:       for (i=0; i<n; i++) {yv[fslots[i]] += xv[i];}
500: #if !defined(PETSC_USE_COMPLEX)
501:     } else  if (addv == MAX_VALUES) {
502:       for (i=0; i<n; i++) {yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[i]);}
503: #endif
504:     } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
505:   }
506:   VecRestoreArray(x,&xv);
507:   if (x != y) {VecRestoreArray(y,&yv);}
508:   return(0);
509: }

513: /* 
514:    Scatter: sequential stride to sequential general 
515: */
516: PetscErrorCode VecScatterBegin_SStoSG(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
517: {
518:   VecScatter_Seq_Stride  *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
519:   VecScatter_Seq_General *gen_to   = (VecScatter_Seq_General*)ctx->todata;
520:   PetscInt               i,n = gen_from->n,*fslots = gen_to->vslots;
521:   PetscErrorCode         ierr;
522:   PetscInt               first = gen_from->first,step = gen_from->step;
523:   PetscScalar            *xv,*yv;
524: 
526:   VecGetArray(x,&xv);
527:   if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}

529:   if (mode & SCATTER_REVERSE){
530:     if (addv == INSERT_VALUES) {
531:       for (i=0; i<n; i++) {yv[first + i*step] = xv[fslots[i]];}
532:     } else  if (addv == ADD_VALUES) {
533:       for (i=0; i<n; i++) {yv[first + i*step] += xv[fslots[i]];}
534: #if !defined(PETSC_USE_COMPLEX)
535:     } else  if (addv == MAX_VALUES) {
536:       for (i=0; i<n; i++) {yv[first + i*step] = PetscMax(yv[first + i*step],xv[fslots[i]]);}
537: #endif
538:     } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
539:   } else {
540:     if (addv == INSERT_VALUES) {
541:       for (i=0; i<n; i++) {yv[fslots[i]] = xv[first + i*step];}
542:     } else  if (addv == ADD_VALUES) {
543:       for (i=0; i<n; i++) {yv[fslots[i]] += xv[first + i*step];}
544: #if !defined(PETSC_USE_COMPLEX)
545:     } else  if (addv == MAX_VALUES) {
546:       for (i=0; i<n; i++) {yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[first + i*step]);}
547: #endif
548:     } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
549:   }
550:   VecRestoreArray(x,&xv);
551:   if (x != y) {VecRestoreArray(y,&yv);}
552:   return(0);
553: }

555: /* 
556:      Scatter: sequential stride to sequential stride 
557: */
560: PetscErrorCode VecScatterBegin_SStoSS(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
561: {
562:   VecScatter_Seq_Stride *gen_to   = (VecScatter_Seq_Stride*)ctx->todata;
563:   VecScatter_Seq_Stride *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
564:   PetscInt              i,n = gen_from->n,to_first = gen_to->first,to_step = gen_to->step;
565:   PetscErrorCode        ierr;
566:   PetscInt              from_first = gen_from->first,from_step = gen_from->step;
567:   PetscScalar           *xv,*yv;
568: 
570:   VecGetArray(x,&xv);
571:   if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}

573:   if (mode & SCATTER_REVERSE){
574:     from_first = gen_to->first;
575:     to_first   = gen_from->first;
576:     from_step  = gen_to->step;
577:     to_step    = gen_from->step;
578:   }

580:   if (addv == INSERT_VALUES) {
581:     if (to_step == 1 && from_step == 1) {
582:       PetscMemcpy(yv+to_first,xv+from_first,n*sizeof(PetscScalar));
583:     } else  {
584:       for (i=0; i<n; i++) {
585:         yv[to_first + i*to_step] = xv[from_first+i*from_step];
586:       }
587:     }
588:   } else if (addv == ADD_VALUES) {
589:     if (to_step == 1 && from_step == 1) {
590:       yv += to_first; xv += from_first;
591:       for (i=0; i<n; i++) {
592:         yv[i] += xv[i];
593:       }
594:     } else {
595:       for (i=0; i<n; i++) {
596:         yv[to_first + i*to_step] += xv[from_first+i*from_step];
597:       }
598:     }
599: #if !defined(PETSC_USE_COMPLEX)
600:   } else if (addv == MAX_VALUES) {
601:     if (to_step == 1 && from_step == 1) {
602:       yv += to_first; xv += from_first;
603:       for (i=0; i<n; i++) {
604:         yv[i] = PetscMax(yv[i],xv[i]);
605:       }
606:     } else {
607:       for (i=0; i<n; i++) {
608:         yv[to_first + i*to_step] = PetscMax(yv[to_first + i*to_step],xv[from_first+i*from_step]);
609:       }
610:     }
611: #endif
612:   } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
613:   VecRestoreArray(x,&xv);
614:   if (x != y) {VecRestoreArray(y,&yv);}
615:   return(0);
616: }

618: /* --------------------------------------------------------------------------*/


623: PetscErrorCode VecScatterCopy_SGToSG(VecScatter in,VecScatter out)
624: {
625:   PetscErrorCode         ierr;
626:   VecScatter_Seq_General *in_to   = (VecScatter_Seq_General*)in->todata,*out_to = PETSC_NULL;
627:   VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from = PETSC_NULL;
628: 
630:   out->begin         = in->begin;
631:   out->end           = in->end;
632:   out->copy          = in->copy;
633:   out->destroy       = in->destroy;
634:   out->view          = in->view;

636:   PetscMalloc4(1,VecScatter_Seq_General,&out_to,in_to->n,PetscInt,&out_to->vslots,1,VecScatter_Seq_General,&out_from,in_from->n,PetscInt,&out_from->vslots);
637:   out_to->n                      = in_to->n;
638:   out_to->type                   = in_to->type;
639:   out_to->nonmatching_computed   = PETSC_FALSE;
640:   out_to->slots_nonmatching      = 0;
641:   out_to->is_copy                = PETSC_FALSE;
642:   PetscMemcpy(out_to->vslots,in_to->vslots,(out_to->n)*sizeof(PetscInt));

644:   out_from->n                    = in_from->n;
645:   out_from->type                 = in_from->type;
646:   out_from->nonmatching_computed = PETSC_FALSE;
647:   out_from->slots_nonmatching    = 0;
648:   out_from->is_copy              = PETSC_FALSE;
649:   PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));

651:   out->todata     = (void*)out_to;
652:   out->fromdata   = (void*)out_from;
653:   return(0);
654: }


659: PetscErrorCode VecScatterCopy_SGToSS(VecScatter in,VecScatter out)
660: {
661:   PetscErrorCode         ierr;
662:   VecScatter_Seq_Stride  *in_to   = (VecScatter_Seq_Stride*)in->todata,*out_to = PETSC_NULL;
663:   VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from = PETSC_NULL;
664: 
666:   out->begin         = in->begin;
667:   out->end           = in->end;
668:   out->copy          = in->copy;
669:   out->destroy       = in->destroy;
670:   out->view          = in->view;

672:   PetscMalloc3(1,VecScatter_Seq_Stride,&out_to,1,VecScatter_Seq_General,&out_from,in_from->n,PetscInt,&out_from->vslots);
673:   out_to->n       = in_to->n;
674:   out_to->type    = in_to->type;
675:   out_to->first   = in_to->first;
676:   out_to->step    = in_to->step;
677:   out_to->type    = in_to->type;

679:   out_from->n                    = in_from->n;
680:   out_from->type                 = in_from->type;
681:   out_from->nonmatching_computed = PETSC_FALSE;
682:   out_from->slots_nonmatching    = 0;
683:   out_from->is_copy              = PETSC_FALSE;
684:   PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));

686:   out->todata     = (void*)out_to;
687:   out->fromdata   = (void*)out_from;
688:   return(0);
689: }

691: /* --------------------------------------------------------------------------*/
692: /* 
693:     Scatter: parallel to sequential vector, sequential strides for both. 
694: */
697: PetscErrorCode VecScatterCopy_SStoSS(VecScatter in,VecScatter out)
698: {
699:   VecScatter_Seq_Stride *in_to   = (VecScatter_Seq_Stride*)in->todata,*out_to = PETSC_NULL;
700:   VecScatter_Seq_Stride *in_from = (VecScatter_Seq_Stride*)in->fromdata,*out_from = PETSC_NULL;
701:   PetscErrorCode        ierr;

704:   out->begin      = in->begin;
705:   out->end        = in->end;
706:   out->copy       = in->copy;
707:   out->destroy    = in->destroy;
708:   out->view       = in->view;

710:   PetscMalloc2(1,VecScatter_Seq_Stride,&out_to,1,VecScatter_Seq_Stride,&out_from);
711:   out_to->n       = in_to->n;
712:   out_to->type    = in_to->type;
713:   out_to->first   = in_to->first;
714:   out_to->step    = in_to->step;
715:   out_to->type    = in_to->type;
716:   out_from->n     = in_from->n;
717:   out_from->type  = in_from->type;
718:   out_from->first = in_from->first;
719:   out_from->step  = in_from->step;
720:   out_from->type  = in_from->type;
721:   out->todata     = (void*)out_to;
722:   out->fromdata   = (void*)out_from;
723:   return(0);
724: }

726: EXTERN PetscErrorCode VecScatterCreate_PtoS(PetscInt,PetscInt *,PetscInt,PetscInt *,Vec,Vec,PetscInt,VecScatter);
727: EXTERN PetscErrorCode VecScatterCreate_PtoP(PetscInt,PetscInt *,PetscInt,PetscInt *,Vec,Vec,VecScatter);
728: EXTERN PetscErrorCode VecScatterCreate_StoP(PetscInt,PetscInt *,PetscInt,PetscInt *,Vec,Vec,PetscInt,VecScatter);

730: /* =======================================================================*/
731: #define VEC_SEQ_ID 0
732: #define VEC_MPI_ID 1

734: /*
735:    Blocksizes we have optimized scatters for 
736: */

738: #define VecScatterOptimizedBS(mbs) ((2 <= mbs && mbs <= 8) || mbs == 12)

740: PetscErrorCode  VecScatterCreateEmpty(MPI_Comm comm,VecScatter *newctx)
741: {
742:   VecScatter     ctx;

746:   PetscHeaderCreate(ctx,_p_VecScatter,int,VEC_SCATTER_COOKIE,0,"VecScatter",comm,VecScatterDestroy,VecScatterView);
747:   ctx->inuse               = PETSC_FALSE;
748:   ctx->beginandendtogether = PETSC_FALSE;
749:   PetscOptionsHasName(PETSC_NULL,"-vecscatter_merge",&ctx->beginandendtogether);
750:   if (ctx->beginandendtogether) {
751:     PetscInfo(ctx,"Using combined (merged) vector scatter begin and end\n");
752:   }
753:   PetscOptionsHasName(PETSC_NULL,"-vecscatter_packtogether",&ctx->packtogether);
754:   if (ctx->packtogether) {
755:     PetscInfo(ctx,"Pack all messages before sending\n");
756:   }
757:   *newctx = ctx;
758:   return(0);
759: }

763: /*@C
764:    VecScatterCreate - Creates a vector scatter context.

766:    Collective on Vec

768:    Input Parameters:
769: +  xin - a vector that defines the shape (parallel data layout of the vector)
770:          of vectors from which we scatter
771: .  yin - a vector that defines the shape (parallel data layout of the vector)
772:          of vectors to which we scatter
773: .  ix - the indices of xin to scatter (if PETSC_NULL scatters all values)
774: -  iy - the indices of yin to hold results (if PETSC_NULL fills entire vector yin)

776:    Output Parameter:
777: .  newctx - location to store the new scatter context

779:    Options Database Keys:
780: +  -vecscatter_merge        - VecScatterBegin() handles all of the communication, VecScatterEnd() is a nop 
781:                               eliminates the chance for overlap of computation and communication 
782: .  -vecscatter_ssend        - Uses MPI_Ssend_init() instead of MPI_Send_init() 
783: .  -vecscatter_sendfirst    - Posts sends before receives 
784: .  -vecscatter_rr           - use ready receiver mode for MPI sends 
785: -  -vecscatter_packtogether - Pack all messages before sending, receive all messages before unpacking
786:                               ONLY implemented for BLOCK SIZE of 4 and 12! (others easily added)

788:     Level: intermediate

790:   Notes:
791:    In calls to VecScatter() you can use different vectors than the xin and 
792:    yin you used above; BUT they must have the same parallel data layout, for example,
793:    they could be obtained from VecDuplicate().
794:    A VecScatter context CANNOT be used in two or more simultaneous scatters;
795:    that is you cannot call a second VecScatterBegin() with the same scatter
796:    context until the VecScatterEnd() has been called on the first VecScatterBegin().
797:    In this case a separate VecScatter is needed for each concurrent scatter.

799:    Concepts: scatter^between vectors
800:    Concepts: gather^between vectors

802: .seealso: VecScatterDestroy(), VecScatterCreateToAll(), VecScatterCreateToZero()
803: @*/
804: PetscErrorCode  VecScatterCreate(Vec xin,IS ix,Vec yin,IS iy,VecScatter *newctx)
805: {
806:   VecScatter     ctx;
808:   PetscMPIInt    size;
809:   PetscInt       totalv,xin_type = VEC_SEQ_ID,yin_type = VEC_SEQ_ID,*range;
810:   MPI_Comm       comm,ycomm;
811:   PetscTruth     ixblock,iyblock,iystride,islocal,cando,flag;
812:   IS             tix = 0,tiy = 0;


816:   /*
817:       Determine if the vectors are "parallel", ie. it shares a comm with other processors, or
818:       sequential (it does not share a comm). The difference is that parallel vectors treat the 
819:       index set as providing indices in the global parallel numbering of the vector, with 
820:       sequential vectors treat the index set as providing indices in the local sequential
821:       numbering
822:   */
823:   PetscObjectGetComm((PetscObject)xin,&comm);
824:   MPI_Comm_size(comm,&size);
825:   if (size > 1) {xin_type = VEC_MPI_ID;}

827:   PetscObjectGetComm((PetscObject)yin,&ycomm);
828:   MPI_Comm_size(ycomm,&size);
829:   if (size > 1) {comm = ycomm; yin_type = VEC_MPI_ID;}
830: 
831:   /* generate the Scatter context */
832:   PetscHeaderCreate(ctx,_p_VecScatter,int,VEC_SCATTER_COOKIE,0,"VecScatter",comm,VecScatterDestroy,VecScatterView);
833:   ctx->inuse               = PETSC_FALSE;

835:   ctx->beginandendtogether = PETSC_FALSE;
836:   PetscOptionsHasName(PETSC_NULL,"-vecscatter_merge",&ctx->beginandendtogether);
837:   if (ctx->beginandendtogether) {
838:     PetscInfo(ctx,"Using combined (merged) vector scatter begin and end\n");
839:   }
840:   PetscOptionsHasName(PETSC_NULL,"-vecscatter_packtogether",&ctx->packtogether);
841:   if (ctx->packtogether) {
842:     PetscInfo(ctx,"Pack all messages before sending\n");
843:   }

845:   VecGetLocalSize(xin,&ctx->from_n);
846:   VecGetLocalSize(yin,&ctx->to_n);

848:   /*
849:       if ix or iy is not included; assume just grabbing entire vector
850:   */
851:   if (!ix && xin_type == VEC_SEQ_ID) {
852:     ISCreateStride(comm,ctx->from_n,0,1,&ix);
853:     tix  = ix;
854:   } else if (!ix && xin_type == VEC_MPI_ID) {
855:     if (yin_type == VEC_MPI_ID) {
856:       PetscInt ntmp, low;
857:       VecGetLocalSize(xin,&ntmp);
858:       VecGetOwnershipRange(xin,&low,PETSC_NULL);
859:       ISCreateStride(comm,ntmp,low,1,&ix);
860:     } else{
861:       PetscInt Ntmp;
862:       VecGetSize(xin,&Ntmp);
863:       ISCreateStride(comm,Ntmp,0,1,&ix);
864:     }
865:     tix  = ix;
866:   } else if (!ix) {
867:     SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"ix not given, but not Seq or MPI vector");
868:   }

870:   if (!iy && yin_type == VEC_SEQ_ID) {
871:     ISCreateStride(comm,ctx->to_n,0,1,&iy);
872:     tiy  = iy;
873:   } else if (!iy && yin_type == VEC_MPI_ID) {
874:     if (xin_type == VEC_MPI_ID) {
875:       PetscInt ntmp, low;
876:       VecGetLocalSize(yin,&ntmp);
877:       VecGetOwnershipRange(yin,&low,PETSC_NULL);
878:       ISCreateStride(comm,ntmp,low,1,&iy);
879:     } else{
880:       PetscInt Ntmp;
881:       VecGetSize(yin,&Ntmp);
882:       ISCreateStride(comm,Ntmp,0,1,&iy);
883:     }
884:     tiy  = iy;
885:   } else if (!iy) {
886:     SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"iy not given, but not Seq or MPI vector");
887:   }

889:   /* ===========================================================================================================
890:         Check for special cases
891:      ==========================================================================================================*/
892:   /* ---------------------------------------------------------------------------*/
893:   if (xin_type == VEC_SEQ_ID && yin_type == VEC_SEQ_ID) {
894:     if (ix->type == IS_GENERAL && iy->type == IS_GENERAL){
895:       PetscInt               nx,ny,*idx,*idy;
896:       VecScatter_Seq_General *to = PETSC_NULL,*from = PETSC_NULL;

898:       ISGetLocalSize(ix,&nx);
899:       ISGetLocalSize(iy,&ny);
900:       if (nx != ny) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
901:       ISGetIndices(ix,&idx);
902:       ISGetIndices(iy,&idy);
903:       PetscMalloc4(1,VecScatter_Seq_General,&to,nx,PetscInt,&to->vslots,1,VecScatter_Seq_General,&from,nx,PetscInt,&from->vslots);
904:       to->n             = nx;
905: #if defined(PETSC_USE_DEBUG)
906:       VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
907: #endif
908:       PetscMemcpy(to->vslots,idy,nx*sizeof(PetscInt));
909:       from->n           = nx;
910: #if defined(PETSC_USE_DEBUG)
911:       VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
912: #endif
913:        PetscMemcpy(from->vslots,idx,nx*sizeof(PetscInt));
914:       to->type          = VEC_SCATTER_SEQ_GENERAL;
915:       from->type        = VEC_SCATTER_SEQ_GENERAL;
916:       ctx->todata       = (void*)to;
917:       ctx->fromdata     = (void*)from;
918:       ctx->begin        = VecScatterBegin_SGtoSG;
919:       ctx->end          = 0;
920:       ctx->destroy      = VecScatterDestroy_SGtoSG;
921:       ctx->copy         = VecScatterCopy_SGToSG;
922:       PetscInfo(xin,"Special case: sequential vector general scatter\n");
923:       goto functionend;
924:     } else if (ix->type == IS_STRIDE &&  iy->type == IS_STRIDE){
925:       PetscInt               nx,ny,to_first,to_step,from_first,from_step;
926:       VecScatter_Seq_Stride  *from8 = PETSC_NULL,*to8 = PETSC_NULL;

928:       ISGetLocalSize(ix,&nx);
929:       ISGetLocalSize(iy,&ny);
930:       if (nx != ny) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
931:       ISStrideGetInfo(iy,&to_first,&to_step);
932:       ISStrideGetInfo(ix,&from_first,&from_step);
933:       PetscMalloc2(1,VecScatter_Seq_Stride,&to8,1,VecScatter_Seq_Stride,&from8);
934:       to8->n             = nx;
935:       to8->first         = to_first;
936:       to8->step          = to_step;
937:       from8->n           = nx;
938:       from8->first       = from_first;
939:       from8->step        = from_step;
940:       to8->type          = VEC_SCATTER_SEQ_STRIDE;
941:       from8->type        = VEC_SCATTER_SEQ_STRIDE;
942:       ctx->todata       = (void*)to8;
943:       ctx->fromdata     = (void*)from8;
944:       ctx->begin        = VecScatterBegin_SStoSS;
945:       ctx->end          = 0;
946:       ctx->destroy      = VecScatterDestroy_SStoSS;
947:       ctx->copy         = VecScatterCopy_SStoSS;
948:       PetscInfo(xin,"Special case: sequential vector stride to stride\n");
949:       goto functionend;
950:     } else if (ix->type == IS_GENERAL && iy->type == IS_STRIDE){
951:       PetscInt               nx,ny,*idx,first,step;
952:       VecScatter_Seq_General *from9 = PETSC_NULL;
953:       VecScatter_Seq_Stride  *to9 = PETSC_NULL;

955:       ISGetLocalSize(ix,&nx);
956:       ISGetIndices(ix,&idx);
957:       ISGetLocalSize(iy,&ny);
958:       ISStrideGetInfo(iy,&first,&step);
959:       if (nx != ny) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
960:       PetscMalloc3(1,VecScatter_Seq_Stride,&to9,1,VecScatter_Seq_General,&from9,nx,PetscInt,&from9->vslots);
961:       to9->n         = nx;
962:       to9->first     = first;
963:       to9->step      = step;
964:       from9->n       = nx;
965: #if defined(PETSC_USE_DEBUG)
966:       VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
967: #endif
968:       PetscMemcpy(from9->vslots,idx,nx*sizeof(PetscInt));
969:       ctx->todata    = (void*)to9; ctx->fromdata = (void*)from9;
970:       if (step == 1)  ctx->begin = VecScatterBegin_SGtoSS_Stride1;
971:       else            ctx->begin = VecScatterBegin_SGtoSS;
972:       ctx->destroy   = VecScatterDestroy_SGtoSS;
973:       ctx->end       = 0;
974:       ctx->copy      = VecScatterCopy_SGToSS;
975:       to9->type      = VEC_SCATTER_SEQ_STRIDE;
976:       from9->type    = VEC_SCATTER_SEQ_GENERAL;
977:       PetscInfo(xin,"Special case: sequential vector general to stride\n");
978:       goto functionend;
979:     } else if (ix->type == IS_STRIDE && iy->type == IS_GENERAL){
980:       PetscInt               nx,ny,*idy,first,step;
981:       VecScatter_Seq_General *to10 = PETSC_NULL;
982:       VecScatter_Seq_Stride  *from10 = PETSC_NULL;

984:       ISGetLocalSize(ix,&nx);
985:       ISGetIndices(iy,&idy);
986:       ISGetLocalSize(iy,&ny);
987:       ISStrideGetInfo(ix,&first,&step);
988:       if (nx != ny) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
989:       PetscMalloc3(1,VecScatter_Seq_General,&to10,nx,PetscInt,&to10->vslots,1,VecScatter_Seq_Stride,&from10);
990:       from10->n         = nx;
991:       from10->first     = first;
992:       from10->step      = step;
993:       to10->n           = nx;
994: #if defined(PETSC_USE_DEBUG)
995:       VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
996: #endif
997:       PetscMemcpy(to10->vslots,idy,nx*sizeof(PetscInt));
998:       ctx->todata     = (void*)to10;
999:       ctx->fromdata   = (void*)from10;
1000:       if (step == 1) ctx->begin = VecScatterBegin_SStoSG_Stride1;
1001:       else           ctx->begin = VecScatterBegin_SStoSG;
1002:       ctx->destroy    = VecScatterDestroy_SStoSG;
1003:       ctx->end        = 0;
1004:       ctx->copy       = 0;
1005:       to10->type      = VEC_SCATTER_SEQ_GENERAL;
1006:       from10->type    = VEC_SCATTER_SEQ_STRIDE;
1007:       PetscInfo(xin,"Special case: sequential vector stride to general\n");
1008:       goto functionend;
1009:     } else {
1010:       PetscInt               nx,ny,*idx,*idy;
1011:       VecScatter_Seq_General *to11 = PETSC_NULL,*from11 = PETSC_NULL;
1012:       PetscTruth             idnx,idny;

1014:       ISGetLocalSize(ix,&nx);
1015:       ISGetLocalSize(iy,&ny);
1016:       if (nx != ny) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");

1018:       ISIdentity(ix,&idnx);
1019:       ISIdentity(iy,&idny);
1020:       if (idnx && idny) {
1021:         VecScatter_Seq_Stride *to13 = PETSC_NULL,*from13 = PETSC_NULL;
1022:         PetscMalloc2(1,VecScatter_Seq_Stride,&to13,1,VecScatter_Seq_Stride,&from13);
1023:         to13->n           = nx;
1024:         to13->first       = 0;
1025:         to13->step        = 1;
1026:         from13->n         = nx;
1027:         from13->first     = 0;
1028:         from13->step      = 1;
1029:         to13->type        = VEC_SCATTER_SEQ_STRIDE;
1030:         from13->type      = VEC_SCATTER_SEQ_STRIDE;
1031:         ctx->todata       = (void*)to13;
1032:         ctx->fromdata     = (void*)from13;
1033:         ctx->begin        = VecScatterBegin_SStoSS;
1034:         ctx->end          = 0;
1035:         ctx->destroy      = VecScatterDestroy_SStoSS;
1036:         ctx->copy         = VecScatterCopy_SStoSS;
1037:         PetscInfo(xin,"Special case: sequential copy\n");
1038:         goto functionend;
1039:       }

1041:       ISGetIndices(iy,&idy);
1042:       ISGetIndices(ix,&idx);
1043:       PetscMalloc4(1,VecScatter_Seq_General,&to11,nx,PetscInt,&to11->vslots,1,VecScatter_Seq_General,&from11,nx,PetscInt,&from11->vslots);
1044:       to11->n           = nx;
1045: #if defined(PETSC_USE_DEBUG)
1046:       VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1047: #endif
1048:       PetscMemcpy(to11->vslots,idy,nx*sizeof(PetscInt));
1049:       from11->n         = nx;
1050: #if defined(PETSC_USE_DEBUG)
1051:       VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1052: #endif
1053:       PetscMemcpy(from11->vslots,idx,nx*sizeof(PetscInt));
1054:       to11->type        = VEC_SCATTER_SEQ_GENERAL;
1055:       from11->type      = VEC_SCATTER_SEQ_GENERAL;
1056:       ctx->todata       = (void*)to11;
1057:       ctx->fromdata     = (void*)from11;
1058:       ctx->begin        = VecScatterBegin_SGtoSG;
1059:       ctx->end          = 0;
1060:       ctx->destroy      = VecScatterDestroy_SGtoSG;
1061:       ctx->copy         = VecScatterCopy_SGToSG;
1062:       ISRestoreIndices(ix,&idx);
1063:       ISRestoreIndices(iy,&idy);
1064:       PetscInfo(xin,"Sequential vector scatter with block indices\n");
1065:       goto functionend;
1066:     }
1067:   }
1068:   /* ---------------------------------------------------------------------------*/
1069:   if (xin_type == VEC_MPI_ID && yin_type == VEC_SEQ_ID) {

1071:   /* ===========================================================================================================
1072:         Check for special cases
1073:      ==========================================================================================================*/
1074:     islocal = PETSC_FALSE;
1075:     /* special case extracting (subset of) local portion */
1076:     if (ix->type == IS_STRIDE && iy->type == IS_STRIDE){
1077:       PetscInt              nx,ny,to_first,to_step,from_first,from_step;
1078:       PetscInt              start,end;
1079:       VecScatter_Seq_Stride *from12 = PETSC_NULL,*to12 = PETSC_NULL;

1081:       VecGetOwnershipRange(xin,&start,&end);
1082:       ISGetLocalSize(ix,&nx);
1083:       ISStrideGetInfo(ix,&from_first,&from_step);
1084:       ISGetLocalSize(iy,&ny);
1085:       ISStrideGetInfo(iy,&to_first,&to_step);
1086:       if (nx != ny) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1087:       if (ix->min >= start && ix->max < end) islocal = PETSC_TRUE; else islocal = PETSC_FALSE;
1088:       MPI_Allreduce(&islocal,&cando,1,MPI_INT,MPI_LAND,xin->comm);
1089:       if (cando) {
1090:         PetscMalloc2(1,VecScatter_Seq_Stride,&to12,1,VecScatter_Seq_Stride,&from12);
1091:         to12->n             = nx;
1092:         to12->first         = to_first;
1093:         to12->step          = to_step;
1094:         from12->n           = nx;
1095:         from12->first       = from_first-start;
1096:         from12->step        = from_step;
1097:         to12->type          = VEC_SCATTER_SEQ_STRIDE;
1098:         from12->type        = VEC_SCATTER_SEQ_STRIDE;
1099:         ctx->todata         = (void*)to12;
1100:         ctx->fromdata       = (void*)from12;
1101:         ctx->begin          = VecScatterBegin_SStoSS;
1102:         ctx->end            = 0;
1103:         ctx->destroy        = VecScatterDestroy_SStoSS;
1104:         ctx->copy           = VecScatterCopy_SStoSS;
1105:         PetscInfo(xin,"Special case: processors only getting local values\n");
1106:         goto functionend;
1107:       }
1108:     } else {
1109:       MPI_Allreduce(&islocal,&cando,1,MPI_INT,MPI_LAND,xin->comm);
1110:     }

1112:     /* test for special case of all processors getting entire vector */
1113:     totalv = 0;
1114:     if (ix->type == IS_STRIDE && iy->type == IS_STRIDE){
1115:       PetscInt             i,nx,ny,to_first,to_step,from_first,from_step,N;
1116:       PetscMPIInt          *count = PETSC_NULL;
1117:       VecScatter_MPI_ToAll *sto = PETSC_NULL;

1119:       ISGetLocalSize(ix,&nx);
1120:       ISStrideGetInfo(ix,&from_first,&from_step);
1121:       ISGetLocalSize(iy,&ny);
1122:       ISStrideGetInfo(iy,&to_first,&to_step);
1123:       if (nx != ny) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1124:       VecGetSize(xin,&N);
1125:       if (nx != N) {
1126:         totalv = 0;
1127:       } else if (from_first == 0        && from_step == 1 &&
1128:                  from_first == to_first && from_step == to_step){
1129:         totalv = 1;
1130:       } else totalv = 0;
1131:       MPI_Allreduce(&totalv,&cando,1,MPI_INT,MPI_LAND,xin->comm);

1133:       if (cando) {

1135:         MPI_Comm_size(ctx->comm,&size);
1136:         PetscMalloc2(1,VecScatter_MPI_ToAll,&sto,size,PetscMPIInt,&count);
1137:         range = xin->map.range;
1138:         for (i=0; i<size; i++) {
1139:           count[i] = range[i+1] - range[i];
1140:         }
1141:         sto->count        = count;
1142:         sto->work1        = 0;
1143:         sto->work2        = 0;
1144:         sto->type         = VEC_SCATTER_MPI_TOALL;
1145:         ctx->todata       = (void*)sto;
1146:         ctx->fromdata     = 0;
1147:         ctx->begin        = VecScatterBegin_MPI_ToAll;
1148:         ctx->end          = 0;
1149:         ctx->destroy      = VecScatterDestroy_MPI_ToAll;
1150:         ctx->copy         = VecScatterCopy_MPI_ToAll;
1151:         PetscInfo(xin,"Special case: all processors get entire parallel vector\n");
1152:         goto functionend;
1153:       }
1154:     } else {
1155:       MPI_Allreduce(&totalv,&cando,1,MPI_INT,MPI_LAND,xin->comm);
1156:     }

1158:     /* test for special case of processor 0 getting entire vector */
1159:     totalv = 0;
1160:     if (ix->type == IS_STRIDE && iy->type == IS_STRIDE){
1161:       PetscInt             i,nx,ny,to_first,to_step,from_first,from_step,N;
1162:       PetscMPIInt          rank,*count = PETSC_NULL;
1163:       VecScatter_MPI_ToAll *sto = PETSC_NULL;

1165:       PetscObjectGetComm((PetscObject)xin,&comm);
1166:       MPI_Comm_rank(comm,&rank);
1167:       ISGetLocalSize(ix,&nx);
1168:       ISStrideGetInfo(ix,&from_first,&from_step);
1169:       ISGetLocalSize(iy,&ny);
1170:       ISStrideGetInfo(iy,&to_first,&to_step);
1171:       if (nx != ny) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1172:       if (!rank) {
1173:         VecGetSize(xin,&N);
1174:         if (nx != N) {
1175:           totalv = 0;
1176:         } else if (from_first == 0        && from_step == 1 &&
1177:                    from_first == to_first && from_step == to_step){
1178:           totalv = 1;
1179:         } else totalv = 0;
1180:       } else {
1181:         if (!nx) totalv = 1;
1182:         else     totalv = 0;
1183:       }
1184:       MPI_Allreduce(&totalv,&cando,1,MPI_INT,MPI_LAND,xin->comm);

1186:       if (cando) {
1187:         MPI_Comm_size(ctx->comm,&size);
1188:         PetscMalloc2(1,VecScatter_MPI_ToAll,&sto,size,PetscMPIInt,&count);
1189:         range = xin->map.range;
1190:         for (i=0; i<size; i++) {
1191:           count[i] = range[i+1] - range[i];
1192:         }
1193:         sto->count        = count;
1194:         sto->work1        = 0;
1195:         sto->work2        = 0;
1196:         sto->type         = VEC_SCATTER_MPI_TOONE;
1197:         ctx->todata       = (void*)sto;
1198:         ctx->fromdata     = 0;
1199:         ctx->begin        = VecScatterBegin_MPI_ToOne;
1200:         ctx->end          = 0;
1201:         ctx->destroy      = VecScatterDestroy_MPI_ToAll;
1202:         ctx->copy         = VecScatterCopy_MPI_ToAll;
1203:         PetscInfo(xin,"Special case: processor zero gets entire parallel vector, rest get none\n");
1204:         goto functionend;
1205:       }
1206:     } else {
1207:       MPI_Allreduce(&totalv,&cando,1,MPI_INT,MPI_LAND,xin->comm);
1208:     }

1210:     ISBlock(ix,&ixblock);
1211:     ISBlock(iy,&iyblock);
1212:     ISStride(iy,&iystride);
1213:     if (ixblock) {
1214:       /* special case block to block */
1215:       if (iyblock) {
1216:         PetscInt nx,ny,*idx,*idy,bsx,bsy;
1217:         ISBlockGetBlockSize(iy,&bsy);
1218:         ISBlockGetBlockSize(ix,&bsx);
1219:         if (bsx == bsy && VecScatterOptimizedBS(bsx)) {
1220:           ISBlockGetSize(ix,&nx);
1221:           ISBlockGetIndices(ix,&idx);
1222:           ISBlockGetSize(iy,&ny);
1223:           ISBlockGetIndices(iy,&idy);
1224:           if (nx != ny) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1225:           VecScatterCreate_PtoS(nx,idx,ny,idy,xin,yin,bsx,ctx);
1226:           ISBlockRestoreIndices(ix,&idx);
1227:           ISBlockRestoreIndices(iy,&idy);
1228:           PetscInfo(xin,"Special case: blocked indices\n");
1229:           goto functionend;
1230:         }
1231:       /* special case block to stride */
1232:       } else if (iystride) {
1233:         PetscInt ystart,ystride,ysize,bsx;
1234:         ISStrideGetInfo(iy,&ystart,&ystride);
1235:         ISGetLocalSize(iy,&ysize);
1236:         ISBlockGetBlockSize(ix,&bsx);
1237:         /* see if stride index set is equivalent to block index set */
1238:         if (VecScatterOptimizedBS(bsx) && ((ystart % bsx) == 0) && (ystride == 1) && ((ysize % bsx) == 0)) {
1239:           PetscInt nx,*idx,*idy,il;
1240:           ISBlockGetSize(ix,&nx); ISBlockGetIndices(ix,&idx);
1241:           if (ysize != bsx*nx) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1242:           PetscMalloc(nx*sizeof(PetscInt),&idy);
1243:           if (nx) {
1244:             idy[0] = ystart;
1245:             for (il=1; il<nx; il++) idy[il] = idy[il-1] + bsx;
1246:           }
1247:           VecScatterCreate_PtoS(nx,idx,nx,idy,xin,yin,bsx,ctx);
1248:           PetscFree(idy);
1249:           ISBlockRestoreIndices(ix,&idx);
1250:           PetscInfo(xin,"Special case: blocked indices to stride\n");
1251:           goto functionend;
1252:         }
1253:       }
1254:     }
1255:     /* left over general case */
1256:     {
1257:       PetscInt nx,ny,*idx,*idy;
1258:       ISGetLocalSize(ix,&nx);
1259:       ISGetIndices(ix,&idx);
1260:       ISGetLocalSize(iy,&ny);
1261:       ISGetIndices(iy,&idy);
1262:       if (nx != ny) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1263:       VecScatterCreate_PtoS(nx,idx,ny,idy,xin,yin,1,ctx);
1264:       ISRestoreIndices(ix,&idx);
1265:       ISRestoreIndices(iy,&idy);
1266:       PetscInfo(xin,"General case: MPI to Seq\n");
1267:       goto functionend;
1268:     }
1269:   }
1270:   /* ---------------------------------------------------------------------------*/
1271:   if (xin_type == VEC_SEQ_ID && yin_type == VEC_MPI_ID) {
1272:   /* ===========================================================================================================
1273:         Check for special cases
1274:      ==========================================================================================================*/
1275:     /* special case local copy portion */
1276:     islocal = PETSC_FALSE;
1277:     if (ix->type == IS_STRIDE && iy->type == IS_STRIDE){
1278:       PetscInt              nx,ny,to_first,to_step,from_step,start,end,from_first;
1279:       VecScatter_Seq_Stride *from = PETSC_NULL,*to = PETSC_NULL;

1281:       VecGetOwnershipRange(yin,&start,&end);
1282:       ISGetLocalSize(ix,&nx);
1283:       ISStrideGetInfo(ix,&from_first,&from_step);
1284:       ISGetLocalSize(iy,&ny);
1285:       ISStrideGetInfo(iy,&to_first,&to_step);
1286:       if (nx != ny) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1287:       if (iy->min >= start && iy->max < end) islocal = PETSC_TRUE; else islocal = PETSC_FALSE;
1288:       MPI_Allreduce(&islocal,&cando,1,MPI_INT,MPI_LAND,yin->comm);
1289:       if (cando) {
1290:         PetscMalloc2(1,VecScatter_Seq_Stride,&to,1,VecScatter_Seq_Stride,&from);
1291:         to->n             = nx;
1292:         to->first         = to_first-start;
1293:         to->step          = to_step;
1294:         from->n           = nx;
1295:         from->first       = from_first;
1296:         from->step        = from_step;
1297:         to->type          = VEC_SCATTER_SEQ_STRIDE;
1298:         from->type        = VEC_SCATTER_SEQ_STRIDE;
1299:         ctx->todata       = (void*)to;
1300:         ctx->fromdata     = (void*)from;
1301:         ctx->begin        = VecScatterBegin_SStoSS;
1302:         ctx->end          = 0;
1303:         ctx->destroy      = VecScatterDestroy_SStoSS;
1304:         ctx->copy         = VecScatterCopy_SStoSS;
1305:         PetscInfo(xin,"Special case: sequential stride to MPI stride\n");
1306:         goto functionend;
1307:       }
1308:     } else {
1309:       MPI_Allreduce(&islocal,&cando,1,MPI_INT,MPI_LAND,yin->comm);
1310:     }
1311:     if (ix->type == IS_BLOCK && iy->type == IS_STRIDE){
1312:       PetscInt ystart,ystride,ysize,bsx;
1313:       ISStrideGetInfo(iy,&ystart,&ystride);
1314:       ISGetLocalSize(iy,&ysize);
1315:       ISBlockGetBlockSize(ix,&bsx);
1316:       /* see if stride index set is equivalent to block index set */
1317:       if (VecScatterOptimizedBS(bsx) && ((ystart % bsx) == 0) && (ystride == 1) && ((ysize % bsx) == 0)) {
1318:         PetscInt nx,*idx,*idy,il;
1319:         ISBlockGetSize(ix,&nx); ISBlockGetIndices(ix,&idx);
1320:         if (ysize != bsx*nx) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1321:         PetscMalloc(nx*sizeof(PetscInt),&idy);
1322:         if (nx) {
1323:           idy[0] = ystart;
1324:           for (il=1; il<nx; il++) idy[il] = idy[il-1] + bsx;
1325:         }
1326:         VecScatterCreate_StoP(nx,idx,nx,idy,xin,yin,bsx,ctx);
1327:         PetscFree(idy);
1328:         ISBlockRestoreIndices(ix,&idx);
1329:         PetscInfo(xin,"Special case: Blocked indices to stride\n");
1330:         goto functionend;
1331:       }
1332:     }

1334:     /* general case */
1335:     {
1336:       PetscInt nx,ny,*idx,*idy;
1337:       ISGetLocalSize(ix,&nx);
1338:       ISGetIndices(ix,&idx);
1339:       ISGetLocalSize(iy,&ny);
1340:       ISGetIndices(iy,&idy);
1341:       if (nx != ny) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1342:       VecScatterCreate_StoP(nx,idx,ny,idy,xin,yin,1,ctx);
1343:       ISRestoreIndices(ix,&idx);
1344:       ISRestoreIndices(iy,&idy);
1345:       PetscInfo(xin,"General case: Seq to MPI\n");
1346:       goto functionend;
1347:     }
1348:   }
1349:   /* ---------------------------------------------------------------------------*/
1350:   if (xin_type == VEC_MPI_ID && yin_type == VEC_MPI_ID) {
1351:     /* no special cases for now */
1352:     PetscInt nx,ny,*idx,*idy;
1353:     ISGetLocalSize(ix,&nx);
1354:     ISGetIndices(ix,&idx);
1355:     ISGetLocalSize(iy,&ny);
1356:     ISGetIndices(iy,&idy);
1357:     if (nx != ny) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1358:     VecScatterCreate_PtoP(nx,idx,ny,idy,xin,yin,ctx);
1359:     ISRestoreIndices(ix,&idx);
1360:     ISRestoreIndices(iy,&idy);
1361:     PetscInfo(xin,"General case: MPI to MPI\n");
1362:     goto functionend;
1363:   }

1365:   functionend:
1366:   *newctx = ctx;
1367:   if (tix) {ISDestroy(tix);}
1368:   if (tiy) {ISDestroy(tiy);}
1369:   PetscOptionsHasName(PETSC_NULL,"-vecscatter_view_info",&flag);
1370:   if (flag) {
1371:     PetscViewer viewer;
1372:     PetscViewerASCIIGetStdout(comm,&viewer);
1373:     PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1374:     VecScatterView(ctx,viewer);
1375:     PetscViewerPopFormat(viewer);
1376:   }
1377:   PetscOptionsHasName(PETSC_NULL,"-vecscatter_view",&flag);
1378:   if (flag) {
1379:     PetscViewer viewer;
1380:     PetscViewerASCIIGetStdout(comm,&viewer);
1381:     VecScatterView(ctx,viewer);
1382:   }
1383:   return(0);
1384: }

1386: /* ------------------------------------------------------------------*/
1389: /*@
1390:    VecScatterGetMerged - Returns true if the scatter is completed in the VecScatterBegin()
1391:       and the VecScatterEnd() does nothing

1393:    Not Collective

1395:    Input Parameter:
1396: .   ctx - scatter context created with VecScatterCreate()

1398:    Output Parameter:
1399: .   flg - PETSC_TRUE if the VecScatterBegin/End() are all done during the VecScatterBegin()

1401:    Level: developer

1403: .seealso: VecScatterCreate(), VecScatterEnd(), VecScatterBegin()
1404: @*/
1405: PetscErrorCode  VecScatterGetMerged(VecScatter ctx,PetscTruth *flg)
1406: {
1409:   *flg = ctx->beginandendtogether;
1410:   return(0);
1411: }

1415: /*@
1416:    VecScatterBegin - Begins a generalized scatter from one vector to
1417:    another. Complete the scattering phase with VecScatterEnd().

1419:    Collective on VecScatter and Vec

1421:    Input Parameters:
1422: +  x - the vector from which we scatter
1423: .  y - the vector to which we scatter
1424: .  addv - either ADD_VALUES or INSERT_VALUES, with INSERT_VALUES mode any location 
1425:           not scattered to retains its old value; i.e. the vector is NOT first zeroed.
1426: .  mode - the scattering mode, usually SCATTER_FORWARD.  The available modes are:
1427:     SCATTER_FORWARD or SCATTER_REVERSE
1428: -  inctx - scatter context generated by VecScatterCreate()

1430:    Level: intermediate

1432:    Options Database
1433: +   -vecscatter_rr   - use ready receiver mode (i.e. receives are post BEFORE sends)
1434: .   -vecscatter_ssend  - use MPI_Ssend() instead of MPI_Send()
1435: .   -vecscatter_packtogether - packs all the message before sending any and receivers all
1436:                                before sending. Default for the alltoall versions.
1437: .   -vecscatter_sendfirst - post ALL sends before posting receives (cannot be used with -vecscatter_rr)
1438: .   -vecscatter_alltoallv - use MPI_Alltoallv() instead of sends and receives
1439: -   -vecscatter_alltoallw  - use MPI_Alltoallw() instead of MPI_Alltoallv() for INSERT_VALUES

1441:    Notes:
1442:    The vectors x and y need not be the same vectors used in the call 
1443:    to VecScatterCreate(), but x must have the same parallel data layout
1444:    as that passed in as the x to VecScatterCreate(), similarly for the y.
1445:    Most likely they have been obtained from VecDuplicate().

1447:    You cannot change the values in the input vector between the calls to VecScatterBegin()
1448:    and VecScatterEnd().

1450:    If you use SCATTER_REVERSE the first two arguments should be reversed, from 
1451:    the SCATTER_FORWARD.
1452:    
1453:    y[iy[i]] = x[ix[i]], for i=0,...,ni-1

1455:    This scatter is far more general than the conventional
1456:    scatter, since it can be a gather or a scatter or a combination,
1457:    depending on the indices ix and iy.  If x is a parallel vector and y
1458:    is sequential, VecScatterBegin() can serve to gather values to a
1459:    single processor.  Similarly, if y is parallel and x sequential, the
1460:    routine can scatter from one processor to many processors.

1462:    Concepts: scatter^between vectors
1463:    Concepts: gather^between vectors

1465: .seealso: VecScatterCreate(), VecScatterEnd()
1466: @*/
1467: PetscErrorCode  VecScatterBegin(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter inctx)
1468: {
1470: #if defined(PETSC_USE_DEBUG)
1471:   PetscInt      to_n,from_n;
1472: #endif

1478:   if (inctx->inuse) SETERRQ(PETSC_ERR_ARG_WRONGSTATE," Scatter ctx already in use");
1479: #if defined(PETSC_USE_DEBUG)
1480:   /*
1481:      Error checking to make sure these vectors match the vectors used
1482:    to create the vector scatter context. -1 in the from_n and to_n indicate the
1483:    vector lengths are unknown (for example with mapped scatters) and thus 
1484:    no error checking is performed.
1485:   */
1486:   if (inctx->from_n >= 0 && inctx->to_n >= 0) {
1487:     VecGetLocalSize(x,&from_n);
1488:     VecGetLocalSize(y,&to_n);
1489:     if (mode & SCATTER_REVERSE) {
1490:       if (to_n != inctx->from_n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter reverse and vector to != ctx from size)",to_n,inctx->from_n);
1491:       if (from_n != inctx->to_n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter reverse and vector from != ctx to size)",from_n,inctx->to_n);
1492:     } else {
1493:       if (to_n != inctx->to_n)     SETERRQ2(PETSC_ERR_ARG_SIZ,"Vector wrong size % for scatter %D (scatter forward and vector to != ctx to size)",to_n,inctx->to_n);
1494:       if (from_n != inctx->from_n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vector wrong size % for scatter %D (scatter forward and vector from != ctx from size)",from_n,inctx->from_n);
1495:     }
1496:   }
1497: #endif

1499:   inctx->inuse = PETSC_TRUE;
1501:   (*inctx->begin)(x,y,addv,mode,inctx);
1502:   if (inctx->beginandendtogether && inctx->end) {
1503:     inctx->inuse = PETSC_FALSE;
1504:     (*inctx->end)(x,y,addv,mode,inctx);
1505:   }
1507:   return(0);
1508: }

1510: /* --------------------------------------------------------------------*/
1513: /*@
1514:    VecScatterEnd - Ends a generalized scatter from one vector to another.  Call
1515:    after first calling VecScatterBegin().

1517:    Collective on VecScatter and Vec

1519:    Input Parameters:
1520: +  x - the vector from which we scatter
1521: .  y - the vector to which we scatter
1522: .  addv - either ADD_VALUES or INSERT_VALUES.
1523: .  mode - the scattering mode, usually SCATTER_FORWARD.  The available modes are:
1524:      SCATTER_FORWARD, SCATTER_REVERSE
1525: -  ctx - scatter context generated by VecScatterCreate()

1527:    Level: intermediate

1529:    Notes:
1530:    If you use SCATTER_REVERSE the first two arguments should be reversed, from 
1531:    the SCATTER_FORWARD.
1532:    y[iy[i]] = x[ix[i]], for i=0,...,ni-1

1534: .seealso: VecScatterBegin(), VecScatterCreate()
1535: @*/
1536: PetscErrorCode  VecScatterEnd(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
1537: {

1544:   ctx->inuse = PETSC_FALSE;
1545:   if (!ctx->end) return(0);
1546:   if (!ctx->beginandendtogether) {
1548:     (*(ctx)->end)(x,y,addv,mode,ctx);
1550:   }
1551:   return(0);
1552: }

1556: /*@
1557:    VecScatterDestroy - Destroys a scatter context created by 
1558:    VecScatterCreate().

1560:    Collective on VecScatter

1562:    Input Parameter:
1563: .  ctx - the scatter context

1565:    Level: intermediate

1567: .seealso: VecScatterCreate(), VecScatterCopy()
1568: @*/
1569: PetscErrorCode  VecScatterDestroy(VecScatter ctx)
1570: {

1575:   if (--ctx->refct > 0) return(0);

1577:   /* if memory was published with AMS then destroy it */
1578:   PetscObjectDepublish(ctx);

1580:   (*ctx->destroy)(ctx);
1581:   return(0);
1582: }

1586: /*@
1587:    VecScatterCopy - Makes a copy of a scatter context.

1589:    Collective on VecScatter

1591:    Input Parameter:
1592: .  sctx - the scatter context

1594:    Output Parameter:
1595: .  ctx - the context copy

1597:    Level: advanced

1599: .seealso: VecScatterCreate(), VecScatterDestroy()
1600: @*/
1601: PetscErrorCode  VecScatterCopy(VecScatter sctx,VecScatter *ctx)
1602: {

1608:   if (!sctx->copy) SETERRQ(PETSC_ERR_SUP,"Cannot copy this type");
1609:   PetscHeaderCreate(*ctx,_p_VecScatter,int,VEC_SCATTER_COOKIE,0,"VecScatter",sctx->comm,VecScatterDestroy,VecScatterView);
1610:   (*ctx)->to_n   = sctx->to_n;
1611:   (*ctx)->from_n = sctx->from_n;
1612:   (*sctx->copy)(sctx,*ctx);
1613:   return(0);
1614: }


1617: /* ------------------------------------------------------------------*/
1620: /*@
1621:    VecScatterView - Views a vector scatter context.

1623:    Collective on VecScatter

1625:    Input Parameters:
1626: +  ctx - the scatter context
1627: -  viewer - the viewer for displaying the context

1629:    Level: intermediate

1631: @*/
1632: PetscErrorCode  VecScatterView(VecScatter ctx,PetscViewer viewer)
1633: {

1638:   if (!viewer) {
1639:     PetscViewerASCIIGetStdout(ctx->comm,&viewer);
1640:   }
1642:   if (!ctx->view) SETERRQ(PETSC_ERR_SUP,"Cannot view this type of scatter context yet");

1644:   (*ctx->view)(ctx,viewer);
1645:   return(0);
1646: }

1650: /*@
1651:    VecScatterRemap - Remaps the "from" and "to" indices in a 
1652:    vector scatter context. FOR EXPERTS ONLY!

1654:    Collective on VecScatter

1656:    Input Parameters:
1657: +  scat - vector scatter context
1658: .  from - remapping for "from" indices (may be PETSC_NULL)
1659: -  to   - remapping for "to" indices (may be PETSC_NULL)

1661:    Level: developer

1663:    Notes: In the parallel case the todata is actually the indices
1664:           from which the data is TAKEN! The from stuff is where the 
1665:           data is finally put. This is VERY VERY confusing!

1667:           In the sequential case the todata is the indices where the 
1668:           data is put and the fromdata is where it is taken from.
1669:           This is backwards from the paralllel case! CRY! CRY! CRY!

1671: @*/
1672: PetscErrorCode  VecScatterRemap(VecScatter scat,PetscInt *rto,PetscInt *rfrom)
1673: {
1674:   VecScatter_Seq_General *to,*from;
1675:   VecScatter_MPI_General *mto;
1676:   PetscInt               i;


1683:   from = (VecScatter_Seq_General *)scat->fromdata;
1684:   mto  = (VecScatter_MPI_General *)scat->todata;

1686:   if (mto->type == VEC_SCATTER_MPI_TOALL) SETERRQ(PETSC_ERR_ARG_SIZ,"Not for to all scatter");

1688:   if (rto) {
1689:     if (mto->type == VEC_SCATTER_MPI_GENERAL) {
1690:       /* handle off processor parts */
1691:       for (i=0; i<mto->starts[mto->n]; i++) {
1692:         mto->indices[i] = rto[mto->indices[i]];
1693:       }
1694:       /* handle local part */
1695:       to = &mto->local;
1696:       for (i=0; i<to->n; i++) {
1697:         to->vslots[i] = rto[to->vslots[i]];
1698:       }
1699:     } else if (from->type == VEC_SCATTER_SEQ_GENERAL) {
1700:       for (i=0; i<from->n; i++) {
1701:         from->vslots[i] = rto[from->vslots[i]];
1702:       }
1703:     } else if (from->type == VEC_SCATTER_SEQ_STRIDE) {
1704:       VecScatter_Seq_Stride *sto = (VecScatter_Seq_Stride*)from;
1705: 
1706:       /* if the remapping is the identity and stride is identity then skip remap */
1707:       if (sto->step == 1 && sto->first == 0) {
1708:         for (i=0; i<sto->n; i++) {
1709:           if (rto[i] != i) {
1710:             SETERRQ(PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1711:           }
1712:         }
1713:       } else SETERRQ(PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1714:     } else SETERRQ(PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1715:   }

1717:   if (rfrom) {
1718:     SETERRQ(PETSC_ERR_SUP,"Unable to remap the FROM in scatters yet");
1719:   }

1721:   /*
1722:      Mark then vector lengths as unknown because we do not know the 
1723:    lengths of the remapped vectors
1724:   */
1725:   scat->from_n = -1;
1726:   scat->to_n   = -1;

1728:   return(0);
1729: }