Actual source code: vscat.c

petsc-3.6.1 2015-07-22
Report Typos and Errors
  2: /*
  3:      Code for creating scatters between vectors. This file
  4:   includes the code for scattering between sequential vectors and
  5:   some special cases for parallel scatters.
  6: */

  8: #include <petsc/private/vecimpl.h>    /*I   "petscvec.h"    I*/

 10: #if defined(PETSC_HAVE_CUSP)
 11: PETSC_INTERN PetscErrorCode VecScatterCUSPIndicesCreate_StoS(PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt*,PetscInt*,PetscCUSPIndices*);
 12: PETSC_INTERN PetscErrorCode VecScatterCUSPIndicesDestroy(PetscCUSPIndices*);
 13: PETSC_INTERN PetscErrorCode VecScatterCUSP_StoS(Vec,Vec,PetscCUSPIndices,InsertMode,ScatterMode);
 14: #endif

 16: /* Logging support */
 17: PetscClassId VEC_SCATTER_CLASSID;

 19: #if defined(PETSC_USE_DEBUG)
 20: /*
 21:      Checks if any indices are less than zero and generates an error
 22: */
 25: static PetscErrorCode VecScatterCheckIndices_Private(PetscInt nmax,PetscInt n,const PetscInt *idx)
 26: {
 27:   PetscInt i;

 30:   for (i=0; i<n; i++) {
 31:     if (idx[i] < 0)     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative index %D at %D location",idx[i],i);
 32:     if (idx[i] >= nmax) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D at %D location greater than max %D",idx[i],i,nmax);
 33:   }
 34:   return(0);
 35: }
 36: #endif

 38: /*
 39:       This is special scatter code for when the entire parallel vector is copied to each processor.

 41:    This code was written by Cameron Cooper, Occidental College, Fall 1995,
 42:    will working at ANL as a SERS student.
 43: */
 46: PetscErrorCode VecScatterBegin_MPI_ToAll(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
 47: {
 49:   PetscInt       yy_n,xx_n;
 50:   PetscScalar    *xv,*yv;

 53:   VecGetLocalSize(y,&yy_n);
 54:   VecGetLocalSize(x,&xx_n);
 55:   VecGetArrayPair(x,y,&xv,&yv);

 57:   if (mode & SCATTER_REVERSE) {
 58:     PetscScalar          *xvt,*xvt2;
 59:     VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
 60:     PetscInt             i;
 61:     PetscMPIInt          *disply = scat->displx;

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

110:     if (addv == INSERT_VALUES) {
111:       MPI_Allgatherv(xv,xx_n,MPIU_SCALAR,yv,scat->count,displx,MPIU_SCALAR,PetscObjectComm((PetscObject)ctx));
112:     } else {
113:       if (scat->work1) yvt = scat->work1;
114:       else {
115:         PetscMalloc1(yy_n,&yvt);
116:         scat->work1 = yvt;
117:       }
118:       MPI_Allgatherv(xv,xx_n,MPIU_SCALAR,yvt,scat->count,displx,MPIU_SCALAR,PetscObjectComm((PetscObject)ctx));
119:       if (addv == ADD_VALUES) {
120:         for (i=0; i<yy_n; i++) yv[i] += yvt[i];
121: #if !defined(PETSC_USE_COMPLEX)
122:       } else if (addv == MAX_VALUES) {
123:         for (i=0; i<yy_n; i++) yv[i] = PetscMax(yv[i],yvt[i]);
124: #endif
125:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
126:     }
127:   }
128:   VecRestoreArrayPair(x,y,&xv,&yv);
129:   return(0);
130: }

134: PetscErrorCode VecScatterView_MPI_ToAll(VecScatter in,PetscViewer viewer)
135: {
137:   PetscBool      isascii;

140:   PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
141:   if (isascii) {
142:     PetscViewerASCIIPrintf(viewer,"Entire parallel vector is copied to each process\n");
143:   }
144:   return(0);
145: }

147: /*
148:       This is special scatter code for when the entire parallel vector is  copied to processor 0.

150: */
153: PetscErrorCode VecScatterBegin_MPI_ToOne(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
154: {
155:   PetscErrorCode    ierr;
156:   PetscMPIInt       rank;
157:   PetscInt          yy_n,xx_n;
158:   PetscScalar       *yv;
159:   const PetscScalar *xv;
160:   MPI_Comm          comm;

163:   VecGetLocalSize(y,&yy_n);
164:   VecGetLocalSize(x,&xx_n);
165:   VecGetArrayRead(x,&xv);
166:   VecGetArray(y,&yv);

168:   PetscObjectGetComm((PetscObject)x,&comm);
169:   MPI_Comm_rank(comm,&rank);

171:   /* --------  Reverse scatter; spread from processor 0 to other processors */
172:   if (mode & SCATTER_REVERSE) {
173:     PetscScalar          *yvt;
174:     VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
175:     PetscInt             i;
176:     PetscMPIInt          *disply = scat->displx;

178:     if (addv == INSERT_VALUES) {
179:       MPI_Scatterv((PetscScalar*)xv,scat->count,disply,MPIU_SCALAR,yv,yy_n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
180:     } else {
181:       if (scat->work2) yvt = scat->work2;
182:       else {
183:         PetscMalloc1(xx_n,&yvt);
184:         scat->work2 = yvt;
185:       }
186:       MPI_Scatterv((PetscScalar*)xv,scat->count,disply,MPIU_SCALAR,yvt,yy_n,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
187:       if (addv == ADD_VALUES) {
188:         for (i=0; i<yy_n; i++) yv[i] += yvt[i];
189: #if !defined(PETSC_USE_COMPLEX)
190:       } else if (addv == MAX_VALUES) {
191:         for (i=0; i<yy_n; i++) yv[i] = PetscMax(yv[i],yvt[i]);
192: #endif
193:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
194:     }
195:   /* ---------  Forward scatter; gather all values onto processor 0 */
196:   } else {
197:     PetscScalar          *yvt  = 0;
198:     VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
199:     PetscInt             i;
200:     PetscMPIInt          *displx = scat->displx;

202:     if (addv == INSERT_VALUES) {
203:       MPI_Gatherv((PetscScalar*)xv,xx_n,MPIU_SCALAR,yv,scat->count,displx,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
204:     } else {
205:       if (!rank) {
206:         if (scat->work1) yvt = scat->work1;
207:         else {
208:           PetscMalloc1(yy_n,&yvt);
209:           scat->work1 = yvt;
210:         }
211:       }
212:       MPI_Gatherv((PetscScalar*)xv,xx_n,MPIU_SCALAR,yvt,scat->count,displx,MPIU_SCALAR,0,PetscObjectComm((PetscObject)ctx));
213:       if (!rank) {
214:         if (addv == ADD_VALUES) {
215:           for (i=0; i<yy_n; i++) yv[i] += yvt[i];
216: #if !defined(PETSC_USE_COMPLEX)
217:         } else if (addv == MAX_VALUES) {
218:           for (i=0; i<yy_n; i++) yv[i] = PetscMax(yv[i],yvt[i]);
219: #endif
220:         }  else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
221:       }
222:     }
223:   }
224:   VecRestoreArrayRead(x,&xv);
225:   VecRestoreArray(y,&yv);
226:   return(0);
227: }

229: /*
230:        The follow to are used for both VecScatterBegin_MPI_ToAll() and VecScatterBegin_MPI_ToOne()
231: */
234: PetscErrorCode VecScatterDestroy_MPI_ToAll(VecScatter ctx)
235: {
236:   VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
237:   PetscErrorCode       ierr;

240:   PetscFree(scat->work1);
241:   PetscFree(scat->work2);
242:   PetscFree3(ctx->todata,scat->count,scat->displx);
243:   return(0);
244: }

248: PetscErrorCode VecScatterDestroy_SGToSG(VecScatter ctx)
249: {

253:   PetscFree2(((VecScatter_Seq_General*)ctx->todata)->vslots,((VecScatter_Seq_General*)ctx->fromdata)->vslots);
254:   PetscFree2(ctx->todata,ctx->fromdata);
255:   return(0);
256: }

260: PetscErrorCode VecScatterDestroy_SGToSS(VecScatter ctx)
261: {

265:   PetscFree(((VecScatter_Seq_General*)ctx->fromdata)->vslots);
266:   PetscFree2(ctx->todata,ctx->fromdata);
267:   return(0);
268: }

272: PetscErrorCode VecScatterDestroy_SSToSG(VecScatter ctx)
273: {

277:   PetscFree(((VecScatter_Seq_General*)ctx->todata)->vslots);
278:   PetscFree2(ctx->todata,ctx->fromdata);
279:   return(0);
280: }

284: PetscErrorCode VecScatterDestroy_SSToSS(VecScatter ctx)
285: {

289:   PetscFree2(ctx->todata,ctx->fromdata);
290:   return(0);
291: }

293: /* -------------------------------------------------------------------------*/
296: PetscErrorCode VecScatterCopy_MPI_ToAll(VecScatter in,VecScatter out)
297: {
298:   VecScatter_MPI_ToAll *in_to = (VecScatter_MPI_ToAll*)in->todata,*sto;
299:   PetscErrorCode       ierr;
300:   PetscMPIInt          size,*count,*displx;
301:   PetscInt             i;

304:   out->ops->begin   = in->ops->begin;
305:   out->ops->end     = in->ops->end;
306:   out->ops->copy    = in->ops->copy;
307:   out->ops->destroy = in->ops->destroy;
308:   out->ops->view    = in->ops->view;

310:   MPI_Comm_size(PetscObjectComm((PetscObject)out),&size);
311:   PetscMalloc3(1,&sto,size,&count,size,&displx);
312:   sto->type   = in_to->type;
313:   sto->count  = count;
314:   sto->displx = displx;
315:   for (i=0; i<size; i++) {
316:     sto->count[i]  = in_to->count[i];
317:     sto->displx[i] = in_to->displx[i];
318:   }
319:   sto->work1    = 0;
320:   sto->work2    = 0;
321:   out->todata   = (void*)sto;
322:   out->fromdata = (void*)0;
323:   return(0);
324: }

326: /* --------------------------------------------------------------------------------------*/
327: /*
328:         Scatter: sequential general to sequential general
329: */
332: PetscErrorCode VecScatterBegin_SGToSG(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
333: {
334:   VecScatter_Seq_General *gen_to   = (VecScatter_Seq_General*)ctx->todata;
335:   VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
336:   PetscErrorCode         ierr;
337:   PetscInt               i,n = gen_from->n,*fslots,*tslots;
338:   PetscScalar            *xv,*yv;

341: #if defined(PETSC_HAVE_CUSP)
342:   if (x->valid_GPU_array == PETSC_CUSP_GPU) {
343:     /* create the scatter indices if not done already */
344:     if (!ctx->spptr) {
345:       PetscInt tofirst = 0,tostep = 0,fromfirst = 0,fromstep = 0;
346:       fslots = gen_from->vslots;
347:       tslots = gen_to->vslots;
348:       VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
349:     }
350:     /* next do the scatter */
351:     VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
352:     return(0);
353:   }
354: #endif

356:   VecGetArrayPair(x,y,&xv,&yv);
357:   if (mode & SCATTER_REVERSE) {
358:     gen_to   = (VecScatter_Seq_General*)ctx->fromdata;
359:     gen_from = (VecScatter_Seq_General*)ctx->todata;
360:   }
361:   fslots = gen_from->vslots;
362:   tslots = gen_to->vslots;

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

377: /*
378:     Scatter: sequential general to sequential stride 1
379: */
382: PetscErrorCode VecScatterBegin_SGToSS_Stride1(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
383: {
384:   VecScatter_Seq_Stride  *gen_to   = (VecScatter_Seq_Stride*)ctx->todata;
385:   VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
386:   PetscInt               i,n = gen_from->n,*fslots = gen_from->vslots;
387:   PetscErrorCode         ierr;
388:   PetscInt               first = gen_to->first;
389:   PetscScalar            *xv,*yv;

392: #if defined(PETSC_HAVE_CUSP)
393:   if (x->valid_GPU_array == PETSC_CUSP_GPU) {
394:     /* create the scatter indices if not done already */
395:     if (!ctx->spptr) {
396:       PetscInt tofirst = first,tostep = 1,fromfirst = 0,fromstep = 0;
397:       PetscInt *tslots = 0;
398:       VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
399:     }
400:     /* next do the scatter */
401:     VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
402:     return(0);
403:   }
404: #endif

406:   VecGetArrayPair(x,y,&xv,&yv);
407:   if (mode & SCATTER_REVERSE) {
408:     xv += first;
409:     if (addv == INSERT_VALUES) {
410:       for (i=0; i<n; i++) yv[fslots[i]] = xv[i];
411:     } else if (addv == ADD_VALUES) {
412:       for (i=0; i<n; i++) yv[fslots[i]] += xv[i];
413: #if !defined(PETSC_USE_COMPLEX)
414:     } else if (addv == MAX_VALUES) {
415:       for (i=0; i<n; i++) yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[i]);
416: #endif
417:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
418:   } else {
419:     yv += first;
420:     if (addv == INSERT_VALUES) {
421:       for (i=0; i<n; i++) yv[i] = xv[fslots[i]];
422:     } else if (addv == ADD_VALUES) {
423:       for (i=0; i<n; i++) yv[i] += xv[fslots[i]];
424: #if !defined(PETSC_USE_COMPLEX)
425:     } else if (addv == MAX_VALUES) {
426:       for (i=0; i<n; i++) yv[i] = PetscMax(yv[i],xv[fslots[i]]);
427: #endif
428:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
429:   }
430:   VecRestoreArrayPair(x,y,&xv,&yv);
431:   return(0);
432: }

434: /*
435:    Scatter: sequential general to sequential stride
436: */
439: PetscErrorCode VecScatterBegin_SGToSS(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
440: {
441:   VecScatter_Seq_Stride  *gen_to   = (VecScatter_Seq_Stride*)ctx->todata;
442:   VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
443:   PetscInt               i,n = gen_from->n,*fslots = gen_from->vslots;
444:   PetscErrorCode         ierr;
445:   PetscInt               first = gen_to->first,step = gen_to->step;
446:   PetscScalar            *xv,*yv;

449: #if defined(PETSC_HAVE_CUSP)
450:   if (x->valid_GPU_array == PETSC_CUSP_GPU) {
451:     /* create the scatter indices if not done already */
452:     if (!ctx->spptr) {
453:       PetscInt tofirst = first,tostep = step,fromfirst = 0,fromstep = 0;
454:       PetscInt * tslots = 0;
455:       VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
456:     }
457:     /* next do the scatter */
458:     VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
459:     return(0);
460:   }
461: #endif

463:   VecGetArrayPair(x,y,&xv,&yv);
464:   if (mode & SCATTER_REVERSE) {
465:     if (addv == INSERT_VALUES) {
466:       for (i=0; i<n; i++) yv[fslots[i]] = xv[first + i*step];
467:     } else if (addv == ADD_VALUES) {
468:       for (i=0; i<n; i++) yv[fslots[i]] += xv[first + i*step];
469: #if !defined(PETSC_USE_COMPLEX)
470:     } else if (addv == MAX_VALUES) {
471:       for (i=0; i<n; i++) yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[first + i*step]);
472: #endif
473:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
474:   } else {
475:     if (addv == INSERT_VALUES) {
476:       for (i=0; i<n; i++) yv[first + i*step] = xv[fslots[i]];
477:     } else if (addv == ADD_VALUES) {
478:       for (i=0; i<n; i++) yv[first + i*step] += xv[fslots[i]];
479: #if !defined(PETSC_USE_COMPLEX)
480:     } else if (addv == MAX_VALUES) {
481:       for (i=0; i<n; i++) yv[first + i*step] = PetscMax(yv[first + i*step],xv[fslots[i]]);
482: #endif
483:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
484:   }
485:   VecRestoreArrayPair(x,y,&xv,&yv);
486:   return(0);
487: }

489: /*
490:     Scatter: sequential stride 1 to sequential general
491: */
494: PetscErrorCode VecScatterBegin_SSToSG_Stride1(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
495: {
496:   VecScatter_Seq_Stride  *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
497:   VecScatter_Seq_General *gen_to   = (VecScatter_Seq_General*)ctx->todata;
498:   PetscInt               i,n = gen_from->n,*tslots = gen_to->vslots;
499:   PetscErrorCode         ierr;
500:   PetscInt               first = gen_from->first;
501:   PetscScalar            *xv,*yv;

504: #if defined(PETSC_HAVE_CUSP)
505:   if (x->valid_GPU_array == PETSC_CUSP_GPU) {
506:     /* create the scatter indices if not done already */
507:     if (!ctx->spptr) {
508:       PetscInt tofirst = 0,tostep = 0,fromfirst = first,fromstep = 1;
509:       PetscInt *fslots = 0;
510:       VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
511:     }
512:     /* next do the scatter */
513:     VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
514:     return(0);
515:   }
516: #endif

518:   VecGetArrayPair(x,y,&xv,&yv);
519:   if (mode & SCATTER_REVERSE) {
520:     yv += first;
521:     if (addv == INSERT_VALUES) {
522:       for (i=0; i<n; i++) yv[i] = xv[tslots[i]];
523:     } else  if (addv == ADD_VALUES) {
524:       for (i=0; i<n; i++) yv[i] += xv[tslots[i]];
525: #if !defined(PETSC_USE_COMPLEX)
526:     } else  if (addv == MAX_VALUES) {
527:       for (i=0; i<n; i++) yv[i] = PetscMax(yv[i],xv[tslots[i]]);
528: #endif
529:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
530:   } else {
531:     xv += first;
532:     if (addv == INSERT_VALUES) {
533:       for (i=0; i<n; i++) yv[tslots[i]] = xv[i];
534:     } else  if (addv == ADD_VALUES) {
535:       for (i=0; i<n; i++) yv[tslots[i]] += xv[i];
536: #if !defined(PETSC_USE_COMPLEX)
537:     } else  if (addv == MAX_VALUES) {
538:       for (i=0; i<n; i++) yv[tslots[i]] = PetscMax(yv[tslots[i]],xv[i]);
539: #endif
540:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
541:   }
542:   VecRestoreArrayPair(x,y,&xv,&yv);
543:   return(0);
544: }

548: /*
549:    Scatter: sequential stride to sequential general
550: */
551: PetscErrorCode VecScatterBegin_SSToSG(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
552: {
553:   VecScatter_Seq_Stride  *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
554:   VecScatter_Seq_General *gen_to   = (VecScatter_Seq_General*)ctx->todata;
555:   PetscInt               i,n = gen_from->n,*tslots = gen_to->vslots;
556:   PetscErrorCode         ierr;
557:   PetscInt               first = gen_from->first,step = gen_from->step;
558:   PetscScalar            *xv,*yv;

561: #if defined(PETSC_HAVE_CUSP)
562:   if (x->valid_GPU_array == PETSC_CUSP_GPU) {
563:     /* create the scatter indices if not done already */
564:     if (!ctx->spptr) {
565:       PetscInt tofirst = 0,tostep = 0,fromfirst = first,fromstep = step;
566:       PetscInt *fslots = 0;
567:       VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
568:     }
569:     /* next do the scatter */
570:     VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
571:     return(0);
572:   }
573: #endif

575:   VecGetArrayPair(x,y,&xv,&yv);
576:   if (mode & SCATTER_REVERSE) {
577:     if (addv == INSERT_VALUES) {
578:       for (i=0; i<n; i++) yv[first + i*step] = xv[tslots[i]];
579:     } else if (addv == ADD_VALUES) {
580:       for (i=0; i<n; i++) yv[first + i*step] += xv[tslots[i]];
581: #if !defined(PETSC_USE_COMPLEX)
582:     } else if (addv == MAX_VALUES) {
583:       for (i=0; i<n; i++) yv[first + i*step] = PetscMax(yv[first + i*step],xv[tslots[i]]);
584: #endif
585:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
586:   } else {
587:     if (addv == INSERT_VALUES) {
588:       for (i=0; i<n; i++) yv[tslots[i]] = xv[first + i*step];
589:     } else if (addv == ADD_VALUES) {
590:       for (i=0; i<n; i++) yv[tslots[i]] += xv[first + i*step];
591: #if !defined(PETSC_USE_COMPLEX)
592:     } else if (addv == MAX_VALUES) {
593:       for (i=0; i<n; i++) yv[tslots[i]] = PetscMax(yv[tslots[i]],xv[first + i*step]);
594: #endif
595:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
596:   }
597:   VecRestoreArrayPair(x,y,&xv,&yv);
598:   return(0);
599: }

603: PetscErrorCode VecScatterView_SSToSG(VecScatter in,PetscViewer viewer)
604: {
605:   PetscErrorCode         ierr;
606:   VecScatter_Seq_Stride  *in_from = (VecScatter_Seq_Stride*)in->todata;
607:   VecScatter_Seq_General *in_to   = (VecScatter_Seq_General*)in->fromdata;
608:   PetscInt               i;
609:   PetscBool              isascii;

612:   PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
613:   if (isascii) {
614:     PetscViewerASCIIPrintf(viewer,"Sequential stride to general scatter\n");
615:     for (i=0; i<in_to->n; i++) {
616:       PetscViewerASCIIPrintf(viewer,"%D to %D\n",in_from->first + in_from->step*i,in_to->vslots[i]);
617:     }
618:   }
619:   return(0);
620: }
621: /*
622:      Scatter: sequential stride to sequential stride
623: */
626: PetscErrorCode VecScatterBegin_SSToSS(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
627: {
628:   VecScatter_Seq_Stride *gen_to   = (VecScatter_Seq_Stride*)ctx->todata;
629:   VecScatter_Seq_Stride *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
630:   PetscInt              i,n = gen_from->n,to_first = gen_to->first,to_step = gen_to->step;
631:   PetscErrorCode        ierr;
632:   PetscInt              from_first = gen_from->first,from_step = gen_from->step;
633:   PetscScalar           *xv,*yv;

636: #if defined(PETSC_HAVE_CUSP)
637:   if (x->valid_GPU_array == PETSC_CUSP_GPU) {
638:     /* create the scatter indices if not done already */
639:     if (!ctx->spptr) {
640:       PetscInt *tslots = 0,*fslots = 0;
641:       VecScatterCUSPIndicesCreate_StoS(n,to_first,from_first,to_step,from_step,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
642:     }
643:     /* next do the scatter */
644:     VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
645:     return(0);
646:   }
647: #endif

649:   VecGetArrayPair(x,y,&xv,&yv);
650:   if (mode & SCATTER_REVERSE) {
651:     from_first = gen_to->first;
652:     to_first   = gen_from->first;
653:     from_step  = gen_to->step;
654:     to_step    = gen_from->step;
655:   }

657:   if (addv == INSERT_VALUES) {
658:     if (to_step == 1 && from_step == 1) {
659:       PetscMemcpy(yv+to_first,xv+from_first,n*sizeof(PetscScalar));
660:     } else  {
661:       for (i=0; i<n; i++) yv[to_first + i*to_step] = xv[from_first+i*from_step];
662:     }
663:   } else if (addv == ADD_VALUES) {
664:     if (to_step == 1 && from_step == 1) {
665:       yv += to_first; xv += from_first;
666:       for (i=0; i<n; i++) yv[i] += xv[i];
667:     } else {
668:       for (i=0; i<n; i++) yv[to_first + i*to_step] += xv[from_first+i*from_step];
669:     }
670: #if !defined(PETSC_USE_COMPLEX)
671:   } else if (addv == MAX_VALUES) {
672:     if (to_step == 1 && from_step == 1) {
673:       yv += to_first; xv += from_first;
674:       for (i=0; i<n; i++) yv[i] = PetscMax(yv[i],xv[i]);
675:     } else {
676:       for (i=0; i<n; i++) yv[to_first + i*to_step] = PetscMax(yv[to_first + i*to_step],xv[from_first+i*from_step]);
677:     }
678: #endif
679:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
680:   VecRestoreArrayPair(x,y,&xv,&yv);
681:   return(0);
682: }

684: /* --------------------------------------------------------------------------*/


689: PetscErrorCode VecScatterCopy_SGToSG(VecScatter in,VecScatter out)
690: {
691:   PetscErrorCode         ierr;
692:   VecScatter_Seq_General *in_to   = (VecScatter_Seq_General*)in->todata,*out_to = NULL;
693:   VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from = NULL;

696:   out->ops->begin   = in->ops->begin;
697:   out->ops->end     = in->ops->end;
698:   out->ops->copy    = in->ops->copy;
699:   out->ops->destroy = in->ops->destroy;
700:   out->ops->view    = in->ops->view;

702:   PetscMalloc2(1,&out_to,1,&out_from);
703:   PetscMalloc2(in_to->n,&out_to->vslots,in_from->n,&out_from->vslots);
704:   out_to->n                    = in_to->n;
705:   out_to->type                 = in_to->type;
706:   out_to->nonmatching_computed = PETSC_FALSE;
707:   out_to->slots_nonmatching    = 0;
708:   out_to->is_copy              = PETSC_FALSE;
709:   PetscMemcpy(out_to->vslots,in_to->vslots,(out_to->n)*sizeof(PetscInt));

711:   out_from->n                    = in_from->n;
712:   out_from->type                 = in_from->type;
713:   out_from->nonmatching_computed = PETSC_FALSE;
714:   out_from->slots_nonmatching    = 0;
715:   out_from->is_copy              = PETSC_FALSE;
716:   PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));

718:   out->todata   = (void*)out_to;
719:   out->fromdata = (void*)out_from;
720:   return(0);
721: }

725: PetscErrorCode VecScatterView_SGToSG(VecScatter in,PetscViewer viewer)
726: {
727:   PetscErrorCode         ierr;
728:   VecScatter_Seq_General *in_to   = (VecScatter_Seq_General*)in->todata;
729:   VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata;
730:   PetscInt               i;
731:   PetscBool              isascii;

734:   PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
735:   if (isascii) {
736:     PetscViewerASCIIPrintf(viewer,"Sequential general scatter\n");
737:     for (i=0; i<in_to->n; i++) {
738:       PetscViewerASCIIPrintf(viewer,"%D to %D\n",in_from->vslots[i],in_to->vslots[i]);
739:     }
740:   }
741:   return(0);
742: }


747: PetscErrorCode VecScatterCopy_SGToSS(VecScatter in,VecScatter out)
748: {
749:   PetscErrorCode         ierr;
750:   VecScatter_Seq_Stride  *in_to   = (VecScatter_Seq_Stride*)in->todata,*out_to = NULL;
751:   VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from = NULL;

754:   out->ops->begin   = in->ops->begin;
755:   out->ops->end     = in->ops->end;
756:   out->ops->copy    = in->ops->copy;
757:   out->ops->destroy = in->ops->destroy;
758:   out->ops->view    = in->ops->view;

760:   PetscMalloc2(1,&out_to,1,&out_from);
761:   PetscMalloc1(in_from->n,&out_from->vslots);
762:   out_to->n     = in_to->n;
763:   out_to->type  = in_to->type;
764:   out_to->first = in_to->first;
765:   out_to->step  = in_to->step;
766:   out_to->type  = in_to->type;

768:   out_from->n                    = in_from->n;
769:   out_from->type                 = in_from->type;
770:   out_from->nonmatching_computed = PETSC_FALSE;
771:   out_from->slots_nonmatching    = 0;
772:   out_from->is_copy              = PETSC_FALSE;
773:   PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));

775:   out->todata   = (void*)out_to;
776:   out->fromdata = (void*)out_from;
777:   return(0);
778: }

782: PetscErrorCode VecScatterView_SGToSS(VecScatter in,PetscViewer viewer)
783: {
784:   PetscErrorCode         ierr;
785:   VecScatter_Seq_Stride  *in_to   = (VecScatter_Seq_Stride*)in->todata;
786:   VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata;
787:   PetscInt               i;
788:   PetscBool              isascii;

791:   PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
792:   if (isascii) {
793:     PetscViewerASCIIPrintf(viewer,"Sequential general scatter to stride\n");
794:     for (i=0; i<in_to->n; i++) {
795:       PetscViewerASCIIPrintf(viewer,"%D to %D\n",in_from->vslots[i],in_to->first + in_to->step*i);
796:     }
797:   }
798:   return(0);
799: }

801: /* --------------------------------------------------------------------------*/
802: /*
803:     Scatter: parallel to sequential vector, sequential strides for both.
804: */
807: PetscErrorCode VecScatterCopy_SSToSS(VecScatter in,VecScatter out)
808: {
809:   VecScatter_Seq_Stride *in_to   = (VecScatter_Seq_Stride*)in->todata,*out_to = NULL;
810:   VecScatter_Seq_Stride *in_from = (VecScatter_Seq_Stride*)in->fromdata,*out_from = NULL;
811:   PetscErrorCode        ierr;

814:   out->ops->begin   = in->ops->begin;
815:   out->ops->end     = in->ops->end;
816:   out->ops->copy    = in->ops->copy;
817:   out->ops->destroy = in->ops->destroy;
818:   out->ops->view    = in->ops->view;

820:   PetscMalloc2(1,&out_to,1,&out_from);
821:   out_to->n       = in_to->n;
822:   out_to->type    = in_to->type;
823:   out_to->first   = in_to->first;
824:   out_to->step    = in_to->step;
825:   out_to->type    = in_to->type;
826:   out_from->n     = in_from->n;
827:   out_from->type  = in_from->type;
828:   out_from->first = in_from->first;
829:   out_from->step  = in_from->step;
830:   out_from->type  = in_from->type;
831:   out->todata     = (void*)out_to;
832:   out->fromdata   = (void*)out_from;
833:   return(0);
834: }

838: PetscErrorCode VecScatterView_SSToSS(VecScatter in,PetscViewer viewer)
839: {
840:   VecScatter_Seq_Stride *in_to   = (VecScatter_Seq_Stride*)in->todata;
841:   VecScatter_Seq_Stride *in_from = (VecScatter_Seq_Stride*)in->fromdata;
842:   PetscErrorCode        ierr;
843:   PetscBool             isascii;

846:   PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
847:   if (isascii) {
848:     PetscViewerASCIIPrintf(viewer,"Sequential stride count %D start %D step to start %D stride %D\n",in_to->n,in_to->first,in_to->step,in_from->first,in_from->step);
849:   }
850:   return(0);
851: }


854: extern PetscErrorCode VecScatterCreate_PtoS(PetscInt,const PetscInt*,PetscInt,const PetscInt*,Vec,Vec,PetscInt,VecScatter);
855: extern PetscErrorCode VecScatterCreate_PtoP(PetscInt,const PetscInt*,PetscInt,const PetscInt*,Vec,Vec,PetscInt,VecScatter);
856: extern PetscErrorCode VecScatterCreate_StoP(PetscInt,const PetscInt*,PetscInt,const PetscInt*,Vec,Vec,PetscInt,VecScatter);

858: /* =======================================================================*/
859: #define VEC_SEQ_ID 0
860: #define VEC_MPI_ID 1
861: #define IS_GENERAL_ID 0
862: #define IS_STRIDE_ID  1
863: #define IS_BLOCK_ID   2

865: /*
866:    Blocksizes we have optimized scatters for
867: */

869: #define VecScatterOptimizedBS(mbs) (2 <= mbs)

871: PetscErrorCode  VecScatterCreateEmpty(MPI_Comm comm,VecScatter *newctx)
872: {
873:   VecScatter     ctx;

877:   PetscHeaderCreate(ctx,VEC_SCATTER_CLASSID,"VecScatter","VecScatter","Vec",comm,VecScatterDestroy,VecScatterView);
878:   ctx->inuse               = PETSC_FALSE;
879:   ctx->beginandendtogether = PETSC_FALSE;
880:   PetscOptionsGetBool(NULL,"-vecscatter_merge",&ctx->beginandendtogether,NULL);
881:   if (ctx->beginandendtogether) {
882:     PetscInfo(ctx,"Using combined (merged) vector scatter begin and end\n");
883:   }
884:   ctx->packtogether = PETSC_FALSE;
885:   PetscOptionsGetBool(NULL,"-vecscatter_packtogether",&ctx->packtogether,NULL);
886:   if (ctx->packtogether) {
887:     PetscInfo(ctx,"Pack all messages before sending\n");
888:   }
889:   *newctx = ctx;
890:   return(0);
891: }

895: /*@C
896:    VecScatterCreate - Creates a vector scatter context.

898:    Collective on Vec

900:    Input Parameters:
901: +  xin - a vector that defines the shape (parallel data layout of the vector)
902:          of vectors from which we scatter
903: .  yin - a vector that defines the shape (parallel data layout of the vector)
904:          of vectors to which we scatter
905: .  ix - the indices of xin to scatter (if NULL scatters all values)
906: -  iy - the indices of yin to hold results (if NULL fills entire vector yin)

908:    Output Parameter:
909: .  newctx - location to store the new scatter context

911:    Options Database Keys: (uses regular MPI_Sends by default)
912: +  -vecscatter_view         - Prints detail of communications
913: .  -vecscatter_view ::ascii_info    - Print less details about communication
914: .  -vecscatter_ssend        - Uses MPI_Ssend_init() instead of MPI_Send_init()
915: .  -vecscatter_rsend           - use ready receiver mode for MPI sends
916: .  -vecscatter_merge        - VecScatterBegin() handles all of the communication, VecScatterEnd() is a nop
917:                               eliminates the chance for overlap of computation and communication
918: .  -vecscatter_sendfirst    - Posts sends before receives
919: .  -vecscatter_packtogether - Pack all messages before sending, receive all messages before unpacking
920: .  -vecscatter_alltoall     - Uses MPI all to all communication for scatter
921: .  -vecscatter_window       - Use MPI 2 window operations to move data
922: .  -vecscatter_nopack       - Avoid packing to work vector when possible (if used with -vecscatter_alltoall then will use MPI_Alltoallw()
923: -  -vecscatter_reproduce    - insure that the order of the communications are done the same for each scatter, this under certain circumstances
924:                               will make the results of scatters deterministic when otherwise they are not (it may be slower also).

926: $
927: $                                                                                    --When packing is used--
928: $                               MPI Datatypes (no packing)  sendfirst   merge        packtogether  persistent*
929: $                                _nopack                   _sendfirst    _merge      _packtogether                -vecscatter_
930: $ ----------------------------------------------------------------------------------------------------------------------------
931: $    Message passing    Send       p                           X            X           X         always
932: $                      Ssend       p                           X            X           X         always          _ssend
933: $                      Rsend       p                        nonsense        X           X         always          _rsend
934: $    AlltoAll  v or w              X                        nonsense     always         X         nonsense        _alltoall
935: $    MPI_Win                       p                        nonsense        p           p         nonsense        _window
936: $
937: $   Since persistent sends and receives require a constant memory address they can only be used when data is packed into the work vector
938: $   because the in and out array may be different for each call to VecScatterBegin/End().
939: $
940: $    p indicates possible, but not implemented. X indicates implemented
941: $

943:     Level: intermediate

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

954:    Currently the MPI_Send(), MPI_Ssend() and MPI_Rsend() all use PERSISTENT versions.
955:    (this unfortunately requires that the same in and out arrays be used for each use, this
956:     is why when not using MPI_alltoallw() we always need to pack the input into the work array before sending
957:     and unpack upon receeving instead of using MPI datatypes to avoid the packing/unpacking).

959:    Both ix and iy cannot be NULL at the same time.

961:    Concepts: scatter^between vectors
962:    Concepts: gather^between vectors

964: .seealso: VecScatterDestroy(), VecScatterCreateToAll(), VecScatterCreateToZero()
965: @*/
966: PetscErrorCode  VecScatterCreate(Vec xin,IS ix,Vec yin,IS iy,VecScatter *newctx)
967: {
968:   VecScatter        ctx;
969:   PetscErrorCode    ierr;
970:   PetscMPIInt       size;
971:   PetscInt          xin_type = VEC_SEQ_ID,yin_type = VEC_SEQ_ID,*range;
972:   PetscInt          ix_type  = IS_GENERAL_ID,iy_type = IS_GENERAL_ID;
973:   MPI_Comm          comm,ycomm;
974:   PetscBool         totalv,ixblock,iyblock,iystride,islocal,cando,flag;
975:   IS                tix = 0,tiy = 0;

978:   if (!ix && !iy) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_SUP,"Cannot pass default in for both input and output indices");

980:   /*
981:       Determine if the vectors are "parallel", ie. it shares a comm with other processors, or
982:       sequential (it does not share a comm). The difference is that parallel vectors treat the
983:       index set as providing indices in the global parallel numbering of the vector, with
984:       sequential vectors treat the index set as providing indices in the local sequential
985:       numbering
986:   */
987:   PetscObjectGetComm((PetscObject)xin,&comm);
988:   MPI_Comm_size(comm,&size);
989:   if (size > 1) xin_type = VEC_MPI_ID;

991:   PetscObjectGetComm((PetscObject)yin,&ycomm);
992:   MPI_Comm_size(ycomm,&size);
993:   if (size > 1) {comm = ycomm; yin_type = VEC_MPI_ID;}

995:   /* generate the Scatter context */
996:   PetscHeaderCreate(ctx,VEC_SCATTER_CLASSID,"VecScatter","VecScatter","Vec",comm,VecScatterDestroy,VecScatterView);
997:   ctx->inuse               = PETSC_FALSE;

999:   ctx->beginandendtogether = PETSC_FALSE;
1000:   PetscOptionsGetBool(NULL,"-vecscatter_merge",&ctx->beginandendtogether,NULL);
1001:   if (ctx->beginandendtogether) {
1002:     PetscInfo(ctx,"Using combined (merged) vector scatter begin and end\n");
1003:   }
1004:   ctx->packtogether = PETSC_FALSE;
1005:   PetscOptionsGetBool(NULL,"-vecscatter_packtogether",&ctx->packtogether,NULL);
1006:   if (ctx->packtogether) {
1007:     PetscInfo(ctx,"Pack all messages before sending\n");
1008:   }

1010:   VecGetLocalSize(xin,&ctx->from_n);
1011:   VecGetLocalSize(yin,&ctx->to_n);

1013:   /*
1014:       if ix or iy is not included; assume just grabbing entire vector
1015:   */
1016:   if (!ix && xin_type == VEC_SEQ_ID) {
1017:     ISCreateStride(comm,ctx->from_n,0,1,&ix);
1018:     tix  = ix;
1019:   } else if (!ix && xin_type == VEC_MPI_ID) {
1020:     if (yin_type == VEC_MPI_ID) {
1021:       PetscInt ntmp, low;
1022:       VecGetLocalSize(xin,&ntmp);
1023:       VecGetOwnershipRange(xin,&low,NULL);
1024:       ISCreateStride(comm,ntmp,low,1,&ix);
1025:     } else {
1026:       PetscInt Ntmp;
1027:       VecGetSize(xin,&Ntmp);
1028:       ISCreateStride(comm,Ntmp,0,1,&ix);
1029:     }
1030:     tix = ix;
1031:   } else if (!ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ix not given, but not Seq or MPI vector");

1033:   if (!iy && yin_type == VEC_SEQ_ID) {
1034:     ISCreateStride(comm,ctx->to_n,0,1,&iy);
1035:     tiy  = iy;
1036:   } else if (!iy && yin_type == VEC_MPI_ID) {
1037:     if (xin_type == VEC_MPI_ID) {
1038:       PetscInt ntmp, low;
1039:       VecGetLocalSize(yin,&ntmp);
1040:       VecGetOwnershipRange(yin,&low,NULL);
1041:       ISCreateStride(comm,ntmp,low,1,&iy);
1042:     } else {
1043:       PetscInt Ntmp;
1044:       VecGetSize(yin,&Ntmp);
1045:       ISCreateStride(comm,Ntmp,0,1,&iy);
1046:     }
1047:     tiy = iy;
1048:   } else if (!iy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"iy not given, but not Seq or MPI vector");

1050:   /*
1051:      Determine types of index sets
1052:   */
1053:   PetscObjectTypeCompare((PetscObject)ix,ISBLOCK,&flag);
1054:   if (flag) ix_type = IS_BLOCK_ID;
1055:   PetscObjectTypeCompare((PetscObject)iy,ISBLOCK,&flag);
1056:   if (flag) iy_type = IS_BLOCK_ID;
1057:   PetscObjectTypeCompare((PetscObject)ix,ISSTRIDE,&flag);
1058:   if (flag) ix_type = IS_STRIDE_ID;
1059:   PetscObjectTypeCompare((PetscObject)iy,ISSTRIDE,&flag);
1060:   if (flag) iy_type = IS_STRIDE_ID;

1062:   /* ===========================================================================================================
1063:         Check for special cases
1064:      ==========================================================================================================*/
1065:   /* ---------------------------------------------------------------------------*/
1066:   if (xin_type == VEC_SEQ_ID && yin_type == VEC_SEQ_ID) {
1067:     if (ix_type == IS_GENERAL_ID && iy_type == IS_GENERAL_ID) {
1068:       PetscInt               nx,ny;
1069:       const PetscInt         *idx,*idy;
1070:       VecScatter_Seq_General *to = NULL,*from = NULL;

1072:       ISGetLocalSize(ix,&nx);
1073:       ISGetLocalSize(iy,&ny);
1074:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1075:       ISGetIndices(ix,&idx);
1076:       ISGetIndices(iy,&idy);
1077:       PetscMalloc2(1,&to,1,&from);
1078:       PetscMalloc2(nx,&to->vslots,nx,&from->vslots);
1079:       to->n = nx;
1080: #if defined(PETSC_USE_DEBUG)
1081:       VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1082: #endif
1083:       PetscMemcpy(to->vslots,idy,nx*sizeof(PetscInt));
1084:       from->n = nx;
1085: #if defined(PETSC_USE_DEBUG)
1086:       VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1087: #endif
1088:        PetscMemcpy(from->vslots,idx,nx*sizeof(PetscInt));
1089:       to->type          = VEC_SCATTER_SEQ_GENERAL;
1090:       from->type        = VEC_SCATTER_SEQ_GENERAL;
1091:       ctx->todata       = (void*)to;
1092:       ctx->fromdata     = (void*)from;
1093:       ctx->ops->begin   = VecScatterBegin_SGToSG;
1094:       ctx->ops->end     = 0;
1095:       ctx->ops->destroy = VecScatterDestroy_SGToSG;
1096:       ctx->ops->copy    = VecScatterCopy_SGToSG;
1097:       ctx->ops->view    = VecScatterView_SGToSG;
1098:       PetscInfo(xin,"Special case: sequential vector general scatter\n");
1099:       goto functionend;
1100:     } else if (ix_type == IS_STRIDE_ID &&  iy_type == IS_STRIDE_ID) {
1101:       PetscInt              nx,ny,to_first,to_step,from_first,from_step;
1102:       VecScatter_Seq_Stride *from8 = NULL,*to8 = NULL;

1104:       ISGetLocalSize(ix,&nx);
1105:       ISGetLocalSize(iy,&ny);
1106:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1107:       ISStrideGetInfo(iy,&to_first,&to_step);
1108:       ISStrideGetInfo(ix,&from_first,&from_step);
1109:       PetscMalloc2(1,&to8,1,&from8);
1110:       to8->n            = nx;
1111:       to8->first        = to_first;
1112:       to8->step         = to_step;
1113:       from8->n          = nx;
1114:       from8->first      = from_first;
1115:       from8->step       = from_step;
1116:       to8->type         = VEC_SCATTER_SEQ_STRIDE;
1117:       from8->type       = VEC_SCATTER_SEQ_STRIDE;
1118:       ctx->todata       = (void*)to8;
1119:       ctx->fromdata     = (void*)from8;
1120:       ctx->ops->begin   = VecScatterBegin_SSToSS;
1121:       ctx->ops->end     = 0;
1122:       ctx->ops->destroy = VecScatterDestroy_SSToSS;
1123:       ctx->ops->copy    = VecScatterCopy_SSToSS;
1124:       ctx->ops->view    = VecScatterView_SSToSS;
1125:       PetscInfo(xin,"Special case: sequential vector stride to stride\n");
1126:       goto functionend;
1127:     } else if (ix_type == IS_GENERAL_ID && iy_type == IS_STRIDE_ID) {
1128:       PetscInt               nx,ny,first,step;
1129:       const PetscInt         *idx;
1130:       VecScatter_Seq_General *from9 = NULL;
1131:       VecScatter_Seq_Stride  *to9   = NULL;

1133:       ISGetLocalSize(ix,&nx);
1134:       ISGetIndices(ix,&idx);
1135:       ISGetLocalSize(iy,&ny);
1136:       ISStrideGetInfo(iy,&first,&step);
1137:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1138:       PetscMalloc2(1,&to9,1,&from9);
1139:       PetscMalloc1(nx,&from9->vslots);
1140:       to9->n     = nx;
1141:       to9->first = first;
1142:       to9->step  = step;
1143:       from9->n   = nx;
1144: #if defined(PETSC_USE_DEBUG)
1145:       VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1146: #endif
1147:       PetscMemcpy(from9->vslots,idx,nx*sizeof(PetscInt));
1148:       ctx->todata    = (void*)to9; ctx->fromdata = (void*)from9;
1149:       if (step == 1)  ctx->ops->begin = VecScatterBegin_SGToSS_Stride1;
1150:       else            ctx->ops->begin = VecScatterBegin_SGToSS;
1151:       ctx->ops->destroy   = VecScatterDestroy_SGToSS;
1152:       ctx->ops->end       = 0;
1153:       ctx->ops->copy      = VecScatterCopy_SGToSS;
1154:       ctx->ops->view      = VecScatterView_SGToSS;
1155:       to9->type      = VEC_SCATTER_SEQ_STRIDE;
1156:       from9->type    = VEC_SCATTER_SEQ_GENERAL;
1157:       PetscInfo(xin,"Special case: sequential vector general to stride\n");
1158:       goto functionend;
1159:     } else if (ix_type == IS_STRIDE_ID && iy_type == IS_GENERAL_ID) {
1160:       PetscInt               nx,ny,first,step;
1161:       const PetscInt         *idy;
1162:       VecScatter_Seq_General *to10 = NULL;
1163:       VecScatter_Seq_Stride  *from10 = NULL;

1165:       ISGetLocalSize(ix,&nx);
1166:       ISGetIndices(iy,&idy);
1167:       ISGetLocalSize(iy,&ny);
1168:       ISStrideGetInfo(ix,&first,&step);
1169:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1170:       PetscMalloc2(1,&to10,1,&from10);
1171:       PetscMalloc1(nx,&to10->vslots);
1172:       from10->n     = nx;
1173:       from10->first = first;
1174:       from10->step  = step;
1175:       to10->n       = nx;
1176: #if defined(PETSC_USE_DEBUG)
1177:       VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1178: #endif
1179:       PetscMemcpy(to10->vslots,idy,nx*sizeof(PetscInt));
1180:       ctx->todata   = (void*)to10;
1181:       ctx->fromdata = (void*)from10;
1182:       if (step == 1) ctx->ops->begin = VecScatterBegin_SSToSG_Stride1;
1183:       else           ctx->ops->begin = VecScatterBegin_SSToSG;
1184:       ctx->ops->destroy = VecScatterDestroy_SSToSG;
1185:       ctx->ops->end     = 0;
1186:       ctx->ops->copy    = 0;
1187:       ctx->ops->view    = VecScatterView_SSToSG;
1188:       to10->type   = VEC_SCATTER_SEQ_GENERAL;
1189:       from10->type = VEC_SCATTER_SEQ_STRIDE;
1190:       PetscInfo(xin,"Special case: sequential vector stride to general\n");
1191:       goto functionend;
1192:     } else {
1193:       PetscInt               nx,ny;
1194:       const PetscInt         *idx,*idy;
1195:       VecScatter_Seq_General *to11 = NULL,*from11 = NULL;
1196:       PetscBool              idnx,idny;

1198:       ISGetLocalSize(ix,&nx);
1199:       ISGetLocalSize(iy,&ny);
1200:       if (nx != ny) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match, in %D out %D",nx,ny);

1202:       ISIdentity(ix,&idnx);
1203:       ISIdentity(iy,&idny);
1204:       if (idnx && idny) {
1205:         VecScatter_Seq_Stride *to13 = NULL,*from13 = NULL;
1206:         PetscMalloc2(1,&to13,1,&from13);
1207:         to13->n       = nx;
1208:         to13->first   = 0;
1209:         to13->step    = 1;
1210:         from13->n     = nx;
1211:         from13->first = 0;
1212:         from13->step  = 1;
1213:         to13->type    = VEC_SCATTER_SEQ_STRIDE;
1214:         from13->type  = VEC_SCATTER_SEQ_STRIDE;
1215:         ctx->todata   = (void*)to13;
1216:         ctx->fromdata = (void*)from13;
1217:         ctx->ops->begin    = VecScatterBegin_SSToSS;
1218:         ctx->ops->end      = 0;
1219:         ctx->ops->destroy  = VecScatterDestroy_SSToSS;
1220:         ctx->ops->copy     = VecScatterCopy_SSToSS;
1221:         ctx->ops->view     = VecScatterView_SSToSS;
1222:         PetscInfo(xin,"Special case: sequential copy\n");
1223:         goto functionend;
1224:       }

1226:       ISGetIndices(iy,&idy);
1227:       ISGetIndices(ix,&idx);
1228:       PetscMalloc2(1,&to11,1,&from11);
1229:       PetscMalloc2(nx,&to11->vslots,nx,&from11->vslots);
1230:       to11->n = nx;
1231: #if defined(PETSC_USE_DEBUG)
1232:       VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1233: #endif
1234:       PetscMemcpy(to11->vslots,idy,nx*sizeof(PetscInt));
1235:       from11->n = nx;
1236: #if defined(PETSC_USE_DEBUG)
1237:       VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1238: #endif
1239:       PetscMemcpy(from11->vslots,idx,nx*sizeof(PetscInt));
1240:       to11->type        = VEC_SCATTER_SEQ_GENERAL;
1241:       from11->type      = VEC_SCATTER_SEQ_GENERAL;
1242:       ctx->todata       = (void*)to11;
1243:       ctx->fromdata     = (void*)from11;
1244:       ctx->ops->begin   = VecScatterBegin_SGToSG;
1245:       ctx->ops->end     = 0;
1246:       ctx->ops->destroy = VecScatterDestroy_SGToSG;
1247:       ctx->ops->copy    = VecScatterCopy_SGToSG;
1248:       ctx->ops->view    = VecScatterView_SGToSG;
1249:       ISRestoreIndices(ix,&idx);
1250:       ISRestoreIndices(iy,&idy);
1251:       PetscInfo(xin,"Sequential vector scatter with block indices\n");
1252:       goto functionend;
1253:     }
1254:   }
1255:   /* ---------------------------------------------------------------------------*/
1256:   if (xin_type == VEC_MPI_ID && yin_type == VEC_SEQ_ID) {

1258:     /* ===========================================================================================================
1259:           Check for special cases
1260:        ==========================================================================================================*/
1261:     islocal = PETSC_FALSE;
1262:     /* special case extracting (subset of) local portion */
1263:     if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1264:       PetscInt              nx,ny,to_first,to_step,from_first,from_step;
1265:       PetscInt              start,end,min,max;
1266:       VecScatter_Seq_Stride *from12 = NULL,*to12 = NULL;

1268:       VecGetOwnershipRange(xin,&start,&end);
1269:       ISGetLocalSize(ix,&nx);
1270:       ISStrideGetInfo(ix,&from_first,&from_step);
1271:       ISGetLocalSize(iy,&ny);
1272:       ISStrideGetInfo(iy,&to_first,&to_step);
1273:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1274:       ISGetMinMax(ix,&min,&max);
1275:       if (min >= start && max < end) islocal = PETSC_TRUE;
1276:       else islocal = PETSC_FALSE;
1277:       MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1278:       if (cando) {
1279:         PetscMalloc2(1,&to12,1,&from12);
1280:         to12->n            = nx;
1281:         to12->first        = to_first;
1282:         to12->step         = to_step;
1283:         from12->n          = nx;
1284:         from12->first      = from_first-start;
1285:         from12->step       = from_step;
1286:         to12->type         = VEC_SCATTER_SEQ_STRIDE;
1287:         from12->type       = VEC_SCATTER_SEQ_STRIDE;
1288:         ctx->todata        = (void*)to12;
1289:         ctx->fromdata      = (void*)from12;
1290:         ctx->ops->begin    = VecScatterBegin_SSToSS;
1291:         ctx->ops->end      = 0;
1292:         ctx->ops->destroy  = VecScatterDestroy_SSToSS;
1293:         ctx->ops->copy     = VecScatterCopy_SSToSS;
1294:         ctx->ops->view     = VecScatterView_SSToSS;
1295:         PetscInfo(xin,"Special case: processors only getting local values\n");
1296:         goto functionend;
1297:       }
1298:     } else {
1299:       MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1300:     }

1302:     /* test for special case of all processors getting entire vector */
1303:     /* contains check that PetscMPIInt can handle the sizes needed */
1304:     totalv = PETSC_FALSE;
1305:     if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1306:       PetscInt             i,nx,ny,to_first,to_step,from_first,from_step,N;
1307:       PetscMPIInt          *count = NULL,*displx;
1308:       VecScatter_MPI_ToAll *sto   = NULL;

1310:       ISGetLocalSize(ix,&nx);
1311:       ISStrideGetInfo(ix,&from_first,&from_step);
1312:       ISGetLocalSize(iy,&ny);
1313:       ISStrideGetInfo(iy,&to_first,&to_step);
1314:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1315:       VecGetSize(xin,&N);
1316:       if (nx != N) totalv = PETSC_FALSE;
1317:       else if (from_first == 0 && from_step == 1 && from_first == to_first && from_step == to_step) totalv = PETSC_TRUE;
1318:       else totalv = PETSC_FALSE;
1319:       MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));

1321: #if defined(PETSC_USE_64BIT_INDICES)
1322:       if (cando && (yin->map->N < PETSC_MPI_INT_MAX)) {
1323: #else
1324:       if (cando) {
1325: #endif
1326:         MPI_Comm_size(PetscObjectComm((PetscObject)ctx),&size);
1327:         PetscMalloc3(1,&sto,size,&count,size,&displx);
1328:         range = xin->map->range;
1329:         for (i=0; i<size; i++) {
1330:           PetscMPIIntCast(range[i+1] - range[i],count+i);
1331:           PetscMPIIntCast(range[i],displx+i);
1332:         }
1333:         sto->count        = count;
1334:         sto->displx       = displx;
1335:         sto->work1        = 0;
1336:         sto->work2        = 0;
1337:         sto->type         = VEC_SCATTER_MPI_TOALL;
1338:         ctx->todata       = (void*)sto;
1339:         ctx->fromdata     = 0;
1340:         ctx->ops->begin   = VecScatterBegin_MPI_ToAll;
1341:         ctx->ops->end     = 0;
1342:         ctx->ops->destroy = VecScatterDestroy_MPI_ToAll;
1343:         ctx->ops->copy    = VecScatterCopy_MPI_ToAll;
1344:         ctx->ops->view    = VecScatterView_MPI_ToAll;
1345:         PetscInfo(xin,"Special case: all processors get entire parallel vector\n");
1346:         goto functionend;
1347:       }
1348:     } else {
1349:       MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1350:     }

1352:     /* test for special case of processor 0 getting entire vector */
1353:     /* contains check that PetscMPIInt can handle the sizes needed */
1354:     totalv = PETSC_FALSE;
1355:     if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1356:       PetscInt             i,nx,ny,to_first,to_step,from_first,from_step,N;
1357:       PetscMPIInt          rank,*count = NULL,*displx;
1358:       VecScatter_MPI_ToAll *sto = NULL;

1360:       PetscObjectGetComm((PetscObject)xin,&comm);
1361:       MPI_Comm_rank(comm,&rank);
1362:       ISGetLocalSize(ix,&nx);
1363:       ISStrideGetInfo(ix,&from_first,&from_step);
1364:       ISGetLocalSize(iy,&ny);
1365:       ISStrideGetInfo(iy,&to_first,&to_step);
1366:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1367:       if (!rank) {
1368:         VecGetSize(xin,&N);
1369:         if (nx != N) totalv = PETSC_FALSE;
1370:         else if (from_first == 0        && from_step == 1 &&
1371:                  from_first == to_first && from_step == to_step) totalv = PETSC_TRUE;
1372:         else totalv = PETSC_FALSE;
1373:       } else {
1374:         if (!nx) totalv = PETSC_TRUE;
1375:         else     totalv = PETSC_FALSE;
1376:       }
1377:       MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));

1379: #if defined(PETSC_USE_64BIT_INDICES)
1380:       if (cando && (yin->map->N < PETSC_MPI_INT_MAX)) {
1381: #else
1382:       if (cando) {
1383: #endif
1384:         MPI_Comm_size(PetscObjectComm((PetscObject)ctx),&size);
1385:         PetscMalloc3(1,&sto,size,&count,size,&displx);
1386:         range = xin->map->range;
1387:         for (i=0; i<size; i++) {
1388:           PetscMPIIntCast(range[i+1] - range[i],count+i);
1389:           PetscMPIIntCast(range[i],displx+i);
1390:         }
1391:         sto->count        = count;
1392:         sto->displx       = displx;
1393:         sto->work1        = 0;
1394:         sto->work2        = 0;
1395:         sto->type         = VEC_SCATTER_MPI_TOONE;
1396:         ctx->todata       = (void*)sto;
1397:         ctx->fromdata     = 0;
1398:         ctx->ops->begin   = VecScatterBegin_MPI_ToOne;
1399:         ctx->ops->end     = 0;
1400:         ctx->ops->destroy = VecScatterDestroy_MPI_ToAll;
1401:         ctx->ops->copy    = VecScatterCopy_MPI_ToAll;
1402:         ctx->ops->view    = VecScatterView_MPI_ToAll;
1403:         PetscInfo(xin,"Special case: processor zero gets entire parallel vector, rest get none\n");
1404:         goto functionend;
1405:       }
1406:     } else {
1407:       MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1408:     }

1410:     PetscObjectTypeCompare((PetscObject)ix,ISBLOCK,&ixblock);
1411:     PetscObjectTypeCompare((PetscObject)iy,ISBLOCK,&iyblock);
1412:     PetscObjectTypeCompare((PetscObject)iy,ISSTRIDE,&iystride);
1413:     if (ixblock) {
1414:       /* special case block to block */
1415:       if (iyblock) {
1416:         PetscInt       nx,ny,bsx,bsy;
1417:         const PetscInt *idx,*idy;
1418:         ISGetBlockSize(iy,&bsy);
1419:         ISGetBlockSize(ix,&bsx);
1420:         if (bsx == bsy && VecScatterOptimizedBS(bsx)) {
1421:           ISBlockGetLocalSize(ix,&nx);
1422:           ISBlockGetIndices(ix,&idx);
1423:           ISBlockGetLocalSize(iy,&ny);
1424:           ISBlockGetIndices(iy,&idy);
1425:           if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1426:           VecScatterCreate_PtoS(nx,idx,ny,idy,xin,yin,bsx,ctx);
1427:           ISBlockRestoreIndices(ix,&idx);
1428:           ISBlockRestoreIndices(iy,&idy);
1429:           PetscInfo(xin,"Special case: blocked indices\n");
1430:           goto functionend;
1431:         }
1432:         /* special case block to stride */
1433:       } else if (iystride) {
1434:         PetscInt ystart,ystride,ysize,bsx;
1435:         ISStrideGetInfo(iy,&ystart,&ystride);
1436:         ISGetLocalSize(iy,&ysize);
1437:         ISGetBlockSize(ix,&bsx);
1438:         /* see if stride index set is equivalent to block index set */
1439:         if (VecScatterOptimizedBS(bsx) && ((ystart % bsx) == 0) && (ystride == 1) && ((ysize % bsx) == 0)) {
1440:           PetscInt       nx,il,*idy;
1441:           const PetscInt *idx;
1442:           ISBlockGetLocalSize(ix,&nx);
1443:           ISBlockGetIndices(ix,&idx);
1444:           if (ysize != bsx*nx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1445:           PetscMalloc1(nx,&idy);
1446:           if (nx) {
1447:             idy[0] = ystart/bsx;
1448:             for (il=1; il<nx; il++) idy[il] = idy[il-1] + 1;
1449:           }
1450:           VecScatterCreate_PtoS(nx,idx,nx,idy,xin,yin,bsx,ctx);
1451:           PetscFree(idy);
1452:           ISBlockRestoreIndices(ix,&idx);
1453:           PetscInfo(xin,"Special case: blocked indices to stride\n");
1454:           goto functionend;
1455:         }
1456:       }
1457:     }
1458:     /* left over general case */
1459:     {
1460:       PetscInt       nx,ny;
1461:       const PetscInt *idx,*idy;
1462:       ISGetLocalSize(ix,&nx);
1463:       ISGetIndices(ix,&idx);
1464:       ISGetLocalSize(iy,&ny);
1465:       ISGetIndices(iy,&idy);
1466:       if (nx != ny) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match (%d %d)",nx,ny);
1467:       VecScatterCreate_PtoS(nx,idx,ny,idy,xin,yin,1,ctx);
1468:       ISRestoreIndices(ix,&idx);
1469:       ISRestoreIndices(iy,&idy);
1470:       PetscInfo(xin,"General case: MPI to Seq\n");
1471:       goto functionend;
1472:     }
1473:   }
1474:   /* ---------------------------------------------------------------------------*/
1475:   if (xin_type == VEC_SEQ_ID && yin_type == VEC_MPI_ID) {
1476:     /* ===========================================================================================================
1477:           Check for special cases
1478:        ==========================================================================================================*/
1479:     /* special case local copy portion */
1480:     islocal = PETSC_FALSE;
1481:     if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1482:       PetscInt              nx,ny,to_first,to_step,from_step,start,end,from_first,min,max;
1483:       VecScatter_Seq_Stride *from = NULL,*to = NULL;

1485:       VecGetOwnershipRange(yin,&start,&end);
1486:       ISGetLocalSize(ix,&nx);
1487:       ISStrideGetInfo(ix,&from_first,&from_step);
1488:       ISGetLocalSize(iy,&ny);
1489:       ISStrideGetInfo(iy,&to_first,&to_step);
1490:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1491:       ISGetMinMax(iy,&min,&max);
1492:       if (min >= start && max < end) islocal = PETSC_TRUE;
1493:       else islocal = PETSC_FALSE;
1494:       MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)yin));
1495:       if (cando) {
1496:         PetscMalloc2(1,&to,1,&from);
1497:         to->n             = nx;
1498:         to->first         = to_first-start;
1499:         to->step          = to_step;
1500:         from->n           = nx;
1501:         from->first       = from_first;
1502:         from->step        = from_step;
1503:         to->type          = VEC_SCATTER_SEQ_STRIDE;
1504:         from->type        = VEC_SCATTER_SEQ_STRIDE;
1505:         ctx->todata       = (void*)to;
1506:         ctx->fromdata     = (void*)from;
1507:         ctx->ops->begin   = VecScatterBegin_SSToSS;
1508:         ctx->ops->end     = 0;
1509:         ctx->ops->destroy = VecScatterDestroy_SSToSS;
1510:         ctx->ops->copy    = VecScatterCopy_SSToSS;
1511:         ctx->ops->view    = VecScatterView_SSToSS;
1512:         PetscInfo(xin,"Special case: sequential stride to MPI stride\n");
1513:         goto functionend;
1514:       }
1515:     } else {
1516:       MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)yin));
1517:     }
1518:     /* special case block to stride */
1519:     if (ix_type == IS_BLOCK_ID && iy_type == IS_STRIDE_ID) {
1520:       PetscInt ystart,ystride,ysize,bsx;
1521:       ISStrideGetInfo(iy,&ystart,&ystride);
1522:       ISGetLocalSize(iy,&ysize);
1523:       ISGetBlockSize(ix,&bsx);
1524:       /* see if stride index set is equivalent to block index set */
1525:       if (VecScatterOptimizedBS(bsx) && ((ystart % bsx) == 0) && (ystride == 1) && ((ysize % bsx) == 0)) {
1526:         PetscInt       nx,il,*idy;
1527:         const PetscInt *idx;
1528:         ISBlockGetLocalSize(ix,&nx);
1529:         ISBlockGetIndices(ix,&idx);
1530:         if (ysize != bsx*nx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1531:         PetscMalloc1(nx,&idy);
1532:         if (nx) {
1533:           idy[0] = ystart/bsx;
1534:           for (il=1; il<nx; il++) idy[il] = idy[il-1] + 1;
1535:         }
1536:         VecScatterCreate_StoP(nx,idx,nx,idy,xin,yin,bsx,ctx);
1537:         PetscFree(idy);
1538:         ISBlockRestoreIndices(ix,&idx);
1539:         PetscInfo(xin,"Special case: Blocked indices to stride\n");
1540:         goto functionend;
1541:       }
1542:     }

1544:     /* general case */
1545:     {
1546:       PetscInt       nx,ny;
1547:       const PetscInt *idx,*idy;
1548:       ISGetLocalSize(ix,&nx);
1549:       ISGetIndices(ix,&idx);
1550:       ISGetLocalSize(iy,&ny);
1551:       ISGetIndices(iy,&idy);
1552:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1553:       VecScatterCreate_StoP(nx,idx,ny,idy,xin,yin,1,ctx);
1554:       ISRestoreIndices(ix,&idx);
1555:       ISRestoreIndices(iy,&idy);
1556:       PetscInfo(xin,"General case: Seq to MPI\n");
1557:       goto functionend;
1558:     }
1559:   }
1560:   /* ---------------------------------------------------------------------------*/
1561:   if (xin_type == VEC_MPI_ID && yin_type == VEC_MPI_ID) {
1562:     /* no special cases for now */
1563:     PetscInt       nx,ny;
1564:     const PetscInt *idx,*idy;
1565:     ISGetLocalSize(ix,&nx);
1566:     ISGetIndices(ix,&idx);
1567:     ISGetLocalSize(iy,&ny);
1568:     ISGetIndices(iy,&idy);
1569:     if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1570:     VecScatterCreate_PtoP(nx,idx,ny,idy,xin,yin,1,ctx);
1571:     ISRestoreIndices(ix,&idx);
1572:     ISRestoreIndices(iy,&idy);
1573:     PetscInfo(xin,"General case: MPI to MPI\n");
1574:     goto functionend;
1575:   }

1577:   functionend:
1578:   *newctx = ctx;
1579:   ISDestroy(&tix);
1580:   ISDestroy(&tiy);
1581:   VecScatterViewFromOptions(ctx,NULL,"-vecscatter_view");
1582:   return(0);
1583: }

1585: /* ------------------------------------------------------------------*/
1588: /*@
1589:    VecScatterGetMerged - Returns true if the scatter is completed in the VecScatterBegin()
1590:       and the VecScatterEnd() does nothing

1592:    Not Collective

1594:    Input Parameter:
1595: .   ctx - scatter context created with VecScatterCreate()

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

1600:    Level: developer

1602: .seealso: VecScatterCreate(), VecScatterEnd(), VecScatterBegin()
1603: @*/
1604: PetscErrorCode  VecScatterGetMerged(VecScatter ctx,PetscBool  *flg)
1605: {
1608:   *flg = ctx->beginandendtogether;
1609:   return(0);
1610: }

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

1618:    Neighbor-wise Collective on VecScatter and Vec

1620:    Input Parameters:
1621: +  inctx - scatter context generated by VecScatterCreate()
1622: .  x - the vector from which we scatter
1623: .  y - the vector to which we scatter
1624: .  addv - either ADD_VALUES or INSERT_VALUES, with INSERT_VALUES mode any location
1625:           not scattered to retains its old value; i.e. the vector is NOT first zeroed.
1626: -  mode - the scattering mode, usually SCATTER_FORWARD.  The available modes are:
1627:     SCATTER_FORWARD or SCATTER_REVERSE


1630:    Level: intermediate

1632:    Options Database: See VecScatterCreate()

1634:    Notes:
1635:    The vectors x and y need not be the same vectors used in the call
1636:    to VecScatterCreate(), but x must have the same parallel data layout
1637:    as that passed in as the x to VecScatterCreate(), similarly for the y.
1638:    Most likely they have been obtained from VecDuplicate().

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

1643:    If you use SCATTER_REVERSE the two arguments x and y should be reversed, from
1644:    the SCATTER_FORWARD.

1646:    y[iy[i]] = x[ix[i]], for i=0,...,ni-1

1648:    This scatter is far more general than the conventional
1649:    scatter, since it can be a gather or a scatter or a combination,
1650:    depending on the indices ix and iy.  If x is a parallel vector and y
1651:    is sequential, VecScatterBegin() can serve to gather values to a
1652:    single processor.  Similarly, if y is parallel and x sequential, the
1653:    routine can scatter from one processor to many processors.

1655:    Concepts: scatter^between vectors
1656:    Concepts: gather^between vectors

1658: .seealso: VecScatterCreate(), VecScatterEnd()
1659: @*/
1660: PetscErrorCode  VecScatterBegin(VecScatter inctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
1661: {
1663: #if defined(PETSC_USE_DEBUG)
1664:   PetscInt       to_n,from_n;
1665: #endif
1670:   if (inctx->inuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE," Scatter ctx already in use");

1672: #if defined(PETSC_USE_DEBUG)
1673:   /*
1674:      Error checking to make sure these vectors match the vectors used
1675:    to create the vector scatter context. -1 in the from_n and to_n indicate the
1676:    vector lengths are unknown (for example with mapped scatters) and thus
1677:    no error checking is performed.
1678:   */
1679:   if (inctx->from_n >= 0 && inctx->to_n >= 0) {
1680:     VecGetLocalSize(x,&from_n);
1681:     VecGetLocalSize(y,&to_n);
1682:     if (mode & SCATTER_REVERSE) {
1683:       if (to_n != inctx->from_n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter reverse and vector to != ctx from size)",to_n,inctx->from_n);
1684:       if (from_n != inctx->to_n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter reverse and vector from != ctx to size)",from_n,inctx->to_n);
1685:     } else {
1686:       if (to_n != inctx->to_n)     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter forward and vector to != ctx to size)",to_n,inctx->to_n);
1687:       if (from_n != inctx->from_n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter forward and vector from != ctx from size)",from_n,inctx->from_n);
1688:     }
1689:   }
1690: #endif

1692:   inctx->inuse = PETSC_TRUE;
1693:   PetscLogEventBarrierBegin(VEC_ScatterBarrier,0,0,0,0,PetscObjectComm((PetscObject)inctx));
1694:   (*inctx->ops->begin)(inctx,x,y,addv,mode);
1695:   if (inctx->beginandendtogether && inctx->ops->end) {
1696:     inctx->inuse = PETSC_FALSE;
1697:     (*inctx->ops->end)(inctx,x,y,addv,mode);
1698:   }
1699:   PetscLogEventBarrierEnd(VEC_ScatterBarrier,0,0,0,0,PetscObjectComm((PetscObject)inctx));
1700:   return(0);
1701: }

1703: /* --------------------------------------------------------------------*/
1706: /*@
1707:    VecScatterEnd - Ends a generalized scatter from one vector to another.  Call
1708:    after first calling VecScatterBegin().

1710:    Neighbor-wise Collective on VecScatter and Vec

1712:    Input Parameters:
1713: +  ctx - scatter context generated by VecScatterCreate()
1714: .  x - the vector from which we scatter
1715: .  y - the vector to which we scatter
1716: .  addv - either ADD_VALUES or INSERT_VALUES.
1717: -  mode - the scattering mode, usually SCATTER_FORWARD.  The available modes are:
1718:      SCATTER_FORWARD, SCATTER_REVERSE

1720:    Level: intermediate

1722:    Notes:
1723:    If you use SCATTER_REVERSE the arguments x and y should be reversed, from the SCATTER_FORWARD.

1725:    y[iy[i]] = x[ix[i]], for i=0,...,ni-1

1727: .seealso: VecScatterBegin(), VecScatterCreate()
1728: @*/
1729: PetscErrorCode  VecScatterEnd(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
1730: {

1737:   ctx->inuse = PETSC_FALSE;
1738:   if (!ctx->ops->end) return(0);
1739:   if (!ctx->beginandendtogether) {
1740:     PetscLogEventBegin(VEC_ScatterEnd,ctx,x,y,0);
1741:     (*(ctx)->ops->end)(ctx,x,y,addv,mode);
1742:     PetscLogEventEnd(VEC_ScatterEnd,ctx,x,y,0);
1743:   }
1744:   return(0);
1745: }

1749: /*@C
1750:    VecScatterDestroy - Destroys a scatter context created by
1751:    VecScatterCreate().

1753:    Collective on VecScatter

1755:    Input Parameter:
1756: .  ctx - the scatter context

1758:    Level: intermediate

1760: .seealso: VecScatterCreate(), VecScatterCopy()
1761: @*/
1762: PetscErrorCode  VecScatterDestroy(VecScatter *ctx)
1763: {

1767:   if (!*ctx) return(0);
1769:   if ((*ctx)->inuse) SETERRQ(((PetscObject)(*ctx))->comm,PETSC_ERR_ARG_WRONGSTATE,"Scatter context is in use");
1770:   if (--((PetscObject)(*ctx))->refct > 0) {*ctx = 0; return(0);}

1772:   /* if memory was published with SAWs then destroy it */
1773:   PetscObjectSAWsViewOff((PetscObject)(*ctx));
1774:   (*(*ctx)->ops->destroy)(*ctx);
1775: #if defined(PETSC_HAVE_CUSP)
1776:   VecScatterCUSPIndicesDestroy((PetscCUSPIndices*)&((*ctx)->spptr));
1777: #endif
1778:   PetscHeaderDestroy(ctx);
1779:   return(0);
1780: }

1784: /*@
1785:    VecScatterCopy - Makes a copy of a scatter context.

1787:    Collective on VecScatter

1789:    Input Parameter:
1790: .  sctx - the scatter context

1792:    Output Parameter:
1793: .  ctx - the context copy

1795:    Level: advanced

1797: .seealso: VecScatterCreate(), VecScatterDestroy()
1798: @*/
1799: PetscErrorCode  VecScatterCopy(VecScatter sctx,VecScatter *ctx)
1800: {

1806:   if (!sctx->ops->copy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot copy this type");
1807:   PetscHeaderCreate(*ctx,VEC_SCATTER_CLASSID,"VecScatter","VecScatter","Vec",PetscObjectComm((PetscObject)sctx),VecScatterDestroy,VecScatterView);
1808:   (*ctx)->to_n   = sctx->to_n;
1809:   (*ctx)->from_n = sctx->from_n;
1810:   (*sctx->ops->copy)(sctx,*ctx);
1811:   return(0);
1812: }


1815: /* ------------------------------------------------------------------*/
1818: /*@
1819:    VecScatterView - Views a vector scatter context.

1821:    Collective on VecScatter

1823:    Input Parameters:
1824: +  ctx - the scatter context
1825: -  viewer - the viewer for displaying the context

1827:    Level: intermediate

1829: @*/
1830: PetscErrorCode  VecScatterView(VecScatter ctx,PetscViewer viewer)
1831: {

1836:   if (!viewer) {
1837:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ctx),&viewer);
1838:   }
1840:   if (ctx->ops->view) {
1841:     (*ctx->ops->view)(ctx,viewer);
1842:   }
1843:   return(0);
1844: }

1848: /*@C
1849:    VecScatterRemap - Remaps the "from" and "to" indices in a
1850:    vector scatter context. FOR EXPERTS ONLY!

1852:    Collective on VecScatter

1854:    Input Parameters:
1855: +  scat - vector scatter context
1856: .  from - remapping for "from" indices (may be NULL)
1857: -  to   - remapping for "to" indices (may be NULL)

1859:    Level: developer

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

1865:           In the sequential case the todata is the indices where the
1866:           data is put and the fromdata is where it is taken from.
1867:           This is backwards from the paralllel case! CRY! CRY! CRY!

1869: @*/
1870: PetscErrorCode  VecScatterRemap(VecScatter scat,PetscInt *rto,PetscInt *rfrom)
1871: {
1872:   VecScatter_Seq_General *to,*from;
1873:   VecScatter_MPI_General *mto;
1874:   PetscInt               i;


1881:   from = (VecScatter_Seq_General*)scat->fromdata;
1882:   mto  = (VecScatter_MPI_General*)scat->todata;

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

1886:   if (rto) {
1887:     if (mto->type == VEC_SCATTER_MPI_GENERAL) {
1888:       /* handle off processor parts */
1889:       for (i=0; i<mto->starts[mto->n]; i++) mto->indices[i] = rto[mto->indices[i]];

1891:       /* handle local part */
1892:       to = &mto->local;
1893:       for (i=0; i<to->n; i++) to->vslots[i] = rto[to->vslots[i]];
1894:     } else if (from->type == VEC_SCATTER_SEQ_GENERAL) {
1895:       for (i=0; i<from->n; i++) from->vslots[i] = rto[from->vslots[i]];
1896:     } else if (from->type == VEC_SCATTER_SEQ_STRIDE) {
1897:       VecScatter_Seq_Stride *sto = (VecScatter_Seq_Stride*)from;

1899:       /* if the remapping is the identity and stride is identity then skip remap */
1900:       if (sto->step == 1 && sto->first == 0) {
1901:         for (i=0; i<sto->n; i++) {
1902:           if (rto[i] != i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1903:         }
1904:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1905:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1906:   }

1908:   if (rfrom) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unable to remap the FROM in scatters yet");

1910:   /*
1911:      Mark then vector lengths as unknown because we do not know the
1912:    lengths of the remapped vectors
1913:   */
1914:   scat->from_n = -1;
1915:   scat->to_n   = -1;
1916:   return(0);
1917: }

1921: /*
1922:  VecScatterGetTypes_Private - Returns the scatter types.

1924:  scatter - The scatter.
1925:  from    - Upon exit this contains the type of the from scatter.
1926:  to      - Upon exit this contains the type of the to scatter.
1927: */
1928: PetscErrorCode VecScatterGetTypes_Private(VecScatter scatter,VecScatterType *from,VecScatterType *to)
1929: {
1930:   VecScatter_Common* fromdata = (VecScatter_Common*)scatter->fromdata;
1931:   VecScatter_Common* todata   = (VecScatter_Common*)scatter->todata;

1934:   *from = fromdata->type;
1935:   *to = todata->type;
1936:   return(0);
1937: }


1942: /*
1943:   VecScatterIsSequential_Private - Returns true if the scatter is sequential.

1945:   scatter - The scatter.
1946:   flag    - Upon exit flag is true if the scatter is of type VecScatter_Seq_General 
1947:             or VecScatter_Seq_Stride; otherwise flag is false.
1948: */
1949: PetscErrorCode VecScatterIsSequential_Private(VecScatter_Common *scatter,PetscBool *flag)
1950: {
1951:   VecScatterType scatterType = scatter->type;

1954:   if (scatterType == VEC_SCATTER_SEQ_GENERAL || scatterType == VEC_SCATTER_SEQ_STRIDE) {
1955:     *flag = PETSC_TRUE;
1956:   } else {
1957:     *flag = PETSC_FALSE;
1958:   }
1959:   return(0);
1960: }