Actual source code: vscat.c

petsc-3.5.2 2014-09-08
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:   VecGetArray(y,&yv);
 56:   if (x != y) {VecGetArray(x,&xv);}

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

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

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

136: PetscErrorCode VecScatterView_MPI_ToAll(VecScatter in,PetscViewer viewer)
137: {
139:   PetscBool      isascii;

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

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

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

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

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

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

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

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

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

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

249: PetscErrorCode VecScatterDestroy_SGToSG(VecScatter ctx)
250: {

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

261: PetscErrorCode VecScatterDestroy_SGToSS(VecScatter ctx)
262: {

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

273: PetscErrorCode VecScatterDestroy_SSToSG(VecScatter ctx)
274: {

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

285: PetscErrorCode VecScatterDestroy_SSToSS(VecScatter ctx)
286: {

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

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

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

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

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

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

357:   VecGetArray(x,&xv);
358:   if (x != y) {
359:     VecGetArray(y,&yv);
360:   } else yv = xv;

362:   if (mode & SCATTER_REVERSE) {
363:     gen_to   = (VecScatter_Seq_General*)ctx->fromdata;
364:     gen_from = (VecScatter_Seq_General*)ctx->todata;
365:   }
366:   fslots = gen_from->vslots;
367:   tslots = gen_to->vslots;

369:   if (addv == INSERT_VALUES) {
370:     for (i=0; i<n; i++) yv[tslots[i]]  = xv[fslots[i]];
371:   } else if (addv == ADD_VALUES) {
372:     for (i=0; i<n; i++) yv[tslots[i]] += xv[fslots[i]];
373: #if !defined(PETSC_USE_COMPLEX)
374:   } else if (addv == MAX_VALUES) {
375:     for (i=0; i<n; i++) yv[tslots[i]] = PetscMax(yv[tslots[i]],xv[fslots[i]]);
376: #endif
377:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
378:   VecRestoreArray(x,&xv);
379:   if (x != y) {VecRestoreArray(y,&yv);}
380:   return(0);
381: }

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

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

412:   VecGetArray(x,&xv);
413:   if (x != y) {
414:     VecGetArray(y,&yv);
415:   } else yv = xv;
416:   if (mode & SCATTER_REVERSE) {
417:     xv += first;
418:     if (addv == INSERT_VALUES) {
419:       for (i=0; i<n; i++) yv[fslots[i]] = xv[i];
420:     } else if (addv == ADD_VALUES) {
421:       for (i=0; i<n; i++) yv[fslots[i]] += xv[i];
422: #if !defined(PETSC_USE_COMPLEX)
423:     } else if (addv == MAX_VALUES) {
424:       for (i=0; i<n; i++) yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[i]);
425: #endif
426:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
427:   } else {
428:     yv += first;
429:     if (addv == INSERT_VALUES) {
430:       for (i=0; i<n; i++) yv[i] = xv[fslots[i]];
431:     } else if (addv == ADD_VALUES) {
432:       for (i=0; i<n; i++) yv[i] += xv[fslots[i]];
433: #if !defined(PETSC_USE_COMPLEX)
434:     } else if (addv == MAX_VALUES) {
435:       for (i=0; i<n; i++) yv[i] = PetscMax(yv[i],xv[fslots[i]]);
436: #endif
437:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
438:   }
439:   VecRestoreArray(x,&xv);
440:   if (x != y) {VecRestoreArray(y,&yv);}
441:   return(0);
442: }

444: /*
445:    Scatter: sequential general to sequential stride
446: */
449: PetscErrorCode VecScatterBegin_SGToSS(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
450: {
451:   VecScatter_Seq_Stride  *gen_to   = (VecScatter_Seq_Stride*)ctx->todata;
452:   VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
453:   PetscInt               i,n = gen_from->n,*fslots = gen_from->vslots;
454:   PetscErrorCode         ierr;
455:   PetscInt               first = gen_to->first,step = gen_to->step;
456:   PetscScalar            *xv,*yv;

459: #if defined(PETSC_HAVE_CUSP)
460:   if (x->valid_GPU_array == PETSC_CUSP_GPU) {
461:     /* create the scatter indices if not done already */
462:     if (!ctx->spptr) {
463:       PetscInt tofirst = first,tostep = step,fromfirst = 0,fromstep = 0;
464:       PetscInt * tslots = 0;
465:       VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
466:     }
467:     /* next do the scatter */
468:     VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
469:     return(0);
470:   }
471: #endif

473:   VecGetArray(x,&xv);
474:   if (x != y) {
475:     VecGetArray(y,&yv);
476:   } else yv = xv;

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

504: /*
505:     Scatter: sequential stride 1 to sequential general
506: */
509: PetscErrorCode VecScatterBegin_SSToSG_Stride1(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
510: {
511:   VecScatter_Seq_Stride  *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
512:   VecScatter_Seq_General *gen_to   = (VecScatter_Seq_General*)ctx->todata;
513:   PetscInt               i,n = gen_from->n,*fslots = gen_to->vslots;
514:   PetscErrorCode         ierr;
515:   PetscInt               first = gen_from->first;
516:   PetscScalar            *xv,*yv;

519: #if defined(PETSC_HAVE_CUSP)
520:   if (x->valid_GPU_array == PETSC_CUSP_GPU) {
521:     /* create the scatter indices if not done already */
522:     if (!ctx->spptr) {
523:       PetscInt tofirst = 0,tostep = 0,fromfirst = first,fromstep = 1;
524:       PetscInt * tslots = 0;
525:       VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
526:     }
527:     /* next do the scatter */
528:     VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
529:     return(0);
530:   }
531: #endif

533:   VecGetArray(x,&xv);
534:   if (x != y) {
535:     VecGetArray(y,&yv);
536:   } else yv = xv;

538:   if (mode & SCATTER_REVERSE) {
539:     yv += first;
540:     if (addv == INSERT_VALUES) {
541:       for (i=0; i<n; i++) yv[i] = xv[fslots[i]];
542:     } else  if (addv == ADD_VALUES) {
543:       for (i=0; i<n; i++) yv[i] += xv[fslots[i]];
544: #if !defined(PETSC_USE_COMPLEX)
545:     } else  if (addv == MAX_VALUES) {
546:       for (i=0; i<n; i++) yv[i] = PetscMax(yv[i],xv[fslots[i]]);
547: #endif
548:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
549:   } else {
550:     xv += first;
551:     if (addv == INSERT_VALUES) {
552:       for (i=0; i<n; i++) yv[fslots[i]] = xv[i];
553:     } else  if (addv == ADD_VALUES) {
554:       for (i=0; i<n; i++) yv[fslots[i]] += xv[i];
555: #if !defined(PETSC_USE_COMPLEX)
556:     } else  if (addv == MAX_VALUES) {
557:       for (i=0; i<n; i++) yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[i]);
558: #endif
559:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
560:   }
561:   VecRestoreArray(x,&xv);
562:   if (x != y) {VecRestoreArray(y,&yv);}
563:   return(0);
564: }

568: /*
569:    Scatter: sequential stride to sequential general
570: */
571: PetscErrorCode VecScatterBegin_SSToSG(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
572: {
573:   VecScatter_Seq_Stride  *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
574:   VecScatter_Seq_General *gen_to   = (VecScatter_Seq_General*)ctx->todata;
575:   PetscInt               i,n = gen_from->n,*fslots = gen_to->vslots;
576:   PetscErrorCode         ierr;
577:   PetscInt               first = gen_from->first,step = gen_from->step;
578:   PetscScalar            *xv,*yv;

581: #if defined(PETSC_HAVE_CUSP)
582:   if (x->valid_GPU_array == PETSC_CUSP_GPU) {
583:     /* create the scatter indices if not done already */
584:     if (!ctx->spptr) {
585:       PetscInt tofirst = 0,tostep = 0,fromfirst = first,fromstep = step;
586:       PetscInt * tslots = 0;
587:       VecScatterCUSPIndicesCreate_StoS(n,tofirst,fromfirst,tostep,fromstep,fslots,tslots,(PetscCUSPIndices*)&(ctx->spptr));
588:     }
589:     /* next do the scatter */
590:     VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
591:     return(0);
592:   }
593: #endif

595:   VecGetArray(x,&xv);
596:   if (x != y) {
597:     VecGetArray(y,&yv);
598:   } else yv = xv;

600:   if (mode & SCATTER_REVERSE) {
601:     if (addv == INSERT_VALUES) {
602:       for (i=0; i<n; i++) yv[first + i*step] = xv[fslots[i]];
603:     } else if (addv == ADD_VALUES) {
604:       for (i=0; i<n; i++) yv[first + i*step] += xv[fslots[i]];
605: #if !defined(PETSC_USE_COMPLEX)
606:     } else if (addv == MAX_VALUES) {
607:       for (i=0; i<n; i++) yv[first + i*step] = PetscMax(yv[first + i*step],xv[fslots[i]]);
608: #endif
609:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
610:   } else {
611:     if (addv == INSERT_VALUES) {
612:       for (i=0; i<n; i++) yv[fslots[i]] = xv[first + i*step];
613:     } else if (addv == ADD_VALUES) {
614:       for (i=0; i<n; i++) yv[fslots[i]] += xv[first + i*step];
615: #if !defined(PETSC_USE_COMPLEX)
616:     } else if (addv == MAX_VALUES) {
617:       for (i=0; i<n; i++) yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[first + i*step]);
618: #endif
619:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
620:   }
621:   VecRestoreArray(x,&xv);
622:   if (x != y) {VecRestoreArray(y,&yv);}
623:   return(0);
624: }

628: PetscErrorCode VecScatterView_SSToSG(VecScatter in,PetscViewer viewer)
629: {
630:   PetscErrorCode         ierr;
631:   VecScatter_Seq_Stride  *in_from = (VecScatter_Seq_Stride*)in->todata;
632:   VecScatter_Seq_General *in_to   = (VecScatter_Seq_General*)in->fromdata;
633:   PetscInt               i;
634:   PetscBool              isascii;

637:   PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
638:   if (isascii) {
639:     PetscViewerASCIIPrintf(viewer,"Sequential stride to general scatter\n");
640:     for (i=0; i<in_to->n; i++) {
641:       PetscViewerASCIIPrintf(viewer,"%D to %D\n",in_from->first + in_from->step*i,in_to->vslots[i]);
642:     }
643:   }
644:   return(0);
645: }
646: /*
647:      Scatter: sequential stride to sequential stride
648: */
651: PetscErrorCode VecScatterBegin_SSToSS(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
652: {
653:   VecScatter_Seq_Stride *gen_to   = (VecScatter_Seq_Stride*)ctx->todata;
654:   VecScatter_Seq_Stride *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
655:   PetscInt              i,n = gen_from->n,to_first = gen_to->first,to_step = gen_to->step;
656:   PetscErrorCode        ierr;
657:   PetscInt              from_first = gen_from->first,from_step = gen_from->step;
658:   PetscScalar           *xv,*yv;

661: #if defined(PETSC_HAVE_CUSP)
662:   if (x->valid_GPU_array == PETSC_CUSP_GPU) {
663:     /* create the scatter indices if not done already */
664:     if (!ctx->spptr) {
665:       PetscInt *tslots = 0,*fslots = 0;
666:       VecScatterCUSPIndicesCreate_StoS(n,to_first,from_first,to_step,from_step,tslots,fslots,(PetscCUSPIndices*)&(ctx->spptr));
667:     }
668:     /* next do the scatter */
669:     VecScatterCUSP_StoS(x,y,(PetscCUSPIndices)ctx->spptr,addv,mode);
670:     return(0);
671:   }
672: #endif

674:   VecGetArrayRead(x,(const PetscScalar **)&xv);
675:   if (x != y) {
676:     VecGetArray(y,&yv);
677:   } else yv = xv;

679:   if (mode & SCATTER_REVERSE) {
680:     from_first = gen_to->first;
681:     to_first   = gen_from->first;
682:     from_step  = gen_to->step;
683:     to_step    = gen_from->step;
684:   }

686:   if (addv == INSERT_VALUES) {
687:     if (to_step == 1 && from_step == 1) {
688:       PetscMemcpy(yv+to_first,xv+from_first,n*sizeof(PetscScalar));
689:     } else  {
690:       for (i=0; i<n; i++) yv[to_first + i*to_step] = xv[from_first+i*from_step];
691:     }
692:   } else if (addv == ADD_VALUES) {
693:     if (to_step == 1 && from_step == 1) {
694:       yv += to_first; xv += from_first;
695:       for (i=0; i<n; i++) yv[i] += xv[i];
696:     } else {
697:       for (i=0; i<n; i++) yv[to_first + i*to_step] += xv[from_first+i*from_step];
698:     }
699: #if !defined(PETSC_USE_COMPLEX)
700:   } else if (addv == MAX_VALUES) {
701:     if (to_step == 1 && from_step == 1) {
702:       yv += to_first; xv += from_first;
703:       for (i=0; i<n; i++) yv[i] = PetscMax(yv[i],xv[i]);
704:     } else {
705:       for (i=0; i<n; i++) yv[to_first + i*to_step] = PetscMax(yv[to_first + i*to_step],xv[from_first+i*from_step]);
706:     }
707: #endif
708:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
709:   VecRestoreArrayRead(x,(const PetscScalar**)&xv);
710:   if (x != y) {VecRestoreArray(y,&yv);}
711:   return(0);
712: }

714: /* --------------------------------------------------------------------------*/


719: PetscErrorCode VecScatterCopy_SGToSG(VecScatter in,VecScatter out)
720: {
721:   PetscErrorCode         ierr;
722:   VecScatter_Seq_General *in_to   = (VecScatter_Seq_General*)in->todata,*out_to = NULL;
723:   VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from = NULL;

726:   out->begin   = in->begin;
727:   out->end     = in->end;
728:   out->copy    = in->copy;
729:   out->destroy = in->destroy;
730:   out->view    = in->view;

732:   PetscMalloc2(1,&out_to,1,&out_from);
733:   PetscMalloc2(in_to->n,&out_to->vslots,in_from->n,&out_from->vslots);
734:   out_to->n                    = in_to->n;
735:   out_to->type                 = in_to->type;
736:   out_to->nonmatching_computed = PETSC_FALSE;
737:   out_to->slots_nonmatching    = 0;
738:   out_to->is_copy              = PETSC_FALSE;
739:   PetscMemcpy(out_to->vslots,in_to->vslots,(out_to->n)*sizeof(PetscInt));

741:   out_from->n                    = in_from->n;
742:   out_from->type                 = in_from->type;
743:   out_from->nonmatching_computed = PETSC_FALSE;
744:   out_from->slots_nonmatching    = 0;
745:   out_from->is_copy              = PETSC_FALSE;
746:   PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));

748:   out->todata   = (void*)out_to;
749:   out->fromdata = (void*)out_from;
750:   return(0);
751: }

755: PetscErrorCode VecScatterView_SGToSG(VecScatter in,PetscViewer viewer)
756: {
757:   PetscErrorCode         ierr;
758:   VecScatter_Seq_General *in_to   = (VecScatter_Seq_General*)in->todata;
759:   VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata;
760:   PetscInt               i;
761:   PetscBool              isascii;

764:   PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
765:   if (isascii) {
766:     PetscViewerASCIIPrintf(viewer,"Sequential general scatter\n");
767:     for (i=0; i<in_to->n; i++) {
768:       PetscViewerASCIIPrintf(viewer,"%D to %D\n",in_from->vslots[i],in_to->vslots[i]);
769:     }
770:   }
771:   return(0);
772: }


777: PetscErrorCode VecScatterCopy_SGToSS(VecScatter in,VecScatter out)
778: {
779:   PetscErrorCode         ierr;
780:   VecScatter_Seq_Stride  *in_to   = (VecScatter_Seq_Stride*)in->todata,*out_to = NULL;
781:   VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from = NULL;

784:   out->begin   = in->begin;
785:   out->end     = in->end;
786:   out->copy    = in->copy;
787:   out->destroy = in->destroy;
788:   out->view    = in->view;

790:   PetscMalloc2(1,&out_to,1,&out_from);
791:   PetscMalloc1(in_from->n,&out_from->vslots);
792:   out_to->n     = in_to->n;
793:   out_to->type  = in_to->type;
794:   out_to->first = in_to->first;
795:   out_to->step  = in_to->step;
796:   out_to->type  = in_to->type;

798:   out_from->n                    = in_from->n;
799:   out_from->type                 = in_from->type;
800:   out_from->nonmatching_computed = PETSC_FALSE;
801:   out_from->slots_nonmatching    = 0;
802:   out_from->is_copy              = PETSC_FALSE;
803:   PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));

805:   out->todata   = (void*)out_to;
806:   out->fromdata = (void*)out_from;
807:   return(0);
808: }

812: PetscErrorCode VecScatterView_SGToSS(VecScatter in,PetscViewer viewer)
813: {
814:   PetscErrorCode         ierr;
815:   VecScatter_Seq_Stride  *in_to   = (VecScatter_Seq_Stride*)in->todata;
816:   VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata;
817:   PetscInt               i;
818:   PetscBool              isascii;

821:   PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
822:   if (isascii) {
823:     PetscViewerASCIIPrintf(viewer,"Sequential general scatter to stride\n");
824:     for (i=0; i<in_to->n; i++) {
825:       PetscViewerASCIIPrintf(viewer,"%D to %D\n",in_from->vslots[i],in_to->first + in_to->step*i);
826:     }
827:   }
828:   return(0);
829: }

831: /* --------------------------------------------------------------------------*/
832: /*
833:     Scatter: parallel to sequential vector, sequential strides for both.
834: */
837: PetscErrorCode VecScatterCopy_SSToSS(VecScatter in,VecScatter out)
838: {
839:   VecScatter_Seq_Stride *in_to   = (VecScatter_Seq_Stride*)in->todata,*out_to = NULL;
840:   VecScatter_Seq_Stride *in_from = (VecScatter_Seq_Stride*)in->fromdata,*out_from = NULL;
841:   PetscErrorCode        ierr;

844:   out->begin   = in->begin;
845:   out->end     = in->end;
846:   out->copy    = in->copy;
847:   out->destroy = in->destroy;
848:   out->view    = in->view;

850:   PetscMalloc2(1,&out_to,1,&out_from);
851:   out_to->n       = in_to->n;
852:   out_to->type    = in_to->type;
853:   out_to->first   = in_to->first;
854:   out_to->step    = in_to->step;
855:   out_to->type    = in_to->type;
856:   out_from->n     = in_from->n;
857:   out_from->type  = in_from->type;
858:   out_from->first = in_from->first;
859:   out_from->step  = in_from->step;
860:   out_from->type  = in_from->type;
861:   out->todata     = (void*)out_to;
862:   out->fromdata   = (void*)out_from;
863:   return(0);
864: }

868: PetscErrorCode VecScatterView_SSToSS(VecScatter in,PetscViewer viewer)
869: {
870:   VecScatter_Seq_Stride *in_to   = (VecScatter_Seq_Stride*)in->todata;
871:   VecScatter_Seq_Stride *in_from = (VecScatter_Seq_Stride*)in->fromdata;
872:   PetscErrorCode        ierr;
873:   PetscBool             isascii;

876:   PetscObjectTypeCompare((PetscObject)in,PETSCVIEWERASCII,&isascii);
877:   if (isascii) {
878:     PetscViewerASCIIPrintf(viewer,"Sequential stride count %D start %D step to start %D stride %D\n",in_to->n,in_to->first,in_to->step,in_from->first,in_from->step);
879:   }
880:   return(0);
881: }


884: extern PetscErrorCode VecScatterCreate_PtoS(PetscInt,const PetscInt*,PetscInt,const PetscInt*,Vec,Vec,PetscInt,VecScatter);
885: extern PetscErrorCode VecScatterCreate_PtoP(PetscInt,const PetscInt*,PetscInt,const PetscInt*,Vec,Vec,PetscInt,VecScatter);
886: extern PetscErrorCode VecScatterCreate_StoP(PetscInt,const PetscInt*,PetscInt,const PetscInt*,Vec,Vec,PetscInt,VecScatter);

888: /* =======================================================================*/
889: #define VEC_SEQ_ID 0
890: #define VEC_MPI_ID 1
891: #define IS_GENERAL_ID 0
892: #define IS_STRIDE_ID  1
893: #define IS_BLOCK_ID   2

895: /*
896:    Blocksizes we have optimized scatters for
897: */

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

901: PetscErrorCode  VecScatterCreateEmpty(MPI_Comm comm,VecScatter *newctx)
902: {
903:   VecScatter     ctx;

907:   PetscHeaderCreate(ctx,_p_VecScatter,int,VEC_SCATTER_CLASSID,"VecScatter","VecScatter","Vec",comm,VecScatterDestroy,VecScatterView);
908:   ctx->inuse               = PETSC_FALSE;
909:   ctx->beginandendtogether = PETSC_FALSE;
910:   PetscOptionsGetBool(NULL,"-vecscatter_merge",&ctx->beginandendtogether,NULL);
911:   if (ctx->beginandendtogether) {
912:     PetscInfo(ctx,"Using combined (merged) vector scatter begin and end\n");
913:   }
914:   ctx->packtogether = PETSC_FALSE;
915:   PetscOptionsGetBool(NULL,"-vecscatter_packtogether",&ctx->packtogether,NULL);
916:   if (ctx->packtogether) {
917:     PetscInfo(ctx,"Pack all messages before sending\n");
918:   }
919:   *newctx = ctx;
920:   return(0);
921: }

925: /*@C
926:    VecScatterCreate - Creates a vector scatter context.

928:    Collective on Vec

930:    Input Parameters:
931: +  xin - a vector that defines the shape (parallel data layout of the vector)
932:          of vectors from which we scatter
933: .  yin - a vector that defines the shape (parallel data layout of the vector)
934:          of vectors to which we scatter
935: .  ix - the indices of xin to scatter (if NULL scatters all values)
936: -  iy - the indices of yin to hold results (if NULL fills entire vector yin)

938:    Output Parameter:
939: .  newctx - location to store the new scatter context

941:    Options Database Keys: (uses regular MPI_Sends by default)
942: +  -vecscatter_view         - Prints detail of communications
943: .  -vecscatter_view ::ascii_info    - Print less details about communication
944: .  -vecscatter_ssend        - Uses MPI_Ssend_init() instead of MPI_Send_init()
945: .  -vecscatter_rsend           - use ready receiver mode for MPI sends
946: .  -vecscatter_merge        - VecScatterBegin() handles all of the communication, VecScatterEnd() is a nop
947:                               eliminates the chance for overlap of computation and communication
948: .  -vecscatter_sendfirst    - Posts sends before receives
949: .  -vecscatter_packtogether - Pack all messages before sending, receive all messages before unpacking
950: .  -vecscatter_alltoall     - Uses MPI all to all communication for scatter
951: .  -vecscatter_window       - Use MPI 2 window operations to move data
952: .  -vecscatter_nopack       - Avoid packing to work vector when possible (if used with -vecscatter_alltoall then will use MPI_Alltoallw()
953: -  -vecscatter_reproduce    - insure that the order of the communications are done the same for each scatter, this under certain circumstances
954:                               will make the results of scatters deterministic when otherwise they are not (it may be slower also).

956: $
957: $                                                                                    --When packing is used--
958: $                               MPI Datatypes (no packing)  sendfirst   merge        packtogether  persistent*
959: $                                _nopack                   _sendfirst    _merge      _packtogether                -vecscatter_
960: $ ----------------------------------------------------------------------------------------------------------------------------
961: $    Message passing    Send       p                           X            X           X         always
962: $                      Ssend       p                           X            X           X         always          _ssend
963: $                      Rsend       p                        nonsense        X           X         always          _rsend
964: $    AlltoAll  v or w              X                        nonsense     always         X         nonsense        _alltoall
965: $    MPI_Win                       p                        nonsense        p           p         nonsense        _window
966: $
967: $   Since persistent sends and receives require a constant memory address they can only be used when data is packed into the work vector
968: $   because the in and out array may be different for each call to VecScatterBegin/End().
969: $
970: $    p indicates possible, but not implemented. X indicates implemented
971: $

973:     Level: intermediate

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

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

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

991:    Concepts: scatter^between vectors
992:    Concepts: gather^between vectors

994: .seealso: VecScatterDestroy(), VecScatterCreateToAll(), VecScatterCreateToZero()
995: @*/
996: PetscErrorCode  VecScatterCreate(Vec xin,IS ix,Vec yin,IS iy,VecScatter *newctx)
997: {
998:   VecScatter        ctx;
999:   PetscErrorCode    ierr;
1000:   PetscMPIInt       size;
1001:   PetscInt          xin_type = VEC_SEQ_ID,yin_type = VEC_SEQ_ID,*range;
1002:   PetscInt          ix_type  = IS_GENERAL_ID,iy_type = IS_GENERAL_ID;
1003:   MPI_Comm          comm,ycomm;
1004:   PetscBool         totalv,ixblock,iyblock,iystride,islocal,cando,flag;
1005:   IS                tix = 0,tiy = 0;

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

1010:   /*
1011:       Determine if the vectors are "parallel", ie. it shares a comm with other processors, or
1012:       sequential (it does not share a comm). The difference is that parallel vectors treat the
1013:       index set as providing indices in the global parallel numbering of the vector, with
1014:       sequential vectors treat the index set as providing indices in the local sequential
1015:       numbering
1016:   */
1017:   PetscObjectGetComm((PetscObject)xin,&comm);
1018:   MPI_Comm_size(comm,&size);
1019:   if (size > 1) xin_type = VEC_MPI_ID;

1021:   PetscObjectGetComm((PetscObject)yin,&ycomm);
1022:   MPI_Comm_size(ycomm,&size);
1023:   if (size > 1) {comm = ycomm; yin_type = VEC_MPI_ID;}

1025:   /* generate the Scatter context */
1026:   PetscHeaderCreate(ctx,_p_VecScatter,int,VEC_SCATTER_CLASSID,"VecScatter","VecScatter","Vec",comm,VecScatterDestroy,VecScatterView);
1027:   ctx->inuse               = PETSC_FALSE;

1029:   ctx->beginandendtogether = PETSC_FALSE;
1030:   PetscOptionsGetBool(NULL,"-vecscatter_merge",&ctx->beginandendtogether,NULL);
1031:   if (ctx->beginandendtogether) {
1032:     PetscInfo(ctx,"Using combined (merged) vector scatter begin and end\n");
1033:   }
1034:   ctx->packtogether = PETSC_FALSE;
1035:   PetscOptionsGetBool(NULL,"-vecscatter_packtogether",&ctx->packtogether,NULL);
1036:   if (ctx->packtogether) {
1037:     PetscInfo(ctx,"Pack all messages before sending\n");
1038:   }

1040:   VecGetLocalSize(xin,&ctx->from_n);
1041:   VecGetLocalSize(yin,&ctx->to_n);

1043:   /*
1044:       if ix or iy is not included; assume just grabbing entire vector
1045:   */
1046:   if (!ix && xin_type == VEC_SEQ_ID) {
1047:     ISCreateStride(comm,ctx->from_n,0,1,&ix);
1048:     tix  = ix;
1049:   } else if (!ix && xin_type == VEC_MPI_ID) {
1050:     if (yin_type == VEC_MPI_ID) {
1051:       PetscInt ntmp, low;
1052:       VecGetLocalSize(xin,&ntmp);
1053:       VecGetOwnershipRange(xin,&low,NULL);
1054:       ISCreateStride(comm,ntmp,low,1,&ix);
1055:     } else {
1056:       PetscInt Ntmp;
1057:       VecGetSize(xin,&Ntmp);
1058:       ISCreateStride(comm,Ntmp,0,1,&ix);
1059:     }
1060:     tix = ix;
1061:   } else if (!ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ix not given, but not Seq or MPI vector");

1063:   if (!iy && yin_type == VEC_SEQ_ID) {
1064:     ISCreateStride(comm,ctx->to_n,0,1,&iy);
1065:     tiy  = iy;
1066:   } else if (!iy && yin_type == VEC_MPI_ID) {
1067:     if (xin_type == VEC_MPI_ID) {
1068:       PetscInt ntmp, low;
1069:       VecGetLocalSize(yin,&ntmp);
1070:       VecGetOwnershipRange(yin,&low,NULL);
1071:       ISCreateStride(comm,ntmp,low,1,&iy);
1072:     } else {
1073:       PetscInt Ntmp;
1074:       VecGetSize(yin,&Ntmp);
1075:       ISCreateStride(comm,Ntmp,0,1,&iy);
1076:     }
1077:     tiy = iy;
1078:   } else if (!iy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"iy not given, but not Seq or MPI vector");

1080:   /*
1081:      Determine types of index sets
1082:   */
1083:   PetscObjectTypeCompare((PetscObject)ix,ISBLOCK,&flag);
1084:   if (flag) ix_type = IS_BLOCK_ID;
1085:   PetscObjectTypeCompare((PetscObject)iy,ISBLOCK,&flag);
1086:   if (flag) iy_type = IS_BLOCK_ID;
1087:   PetscObjectTypeCompare((PetscObject)ix,ISSTRIDE,&flag);
1088:   if (flag) ix_type = IS_STRIDE_ID;
1089:   PetscObjectTypeCompare((PetscObject)iy,ISSTRIDE,&flag);
1090:   if (flag) iy_type = IS_STRIDE_ID;

1092:   /* ===========================================================================================================
1093:         Check for special cases
1094:      ==========================================================================================================*/
1095:   /* ---------------------------------------------------------------------------*/
1096:   if (xin_type == VEC_SEQ_ID && yin_type == VEC_SEQ_ID) {
1097:     if (ix_type == IS_GENERAL_ID && iy_type == IS_GENERAL_ID) {
1098:       PetscInt               nx,ny;
1099:       const PetscInt         *idx,*idy;
1100:       VecScatter_Seq_General *to = NULL,*from = NULL;

1102:       ISGetLocalSize(ix,&nx);
1103:       ISGetLocalSize(iy,&ny);
1104:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1105:       ISGetIndices(ix,&idx);
1106:       ISGetIndices(iy,&idy);
1107:       PetscMalloc2(1,&to,1,&from);
1108:       PetscMalloc2(nx,&to->vslots,nx,&from->vslots);
1109:       to->n = nx;
1110: #if defined(PETSC_USE_DEBUG)
1111:       VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1112: #endif
1113:       PetscMemcpy(to->vslots,idy,nx*sizeof(PetscInt));
1114:       from->n = nx;
1115: #if defined(PETSC_USE_DEBUG)
1116:       VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1117: #endif
1118:        PetscMemcpy(from->vslots,idx,nx*sizeof(PetscInt));
1119:       to->type      = VEC_SCATTER_SEQ_GENERAL;
1120:       from->type    = VEC_SCATTER_SEQ_GENERAL;
1121:       ctx->todata   = (void*)to;
1122:       ctx->fromdata = (void*)from;
1123:       ctx->begin    = VecScatterBegin_SGToSG;
1124:       ctx->end      = 0;
1125:       ctx->destroy  = VecScatterDestroy_SGToSG;
1126:       ctx->copy     = VecScatterCopy_SGToSG;
1127:       ctx->view     = VecScatterView_SGToSG;
1128:       PetscInfo(xin,"Special case: sequential vector general scatter\n");
1129:       goto functionend;
1130:     } else if (ix_type == IS_STRIDE_ID &&  iy_type == IS_STRIDE_ID) {
1131:       PetscInt              nx,ny,to_first,to_step,from_first,from_step;
1132:       VecScatter_Seq_Stride *from8 = NULL,*to8 = NULL;

1134:       ISGetLocalSize(ix,&nx);
1135:       ISGetLocalSize(iy,&ny);
1136:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1137:       ISStrideGetInfo(iy,&to_first,&to_step);
1138:       ISStrideGetInfo(ix,&from_first,&from_step);
1139:       PetscMalloc2(1,&to8,1,&from8);
1140:       to8->n        = nx;
1141:       to8->first    = to_first;
1142:       to8->step     = to_step;
1143:       from8->n      = nx;
1144:       from8->first  = from_first;
1145:       from8->step   = from_step;
1146:       to8->type     = VEC_SCATTER_SEQ_STRIDE;
1147:       from8->type   = VEC_SCATTER_SEQ_STRIDE;
1148:       ctx->todata   = (void*)to8;
1149:       ctx->fromdata = (void*)from8;
1150:       ctx->begin    = VecScatterBegin_SSToSS;
1151:       ctx->end      = 0;
1152:       ctx->destroy  = VecScatterDestroy_SSToSS;
1153:       ctx->copy     = VecScatterCopy_SSToSS;
1154:       ctx->view     = VecScatterView_SSToSS;
1155:       PetscInfo(xin,"Special case: sequential vector stride to stride\n");
1156:       goto functionend;
1157:     } else if (ix_type == IS_GENERAL_ID && iy_type == IS_STRIDE_ID) {
1158:       PetscInt               nx,ny,first,step;
1159:       const PetscInt         *idx;
1160:       VecScatter_Seq_General *from9 = NULL;
1161:       VecScatter_Seq_Stride  *to9   = NULL;

1163:       ISGetLocalSize(ix,&nx);
1164:       ISGetIndices(ix,&idx);
1165:       ISGetLocalSize(iy,&ny);
1166:       ISStrideGetInfo(iy,&first,&step);
1167:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1168:       PetscMalloc2(1,&to9,1,&from9);
1169:       PetscMalloc1(nx,&from9->vslots);
1170:       to9->n     = nx;
1171:       to9->first = first;
1172:       to9->step  = step;
1173:       from9->n   = nx;
1174: #if defined(PETSC_USE_DEBUG)
1175:       VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1176: #endif
1177:       PetscMemcpy(from9->vslots,idx,nx*sizeof(PetscInt));
1178:       ctx->todata    = (void*)to9; ctx->fromdata = (void*)from9;
1179:       if (step == 1)  ctx->begin = VecScatterBegin_SGToSS_Stride1;
1180:       else            ctx->begin = VecScatterBegin_SGToSS;
1181:       ctx->destroy   = VecScatterDestroy_SGToSS;
1182:       ctx->end       = 0;
1183:       ctx->copy      = VecScatterCopy_SGToSS;
1184:       ctx->view      = VecScatterView_SGToSS;
1185:       to9->type      = VEC_SCATTER_SEQ_STRIDE;
1186:       from9->type    = VEC_SCATTER_SEQ_GENERAL;
1187:       PetscInfo(xin,"Special case: sequential vector general to stride\n");
1188:       goto functionend;
1189:     } else if (ix_type == IS_STRIDE_ID && iy_type == IS_GENERAL_ID) {
1190:       PetscInt               nx,ny,first,step;
1191:       const PetscInt         *idy;
1192:       VecScatter_Seq_General *to10 = NULL;
1193:       VecScatter_Seq_Stride  *from10 = NULL;

1195:       ISGetLocalSize(ix,&nx);
1196:       ISGetIndices(iy,&idy);
1197:       ISGetLocalSize(iy,&ny);
1198:       ISStrideGetInfo(ix,&first,&step);
1199:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1200:       PetscMalloc2(1,&to10,1,&from10);
1201:       PetscMalloc1(nx,&to10->vslots);
1202:       from10->n     = nx;
1203:       from10->first = first;
1204:       from10->step  = step;
1205:       to10->n       = nx;
1206: #if defined(PETSC_USE_DEBUG)
1207:       VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1208: #endif
1209:       PetscMemcpy(to10->vslots,idy,nx*sizeof(PetscInt));
1210:       ctx->todata   = (void*)to10;
1211:       ctx->fromdata = (void*)from10;
1212:       if (step == 1) ctx->begin = VecScatterBegin_SSToSG_Stride1;
1213:       else           ctx->begin = VecScatterBegin_SSToSG;
1214:       ctx->destroy = VecScatterDestroy_SSToSG;
1215:       ctx->end     = 0;
1216:       ctx->copy    = 0;
1217:       ctx->view    = VecScatterView_SSToSG;
1218:       to10->type   = VEC_SCATTER_SEQ_GENERAL;
1219:       from10->type = VEC_SCATTER_SEQ_STRIDE;
1220:       PetscInfo(xin,"Special case: sequential vector stride to general\n");
1221:       goto functionend;
1222:     } else {
1223:       PetscInt               nx,ny;
1224:       const PetscInt         *idx,*idy;
1225:       VecScatter_Seq_General *to11 = NULL,*from11 = NULL;
1226:       PetscBool              idnx,idny;

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

1232:       ISIdentity(ix,&idnx);
1233:       ISIdentity(iy,&idny);
1234:       if (idnx && idny) {
1235:         VecScatter_Seq_Stride *to13 = NULL,*from13 = NULL;
1236:         PetscMalloc2(1,&to13,1,&from13);
1237:         to13->n       = nx;
1238:         to13->first   = 0;
1239:         to13->step    = 1;
1240:         from13->n     = nx;
1241:         from13->first = 0;
1242:         from13->step  = 1;
1243:         to13->type    = VEC_SCATTER_SEQ_STRIDE;
1244:         from13->type  = VEC_SCATTER_SEQ_STRIDE;
1245:         ctx->todata   = (void*)to13;
1246:         ctx->fromdata = (void*)from13;
1247:         ctx->begin    = VecScatterBegin_SSToSS;
1248:         ctx->end      = 0;
1249:         ctx->destroy  = VecScatterDestroy_SSToSS;
1250:         ctx->copy     = VecScatterCopy_SSToSS;
1251:         ctx->view     = VecScatterView_SSToSS;
1252:         PetscInfo(xin,"Special case: sequential copy\n");
1253:         goto functionend;
1254:       }

1256:       ISGetIndices(iy,&idy);
1257:       ISGetIndices(ix,&idx);
1258:       PetscMalloc2(1,&to11,1,&from11);
1259:       PetscMalloc2(nx,&to11->vslots,nx,&from11->vslots);
1260:       to11->n = nx;
1261: #if defined(PETSC_USE_DEBUG)
1262:       VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1263: #endif
1264:       PetscMemcpy(to11->vslots,idy,nx*sizeof(PetscInt));
1265:       from11->n = nx;
1266: #if defined(PETSC_USE_DEBUG)
1267:       VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1268: #endif
1269:       PetscMemcpy(from11->vslots,idx,nx*sizeof(PetscInt));
1270:       to11->type    = VEC_SCATTER_SEQ_GENERAL;
1271:       from11->type  = VEC_SCATTER_SEQ_GENERAL;
1272:       ctx->todata   = (void*)to11;
1273:       ctx->fromdata = (void*)from11;
1274:       ctx->begin    = VecScatterBegin_SGToSG;
1275:       ctx->end      = 0;
1276:       ctx->destroy  = VecScatterDestroy_SGToSG;
1277:       ctx->copy     = VecScatterCopy_SGToSG;
1278:       ctx->view     = VecScatterView_SGToSG;
1279:       ISRestoreIndices(ix,&idx);
1280:       ISRestoreIndices(iy,&idy);
1281:       PetscInfo(xin,"Sequential vector scatter with block indices\n");
1282:       goto functionend;
1283:     }
1284:   }
1285:   /* ---------------------------------------------------------------------------*/
1286:   if (xin_type == VEC_MPI_ID && yin_type == VEC_SEQ_ID) {

1288:     /* ===========================================================================================================
1289:           Check for special cases
1290:        ==========================================================================================================*/
1291:     islocal = PETSC_FALSE;
1292:     /* special case extracting (subset of) local portion */
1293:     if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1294:       PetscInt              nx,ny,to_first,to_step,from_first,from_step;
1295:       PetscInt              start,end,min,max;
1296:       VecScatter_Seq_Stride *from12 = NULL,*to12 = NULL;

1298:       VecGetOwnershipRange(xin,&start,&end);
1299:       ISGetLocalSize(ix,&nx);
1300:       ISStrideGetInfo(ix,&from_first,&from_step);
1301:       ISGetLocalSize(iy,&ny);
1302:       ISStrideGetInfo(iy,&to_first,&to_step);
1303:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1304:       ISGetMinMax(ix,&min,&max);
1305:       if (min >= start && max < end) islocal = PETSC_TRUE;
1306:       else islocal = PETSC_FALSE;
1307:       MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1308:       if (cando) {
1309:         PetscMalloc2(1,&to12,1,&from12);
1310:         to12->n       = nx;
1311:         to12->first   = to_first;
1312:         to12->step    = to_step;
1313:         from12->n     = nx;
1314:         from12->first = from_first-start;
1315:         from12->step  = from_step;
1316:         to12->type    = VEC_SCATTER_SEQ_STRIDE;
1317:         from12->type  = VEC_SCATTER_SEQ_STRIDE;
1318:         ctx->todata   = (void*)to12;
1319:         ctx->fromdata = (void*)from12;
1320:         ctx->begin    = VecScatterBegin_SSToSS;
1321:         ctx->end      = 0;
1322:         ctx->destroy  = VecScatterDestroy_SSToSS;
1323:         ctx->copy     = VecScatterCopy_SSToSS;
1324:         ctx->view     = VecScatterView_SSToSS;
1325:         PetscInfo(xin,"Special case: processors only getting local values\n");
1326:         goto functionend;
1327:       }
1328:     } else {
1329:       MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1330:     }

1332:     /* test for special case of all processors getting entire vector */
1333:     /* contains check that PetscMPIInt can handle the sizes needed */
1334:     totalv = PETSC_FALSE;
1335:     if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1336:       PetscInt             i,nx,ny,to_first,to_step,from_first,from_step,N;
1337:       PetscMPIInt          *count = NULL,*displx;
1338:       VecScatter_MPI_ToAll *sto   = NULL;

1340:       ISGetLocalSize(ix,&nx);
1341:       ISStrideGetInfo(ix,&from_first,&from_step);
1342:       ISGetLocalSize(iy,&ny);
1343:       ISStrideGetInfo(iy,&to_first,&to_step);
1344:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1345:       VecGetSize(xin,&N);
1346:       if (nx != N) totalv = PETSC_FALSE;
1347:       else if (from_first == 0 && from_step == 1 && from_first == to_first && from_step == to_step) totalv = PETSC_TRUE;
1348:       else totalv = PETSC_FALSE;
1349:       MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));

1351: #if defined(PETSC_USE_64BIT_INDICES)
1352:       if (cando && (yin->map->N < PETSC_MPI_INT_MAX)) {
1353: #else
1354:       if (cando) {
1355: #endif
1356:         MPI_Comm_size(PetscObjectComm((PetscObject)ctx),&size);
1357:         PetscMalloc3(1,&sto,size,&count,size,&displx);
1358:         range = xin->map->range;
1359:         for (i=0; i<size; i++) {
1360:           PetscMPIIntCast(range[i+1] - range[i],count+i);
1361:           PetscMPIIntCast(range[i],displx+i);
1362:         }
1363:         sto->count    = count;
1364:         sto->displx   = displx;
1365:         sto->work1    = 0;
1366:         sto->work2    = 0;
1367:         sto->type     = VEC_SCATTER_MPI_TOALL;
1368:         ctx->todata   = (void*)sto;
1369:         ctx->fromdata = 0;
1370:         ctx->begin    = VecScatterBegin_MPI_ToAll;
1371:         ctx->end      = 0;
1372:         ctx->destroy  = VecScatterDestroy_MPI_ToAll;
1373:         ctx->copy     = VecScatterCopy_MPI_ToAll;
1374:         ctx->view     = VecScatterView_MPI_ToAll;
1375:         PetscInfo(xin,"Special case: all processors get entire parallel vector\n");
1376:         goto functionend;
1377:       }
1378:     } else {
1379:       MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1380:     }

1382:     /* test for special case of processor 0 getting entire vector */
1383:     /* contains check that PetscMPIInt can handle the sizes needed */
1384:     totalv = PETSC_FALSE;
1385:     if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1386:       PetscInt             i,nx,ny,to_first,to_step,from_first,from_step,N;
1387:       PetscMPIInt          rank,*count = NULL,*displx;
1388:       VecScatter_MPI_ToAll *sto = NULL;

1390:       PetscObjectGetComm((PetscObject)xin,&comm);
1391:       MPI_Comm_rank(comm,&rank);
1392:       ISGetLocalSize(ix,&nx);
1393:       ISStrideGetInfo(ix,&from_first,&from_step);
1394:       ISGetLocalSize(iy,&ny);
1395:       ISStrideGetInfo(iy,&to_first,&to_step);
1396:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1397:       if (!rank) {
1398:         VecGetSize(xin,&N);
1399:         if (nx != N) totalv = PETSC_FALSE;
1400:         else if (from_first == 0        && from_step == 1 &&
1401:                  from_first == to_first && from_step == to_step) totalv = PETSC_TRUE;
1402:         else totalv = PETSC_FALSE;
1403:       } else {
1404:         if (!nx) totalv = PETSC_TRUE;
1405:         else     totalv = PETSC_FALSE;
1406:       }
1407:       MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));

1409: #if defined(PETSC_USE_64BIT_INDICES)
1410:       if (cando && (yin->map->N < PETSC_MPI_INT_MAX)) {
1411: #else
1412:       if (cando) {
1413: #endif
1414:         MPI_Comm_size(PetscObjectComm((PetscObject)ctx),&size);
1415:         PetscMalloc3(1,&sto,size,&count,size,&displx);
1416:         range = xin->map->range;
1417:         for (i=0; i<size; i++) {
1418:           PetscMPIIntCast(range[i+1] - range[i],count+i);
1419:           PetscMPIIntCast(range[i],displx+i);
1420:         }
1421:         sto->count    = count;
1422:         sto->displx   = displx;
1423:         sto->work1    = 0;
1424:         sto->work2    = 0;
1425:         sto->type     = VEC_SCATTER_MPI_TOONE;
1426:         ctx->todata   = (void*)sto;
1427:         ctx->fromdata = 0;
1428:         ctx->begin    = VecScatterBegin_MPI_ToOne;
1429:         ctx->end      = 0;
1430:         ctx->destroy  = VecScatterDestroy_MPI_ToAll;
1431:         ctx->copy     = VecScatterCopy_MPI_ToAll;
1432:         ctx->view     = VecScatterView_MPI_ToAll;
1433:         PetscInfo(xin,"Special case: processor zero gets entire parallel vector, rest get none\n");
1434:         goto functionend;
1435:       }
1436:     } else {
1437:       MPI_Allreduce(&totalv,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)xin));
1438:     }

1440:     PetscObjectTypeCompare((PetscObject)ix,ISBLOCK,&ixblock);
1441:     PetscObjectTypeCompare((PetscObject)iy,ISBLOCK,&iyblock);
1442:     PetscObjectTypeCompare((PetscObject)iy,ISSTRIDE,&iystride);
1443:     if (ixblock) {
1444:       /* special case block to block */
1445:       if (iyblock) {
1446:         PetscInt       nx,ny,bsx,bsy;
1447:         const PetscInt *idx,*idy;
1448:         ISGetBlockSize(iy,&bsy);
1449:         ISGetBlockSize(ix,&bsx);
1450:         if (bsx == bsy && VecScatterOptimizedBS(bsx)) {
1451:           ISBlockGetLocalSize(ix,&nx);
1452:           ISBlockGetIndices(ix,&idx);
1453:           ISBlockGetLocalSize(iy,&ny);
1454:           ISBlockGetIndices(iy,&idy);
1455:           if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1456:           VecScatterCreate_PtoS(nx,idx,ny,idy,xin,yin,bsx,ctx);
1457:           ISBlockRestoreIndices(ix,&idx);
1458:           ISBlockRestoreIndices(iy,&idy);
1459:           PetscInfo(xin,"Special case: blocked indices\n");
1460:           goto functionend;
1461:         }
1462:         /* special case block to stride */
1463:       } else if (iystride) {
1464:         PetscInt ystart,ystride,ysize,bsx;
1465:         ISStrideGetInfo(iy,&ystart,&ystride);
1466:         ISGetLocalSize(iy,&ysize);
1467:         ISGetBlockSize(ix,&bsx);
1468:         /* see if stride index set is equivalent to block index set */
1469:         if (VecScatterOptimizedBS(bsx) && ((ystart % bsx) == 0) && (ystride == 1) && ((ysize % bsx) == 0)) {
1470:           PetscInt       nx,il,*idy;
1471:           const PetscInt *idx;
1472:           ISBlockGetLocalSize(ix,&nx);
1473:           ISBlockGetIndices(ix,&idx);
1474:           if (ysize != bsx*nx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1475:           PetscMalloc1(nx,&idy);
1476:           if (nx) {
1477:             idy[0] = ystart/bsx;
1478:             for (il=1; il<nx; il++) idy[il] = idy[il-1] + 1;
1479:           }
1480:           VecScatterCreate_PtoS(nx,idx,nx,idy,xin,yin,bsx,ctx);
1481:           PetscFree(idy);
1482:           ISBlockRestoreIndices(ix,&idx);
1483:           PetscInfo(xin,"Special case: blocked indices to stride\n");
1484:           goto functionend;
1485:         }
1486:       }
1487:     }
1488:     /* left over general case */
1489:     {
1490:       PetscInt       nx,ny;
1491:       const PetscInt *idx,*idy;
1492:       ISGetLocalSize(ix,&nx);
1493:       ISGetIndices(ix,&idx);
1494:       ISGetLocalSize(iy,&ny);
1495:       ISGetIndices(iy,&idy);
1496:       if (nx != ny) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match (%d %d)",nx,ny);
1497:       VecScatterCreate_PtoS(nx,idx,ny,idy,xin,yin,1,ctx);
1498:       ISRestoreIndices(ix,&idx);
1499:       ISRestoreIndices(iy,&idy);
1500:       PetscInfo(xin,"General case: MPI to Seq\n");
1501:       goto functionend;
1502:     }
1503:   }
1504:   /* ---------------------------------------------------------------------------*/
1505:   if (xin_type == VEC_SEQ_ID && yin_type == VEC_MPI_ID) {
1506:     /* ===========================================================================================================
1507:           Check for special cases
1508:        ==========================================================================================================*/
1509:     /* special case local copy portion */
1510:     islocal = PETSC_FALSE;
1511:     if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID) {
1512:       PetscInt              nx,ny,to_first,to_step,from_step,start,end,from_first,min,max;
1513:       VecScatter_Seq_Stride *from = NULL,*to = NULL;

1515:       VecGetOwnershipRange(yin,&start,&end);
1516:       ISGetLocalSize(ix,&nx);
1517:       ISStrideGetInfo(ix,&from_first,&from_step);
1518:       ISGetLocalSize(iy,&ny);
1519:       ISStrideGetInfo(iy,&to_first,&to_step);
1520:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1521:       ISGetMinMax(iy,&min,&max);
1522:       if (min >= start && max < end) islocal = PETSC_TRUE;
1523:       else islocal = PETSC_FALSE;
1524:       MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)yin));
1525:       if (cando) {
1526:         PetscMalloc2(1,&to,1,&from);
1527:         to->n         = nx;
1528:         to->first     = to_first-start;
1529:         to->step      = to_step;
1530:         from->n       = nx;
1531:         from->first   = from_first;
1532:         from->step    = from_step;
1533:         to->type      = VEC_SCATTER_SEQ_STRIDE;
1534:         from->type    = VEC_SCATTER_SEQ_STRIDE;
1535:         ctx->todata   = (void*)to;
1536:         ctx->fromdata = (void*)from;
1537:         ctx->begin    = VecScatterBegin_SSToSS;
1538:         ctx->end      = 0;
1539:         ctx->destroy  = VecScatterDestroy_SSToSS;
1540:         ctx->copy     = VecScatterCopy_SSToSS;
1541:         ctx->view     = VecScatterView_SSToSS;
1542:         PetscInfo(xin,"Special case: sequential stride to MPI stride\n");
1543:         goto functionend;
1544:       }
1545:     } else {
1546:       MPI_Allreduce(&islocal,&cando,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)yin));
1547:     }
1548:     /* special case block to stride */
1549:     if (ix_type == IS_BLOCK_ID && iy_type == IS_STRIDE_ID) {
1550:       PetscInt ystart,ystride,ysize,bsx;
1551:       ISStrideGetInfo(iy,&ystart,&ystride);
1552:       ISGetLocalSize(iy,&ysize);
1553:       ISGetBlockSize(ix,&bsx);
1554:       /* see if stride index set is equivalent to block index set */
1555:       if (VecScatterOptimizedBS(bsx) && ((ystart % bsx) == 0) && (ystride == 1) && ((ysize % bsx) == 0)) {
1556:         PetscInt       nx,il,*idy;
1557:         const PetscInt *idx;
1558:         ISBlockGetLocalSize(ix,&nx);
1559:         ISBlockGetIndices(ix,&idx);
1560:         if (ysize != bsx*nx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1561:         PetscMalloc1(nx,&idy);
1562:         if (nx) {
1563:           idy[0] = ystart/bsx;
1564:           for (il=1; il<nx; il++) idy[il] = idy[il-1] + 1;
1565:         }
1566:         VecScatterCreate_StoP(nx,idx,nx,idy,xin,yin,bsx,ctx);
1567:         PetscFree(idy);
1568:         ISBlockRestoreIndices(ix,&idx);
1569:         PetscInfo(xin,"Special case: Blocked indices to stride\n");
1570:         goto functionend;
1571:       }
1572:     }

1574:     /* general case */
1575:     {
1576:       PetscInt       nx,ny;
1577:       const PetscInt *idx,*idy;
1578:       ISGetLocalSize(ix,&nx);
1579:       ISGetIndices(ix,&idx);
1580:       ISGetLocalSize(iy,&ny);
1581:       ISGetIndices(iy,&idy);
1582:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1583:       VecScatterCreate_StoP(nx,idx,ny,idy,xin,yin,1,ctx);
1584:       ISRestoreIndices(ix,&idx);
1585:       ISRestoreIndices(iy,&idy);
1586:       PetscInfo(xin,"General case: Seq to MPI\n");
1587:       goto functionend;
1588:     }
1589:   }
1590:   /* ---------------------------------------------------------------------------*/
1591:   if (xin_type == VEC_MPI_ID && yin_type == VEC_MPI_ID) {
1592:     /* no special cases for now */
1593:     PetscInt       nx,ny;
1594:     const PetscInt *idx,*idy;
1595:     ISGetLocalSize(ix,&nx);
1596:     ISGetIndices(ix,&idx);
1597:     ISGetLocalSize(iy,&ny);
1598:     ISGetIndices(iy,&idy);
1599:     if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1600:     VecScatterCreate_PtoP(nx,idx,ny,idy,xin,yin,1,ctx);
1601:     ISRestoreIndices(ix,&idx);
1602:     ISRestoreIndices(iy,&idy);
1603:     PetscInfo(xin,"General case: MPI to MPI\n");
1604:     goto functionend;
1605:   }

1607:   functionend:
1608:   *newctx = ctx;
1609:   ISDestroy(&tix);
1610:   ISDestroy(&tiy);
1611:   VecScatterViewFromOptions(ctx,NULL,"-vecscatter_view");
1612:   return(0);
1613: }

1615: /* ------------------------------------------------------------------*/
1618: /*@
1619:    VecScatterGetMerged - Returns true if the scatter is completed in the VecScatterBegin()
1620:       and the VecScatterEnd() does nothing

1622:    Not Collective

1624:    Input Parameter:
1625: .   ctx - scatter context created with VecScatterCreate()

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

1630:    Level: developer

1632: .seealso: VecScatterCreate(), VecScatterEnd(), VecScatterBegin()
1633: @*/
1634: PetscErrorCode  VecScatterGetMerged(VecScatter ctx,PetscBool  *flg)
1635: {
1638:   *flg = ctx->beginandendtogether;
1639:   return(0);
1640: }

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

1648:    Neighbor-wise Collective on VecScatter and Vec

1650:    Input Parameters:
1651: +  inctx - scatter context generated by VecScatterCreate()
1652: .  x - the vector from which we scatter
1653: .  y - the vector to which we scatter
1654: .  addv - either ADD_VALUES or INSERT_VALUES, with INSERT_VALUES mode any location
1655:           not scattered to retains its old value; i.e. the vector is NOT first zeroed.
1656: -  mode - the scattering mode, usually SCATTER_FORWARD.  The available modes are:
1657:     SCATTER_FORWARD or SCATTER_REVERSE


1660:    Level: intermediate

1662:    Options Database: See VecScatterCreate()

1664:    Notes:
1665:    The vectors x and y need not be the same vectors used in the call
1666:    to VecScatterCreate(), but x must have the same parallel data layout
1667:    as that passed in as the x to VecScatterCreate(), similarly for the y.
1668:    Most likely they have been obtained from VecDuplicate().

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

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

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

1678:    This scatter is far more general than the conventional
1679:    scatter, since it can be a gather or a scatter or a combination,
1680:    depending on the indices ix and iy.  If x is a parallel vector and y
1681:    is sequential, VecScatterBegin() can serve to gather values to a
1682:    single processor.  Similarly, if y is parallel and x sequential, the
1683:    routine can scatter from one processor to many processors.

1685:    Concepts: scatter^between vectors
1686:    Concepts: gather^between vectors

1688: .seealso: VecScatterCreate(), VecScatterEnd()
1689: @*/
1690: PetscErrorCode  VecScatterBegin(VecScatter inctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
1691: {
1693: #if defined(PETSC_USE_DEBUG)
1694:   PetscInt       to_n,from_n;
1695: #endif
1700:   if (inctx->inuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE," Scatter ctx already in use");

1702: #if defined(PETSC_USE_DEBUG)
1703:   /*
1704:      Error checking to make sure these vectors match the vectors used
1705:    to create the vector scatter context. -1 in the from_n and to_n indicate the
1706:    vector lengths are unknown (for example with mapped scatters) and thus
1707:    no error checking is performed.
1708:   */
1709:   if (inctx->from_n >= 0 && inctx->to_n >= 0) {
1710:     VecGetLocalSize(x,&from_n);
1711:     VecGetLocalSize(y,&to_n);
1712:     if (mode & SCATTER_REVERSE) {
1713:       if (to_n != inctx->from_n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter reverse and vector to != ctx from size)",to_n,inctx->from_n);
1714:       if (from_n != inctx->to_n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter reverse and vector from != ctx to size)",from_n,inctx->to_n);
1715:     } else {
1716:       if (to_n != inctx->to_n)     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter forward and vector to != ctx to size)",to_n,inctx->to_n);
1717:       if (from_n != inctx->from_n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter forward and vector from != ctx from size)",from_n,inctx->from_n);
1718:     }
1719:   }
1720: #endif

1722:   inctx->inuse = PETSC_TRUE;
1723:   PetscLogEventBarrierBegin(VEC_ScatterBarrier,0,0,0,0,PetscObjectComm((PetscObject)inctx));
1724:   (*inctx->begin)(inctx,x,y,addv,mode);
1725:   if (inctx->beginandendtogether && inctx->end) {
1726:     inctx->inuse = PETSC_FALSE;
1727:     (*inctx->end)(inctx,x,y,addv,mode);
1728:   }
1729:   PetscLogEventBarrierEnd(VEC_ScatterBarrier,0,0,0,0,PetscObjectComm((PetscObject)inctx));
1730:   return(0);
1731: }

1733: /* --------------------------------------------------------------------*/
1736: /*@
1737:    VecScatterEnd - Ends a generalized scatter from one vector to another.  Call
1738:    after first calling VecScatterBegin().

1740:    Neighbor-wise Collective on VecScatter and Vec

1742:    Input Parameters:
1743: +  ctx - scatter context generated by VecScatterCreate()
1744: .  x - the vector from which we scatter
1745: .  y - the vector to which we scatter
1746: .  addv - either ADD_VALUES or INSERT_VALUES.
1747: -  mode - the scattering mode, usually SCATTER_FORWARD.  The available modes are:
1748:      SCATTER_FORWARD, SCATTER_REVERSE

1750:    Level: intermediate

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

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

1757: .seealso: VecScatterBegin(), VecScatterCreate()
1758: @*/
1759: PetscErrorCode  VecScatterEnd(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
1760: {

1767:   ctx->inuse = PETSC_FALSE;
1768:   if (!ctx->end) return(0);
1769:   if (!ctx->beginandendtogether) {
1770:     PetscLogEventBegin(VEC_ScatterEnd,ctx,x,y,0);
1771:     (*(ctx)->end)(ctx,x,y,addv,mode);
1772:     PetscLogEventEnd(VEC_ScatterEnd,ctx,x,y,0);
1773:   }
1774:   return(0);
1775: }

1779: /*@C
1780:    VecScatterDestroy - Destroys a scatter context created by
1781:    VecScatterCreate().

1783:    Collective on VecScatter

1785:    Input Parameter:
1786: .  ctx - the scatter context

1788:    Level: intermediate

1790: .seealso: VecScatterCreate(), VecScatterCopy()
1791: @*/
1792: PetscErrorCode  VecScatterDestroy(VecScatter *ctx)
1793: {

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

1802:   /* if memory was published with SAWs then destroy it */
1803:   PetscObjectSAWsViewOff((PetscObject)(*ctx));
1804:   (*(*ctx)->destroy)(*ctx);
1805: #if defined(PETSC_HAVE_CUSP)
1806:   VecScatterCUSPIndicesDestroy((PetscCUSPIndices*)&((*ctx)->spptr));
1807: #endif
1808:   PetscHeaderDestroy(ctx);
1809:   return(0);
1810: }

1814: /*@
1815:    VecScatterCopy - Makes a copy of a scatter context.

1817:    Collective on VecScatter

1819:    Input Parameter:
1820: .  sctx - the scatter context

1822:    Output Parameter:
1823: .  ctx - the context copy

1825:    Level: advanced

1827: .seealso: VecScatterCreate(), VecScatterDestroy()
1828: @*/
1829: PetscErrorCode  VecScatterCopy(VecScatter sctx,VecScatter *ctx)
1830: {

1836:   if (!sctx->copy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot copy this type");
1837:   PetscHeaderCreate(*ctx,_p_VecScatter,int,VEC_SCATTER_CLASSID,"VecScatter","VecScatter","Vec",PetscObjectComm((PetscObject)sctx),VecScatterDestroy,VecScatterView);
1838:   (*ctx)->to_n   = sctx->to_n;
1839:   (*ctx)->from_n = sctx->from_n;
1840:   (*sctx->copy)(sctx,*ctx);
1841:   return(0);
1842: }


1845: /* ------------------------------------------------------------------*/
1848: /*@
1849:    VecScatterView - Views a vector scatter context.

1851:    Collective on VecScatter

1853:    Input Parameters:
1854: +  ctx - the scatter context
1855: -  viewer - the viewer for displaying the context

1857:    Level: intermediate

1859: @*/
1860: PetscErrorCode  VecScatterView(VecScatter ctx,PetscViewer viewer)
1861: {

1866:   if (!viewer) {
1867:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ctx),&viewer);
1868:   }
1870:   if (ctx->view) {
1871:     (*ctx->view)(ctx,viewer);
1872:   }
1873:   return(0);
1874: }

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

1882:    Collective on VecScatter

1884:    Input Parameters:
1885: +  scat - vector scatter context
1886: .  from - remapping for "from" indices (may be NULL)
1887: -  to   - remapping for "to" indices (may be NULL)

1889:    Level: developer

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

1895:           In the sequential case the todata is the indices where the
1896:           data is put and the fromdata is where it is taken from.
1897:           This is backwards from the paralllel case! CRY! CRY! CRY!

1899: @*/
1900: PetscErrorCode  VecScatterRemap(VecScatter scat,PetscInt *rto,PetscInt *rfrom)
1901: {
1902:   VecScatter_Seq_General *to,*from;
1903:   VecScatter_MPI_General *mto;
1904:   PetscInt               i;


1911:   from = (VecScatter_Seq_General*)scat->fromdata;
1912:   mto  = (VecScatter_MPI_General*)scat->todata;

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

1916:   if (rto) {
1917:     if (mto->type == VEC_SCATTER_MPI_GENERAL) {
1918:       /* handle off processor parts */
1919:       for (i=0; i<mto->starts[mto->n]; i++) mto->indices[i] = rto[mto->indices[i]];

1921:       /* handle local part */
1922:       to = &mto->local;
1923:       for (i=0; i<to->n; i++) to->vslots[i] = rto[to->vslots[i]];
1924:     } else if (from->type == VEC_SCATTER_SEQ_GENERAL) {
1925:       for (i=0; i<from->n; i++) from->vslots[i] = rto[from->vslots[i]];
1926:     } else if (from->type == VEC_SCATTER_SEQ_STRIDE) {
1927:       VecScatter_Seq_Stride *sto = (VecScatter_Seq_Stride*)from;

1929:       /* if the remapping is the identity and stride is identity then skip remap */
1930:       if (sto->step == 1 && sto->first == 0) {
1931:         for (i=0; i<sto->n; i++) {
1932:           if (rto[i] != i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1933:         }
1934:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1935:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1936:   }

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

1940:   /*
1941:      Mark then vector lengths as unknown because we do not know the
1942:    lengths of the remapped vectors
1943:   */
1944:   scat->from_n = -1;
1945:   scat->to_n   = -1;
1946:   return(0);
1947: }