Actual source code: vscat.c

petsc-3.3-p7 2013-05-11
  2: /*
  3:      Code for creating scatters between vectors. This file 
  4:   includes the code for scattering between sequential vectors and
  5:   some special cases for parallel scatters.
  6: */

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

 11: /* Logging support */
 12: PetscClassId  VEC_SCATTER_CLASSID;

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

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

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

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

 48:   VecGetLocalSize(y,&yy_n);
 49:   VecGetLocalSize(x,&xx_n);
 50:   VecGetArray(y,&yv);
 51:   if (x != y) {VecGetArray(x,&xv);}

 53:   if (mode & SCATTER_REVERSE) {
 54:     PetscScalar          *xvt,*xvt2;
 55:     VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
 56:     PetscInt             i;
 57:     PetscMPIInt          *disply = scat->displx;

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

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

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

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

160:   VecGetLocalSize(y,&yy_n);
161:   VecGetLocalSize(x,&xx_n);
162:   VecGetArray(x,&xv);
163:   VecGetArray(y,&yv);

165:   PetscObjectGetComm((PetscObject)x,&comm);
166:   MPI_Comm_rank(comm,&rank);

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

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

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

234: /*
235:        The follow to are used for both VecScatterBegin_MPI_ToAll() and VecScatterBegin_MPI_ToOne()
236: */
239: PetscErrorCode VecScatterDestroy_MPI_ToAll(VecScatter ctx)
240: {
241:   VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
242:   PetscErrorCode       ierr;

245:   PetscFree(scat->work1);
246:   PetscFree(scat->work2);
247:   PetscFree3(ctx->todata,scat->count,scat->displx);
248:   return(0);
249: }

253: PetscErrorCode VecScatterDestroy_SGtoSG(VecScatter ctx)
254: {
255:   PetscErrorCode         ierr;

258:   PetscFree2(((VecScatter_Seq_General*)ctx->todata)->vslots,((VecScatter_Seq_General*)ctx->fromdata)->vslots);
259:   PetscFree2(ctx->todata,ctx->fromdata);
260:   return(0);
261: }

265: PetscErrorCode VecScatterDestroy_SGtoSS(VecScatter ctx)
266: {
267:   PetscErrorCode         ierr;

270:   PetscFree(((VecScatter_Seq_General*)ctx->fromdata)->vslots);
271:   PetscFree2(ctx->todata,ctx->fromdata);
272:   return(0);
273: }

277: PetscErrorCode VecScatterDestroy_SStoSG(VecScatter ctx)
278: {
279:   PetscErrorCode         ierr;

282:   PetscFree(((VecScatter_Seq_General*)ctx->todata)->vslots);
283:   PetscFree2(ctx->todata,ctx->fromdata);
284:   return(0);
285: }

289: PetscErrorCode VecScatterDestroy_SStoSS(VecScatter ctx)
290: {

294:   PetscFree2(ctx->todata,ctx->fromdata);
295:   return(0);
296: }

298: /* -------------------------------------------------------------------------*/
301: PetscErrorCode VecScatterCopy_MPI_ToAll(VecScatter in,VecScatter out)
302: {
303:   VecScatter_MPI_ToAll *in_to = (VecScatter_MPI_ToAll*)in->todata,*sto;
304:   PetscErrorCode       ierr;
305:   PetscMPIInt          size,*count,*displx;
306:   PetscInt             i;

309:   out->begin          = in->begin;
310:   out->end            = in->end;
311:   out->copy           = in->copy;
312:   out->destroy        = in->destroy;
313:   out->view           = in->view;

315:   MPI_Comm_size(((PetscObject)out)->comm,&size);
316:   PetscMalloc3(1,VecScatter_MPI_ToAll,&sto,size,PetscMPIInt,&count,size,PetscMPIInt,&displx);
317:   sto->type           = in_to->type;
318:   sto->count          = count;
319:   sto->displx         = displx;
320:   for (i=0; i<size; i++) {
321:     sto->count[i]  = in_to->count[i];
322:     sto->displx[i] = in_to->displx[i];
323:   }
324:   sto->work1         = 0;
325:   sto->work2         = 0;
326:   out->todata        = (void*)sto;
327:   out->fromdata      = (void*)0;
328:   return(0);
329: }

331: /* --------------------------------------------------------------------------------------*/
332: /* 
333:         Scatter: sequential general to sequential general 
334: */
337: PetscErrorCode VecScatterBegin_SGtoSG(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
338: {
339:   VecScatter_Seq_General *gen_to = (VecScatter_Seq_General*)ctx->todata;
340:   VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
341:   PetscErrorCode         ierr;
342:   PetscInt               i,n = gen_from->n,*fslots,*tslots;
343:   PetscScalar            *xv,*yv;
344: 
346:   VecGetArray(x,&xv);
347:   if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}
348:   if (mode & SCATTER_REVERSE){
349:     gen_to   = (VecScatter_Seq_General*)ctx->fromdata;
350:     gen_from = (VecScatter_Seq_General*)ctx->todata;
351:   }
352:   fslots = gen_from->vslots;
353:   tslots = gen_to->vslots;

355:   if (addv == INSERT_VALUES) {
356:     for (i=0; i<n; i++) {yv[tslots[i]] = xv[fslots[i]];}
357:   } else if (addv == ADD_VALUES) {
358:     for (i=0; i<n; i++) {yv[tslots[i]] += xv[fslots[i]];}
359: #if !defined(PETSC_USE_COMPLEX)
360:   } else  if (addv == MAX_VALUES) {
361:     for (i=0; i<n; i++) {yv[tslots[i]] = PetscMax(yv[tslots[i]],xv[fslots[i]]);}
362: #endif
363:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");
364:   VecRestoreArray(x,&xv);
365:   if (x != y) {VecRestoreArray(y,&yv);}
366:   return(0);
367: }

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

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

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

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

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

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

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

548: /* 
549:      Scatter: sequential stride to sequential stride 
550: */
553: PetscErrorCode VecScatterBegin_SStoSS(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
554: {
555:   VecScatter_Seq_Stride *gen_to   = (VecScatter_Seq_Stride*)ctx->todata;
556:   VecScatter_Seq_Stride *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
557:   PetscInt              i,n = gen_from->n,to_first = gen_to->first,to_step = gen_to->step;
558:   PetscErrorCode        ierr;
559:   PetscInt              from_first = gen_from->first,from_step = gen_from->step;
560:   PetscScalar           *xv,*yv;
561: 
563:   VecGetArrayRead(x,(const PetscScalar **)&xv);
564:   if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}

566:   if (mode & SCATTER_REVERSE){
567:     from_first = gen_to->first;
568:     to_first   = gen_from->first;
569:     from_step  = gen_to->step;
570:     to_step    = gen_from->step;
571:   }

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

611: /* --------------------------------------------------------------------------*/


616: PetscErrorCode VecScatterCopy_SGToSG(VecScatter in,VecScatter out)
617: {
618:   PetscErrorCode         ierr;
619:   VecScatter_Seq_General *in_to   = (VecScatter_Seq_General*)in->todata,*out_to = PETSC_NULL;
620:   VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from = PETSC_NULL;
621: 
623:   out->begin         = in->begin;
624:   out->end           = in->end;
625:   out->copy          = in->copy;
626:   out->destroy       = in->destroy;
627:   out->view          = in->view;

629:   PetscMalloc2(1,VecScatter_Seq_General,&out_to,1,VecScatter_Seq_General,&out_from);
630:   PetscMalloc2(in_to->n,PetscInt,&out_to->vslots,in_from->n,PetscInt,&out_from->vslots);
631:   out_to->n                      = in_to->n;
632:   out_to->type                   = in_to->type;
633:   out_to->nonmatching_computed   = PETSC_FALSE;
634:   out_to->slots_nonmatching      = 0;
635:   out_to->is_copy                = PETSC_FALSE;
636:   PetscMemcpy(out_to->vslots,in_to->vslots,(out_to->n)*sizeof(PetscInt));

638:   out_from->n                    = in_from->n;
639:   out_from->type                 = in_from->type;
640:   out_from->nonmatching_computed = PETSC_FALSE;
641:   out_from->slots_nonmatching    = 0;
642:   out_from->is_copy              = PETSC_FALSE;
643:   PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));

645:   out->todata     = (void*)out_to;
646:   out->fromdata   = (void*)out_from;
647:   return(0);
648: }


653: PetscErrorCode VecScatterCopy_SGToSS(VecScatter in,VecScatter out)
654: {
655:   PetscErrorCode         ierr;
656:   VecScatter_Seq_Stride  *in_to   = (VecScatter_Seq_Stride*)in->todata,*out_to = PETSC_NULL;
657:   VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from = PETSC_NULL;
658: 
660:   out->begin         = in->begin;
661:   out->end           = in->end;
662:   out->copy          = in->copy;
663:   out->destroy       = in->destroy;
664:   out->view          = in->view;

666:   PetscMalloc2(1,VecScatter_Seq_Stride,&out_to,1,VecScatter_Seq_General,&out_from);
667:   PetscMalloc(in_from->n*sizeof(PetscInt),&out_from->vslots);
668:   out_to->n       = in_to->n;
669:   out_to->type    = in_to->type;
670:   out_to->first   = in_to->first;
671:   out_to->step    = in_to->step;
672:   out_to->type    = in_to->type;

674:   out_from->n                    = in_from->n;
675:   out_from->type                 = in_from->type;
676:   out_from->nonmatching_computed = PETSC_FALSE;
677:   out_from->slots_nonmatching    = 0;
678:   out_from->is_copy              = PETSC_FALSE;
679:   PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));

681:   out->todata     = (void*)out_to;
682:   out->fromdata   = (void*)out_from;
683:   return(0);
684: }

686: /* --------------------------------------------------------------------------*/
687: /* 
688:     Scatter: parallel to sequential vector, sequential strides for both. 
689: */
692: PetscErrorCode VecScatterCopy_SStoSS(VecScatter in,VecScatter out)
693: {
694:   VecScatter_Seq_Stride *in_to   = (VecScatter_Seq_Stride*)in->todata,*out_to = PETSC_NULL;
695:   VecScatter_Seq_Stride *in_from = (VecScatter_Seq_Stride*)in->fromdata,*out_from = PETSC_NULL;
696:   PetscErrorCode        ierr;

699:   out->begin      = in->begin;
700:   out->end        = in->end;
701:   out->copy       = in->copy;
702:   out->destroy    = in->destroy;
703:   out->view       = in->view;

705:   PetscMalloc2(1,VecScatter_Seq_Stride,&out_to,1,VecScatter_Seq_Stride,&out_from);
706:   out_to->n       = in_to->n;
707:   out_to->type    = in_to->type;
708:   out_to->first   = in_to->first;
709:   out_to->step    = in_to->step;
710:   out_to->type    = in_to->type;
711:   out_from->n     = in_from->n;
712:   out_from->type  = in_from->type;
713:   out_from->first = in_from->first;
714:   out_from->step  = in_from->step;
715:   out_from->type  = in_from->type;
716:   out->todata     = (void*)out_to;
717:   out->fromdata   = (void*)out_from;
718:   return(0);
719: }

721: extern PetscErrorCode VecScatterCreate_PtoS(PetscInt,const PetscInt *,PetscInt,const PetscInt *,Vec,Vec,PetscInt,VecScatter);
722: extern PetscErrorCode VecScatterCreate_PtoP(PetscInt,const PetscInt *,PetscInt,const PetscInt *,Vec,Vec,PetscInt,VecScatter);
723: extern PetscErrorCode VecScatterCreate_StoP(PetscInt,const PetscInt *,PetscInt,const PetscInt *,Vec,Vec,PetscInt,VecScatter);

725: /* =======================================================================*/
726: #define VEC_SEQ_ID 0
727: #define VEC_MPI_ID 1
728: #define IS_GENERAL_ID 0
729: #define IS_STRIDE_ID  1
730: #define IS_BLOCK_ID   2

732: /*
733:    Blocksizes we have optimized scatters for 
734: */

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

738: PetscErrorCode  VecScatterCreateEmpty(MPI_Comm comm,VecScatter *newctx)
739: {
740:   VecScatter     ctx;

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

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

765:    Collective on Vec

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

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

778:    Options Database Keys: (uses regular MPI_Sends by default)
779: +  -vecscatter_view         - Prints detail of communications
780: .  -vecscatter_view_info    - Print less details about communication
781: .  -vecscatter_ssend        - Uses MPI_Ssend_init() instead of MPI_Send_init() 
782: .  -vecscatter_rsend           - use ready receiver mode for MPI sends 
783: .  -vecscatter_merge        - VecScatterBegin() handles all of the communication, VecScatterEnd() is a nop 
784:                               eliminates the chance for overlap of computation and communication 
785: .  -vecscatter_sendfirst    - Posts sends before receives 
786: .  -vecscatter_packtogether - Pack all messages before sending, receive all messages before unpacking
787: .  -vecscatter_alltoall     - Uses MPI all to all communication for scatter
788: .  -vecscatter_window       - Use MPI 2 window operations to move data
789: .  -vecscatter_nopack       - Avoid packing to work vector when possible (if used with -vecscatter_alltoall then will use MPI_Alltoallw()
790: -  -vecscatter_reproduce    - insure that the order of the communications are done the same for each scatter, this under certain circumstances
791:                               will make the results of scatters deterministic when otherwise they are not (it may be slower also).

793: $
794: $                                                                                    --When packing is used--
795: $                               MPI Datatypes (no packing)  sendfirst   merge        packtogether  persistent*    
796: $                                _nopack                   _sendfirst    _merge      _packtogether                -vecscatter_
797: $ ----------------------------------------------------------------------------------------------------------------------------
798: $    Message passing    Send       p                           X            X           X         always
799: $                      Ssend       p                           X            X           X         always          _ssend
800: $                      Rsend       p                        nonsense        X           X         always          _rsend
801: $    AlltoAll  v or w              X                        nonsense     always         X         nonsense        _alltoall
802: $    MPI_Win                       p                        nonsense        p           p         nonsense        _window
803: $                              
804: $   Since persistent sends and receives require a constant memory address they can only be used when data is packed into the work vector
805: $   because the in and out array may be different for each call to VecScatterBegin/End().
806: $
807: $    p indicates possible, but not implemented. X indicates implemented
808: $

810:     Level: intermediate

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

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

826:    Both ix and iy cannot be PETSC_NULL at the same time.

828:    Concepts: scatter^between vectors
829:    Concepts: gather^between vectors

831: .seealso: VecScatterDestroy(), VecScatterCreateToAll(), VecScatterCreateToZero()
832: @*/
833: PetscErrorCode  VecScatterCreate(Vec xin,IS ix,Vec yin,IS iy,VecScatter *newctx)
834: {
835:   VecScatter     ctx;
837:   PetscMPIInt    size;
838:   PetscInt       totalv,xin_type = VEC_SEQ_ID,yin_type = VEC_SEQ_ID,*range;
839:   PetscInt       ix_type = IS_GENERAL_ID,iy_type = IS_GENERAL_ID;
840:   MPI_Comm       comm,ycomm;
841:   PetscBool      ixblock,iyblock,iystride,islocal,cando,flag;
842:   IS             tix = 0,tiy = 0;

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

847:   /*
848:       Determine if the vectors are "parallel", ie. it shares a comm with other processors, or
849:       sequential (it does not share a comm). The difference is that parallel vectors treat the 
850:       index set as providing indices in the global parallel numbering of the vector, with 
851:       sequential vectors treat the index set as providing indices in the local sequential
852:       numbering
853:   */
854:   PetscObjectGetComm((PetscObject)xin,&comm);
855:   MPI_Comm_size(comm,&size);
856:   if (size > 1) {xin_type = VEC_MPI_ID;}

858:   PetscObjectGetComm((PetscObject)yin,&ycomm);
859:   MPI_Comm_size(ycomm,&size);
860:   if (size > 1) {comm = ycomm; yin_type = VEC_MPI_ID;}

862:   /* generate the Scatter context */
863:   PetscHeaderCreate(ctx,_p_VecScatter,int,VEC_SCATTER_CLASSID,0,"VecScatter","VecScatter","Vec",comm,VecScatterDestroy,VecScatterView);
864:   ctx->inuse               = PETSC_FALSE;

866:   ctx->beginandendtogether = PETSC_FALSE;
867:   PetscOptionsGetBool(PETSC_NULL,"-vecscatter_merge",&ctx->beginandendtogether,PETSC_NULL);
868:   if (ctx->beginandendtogether) {
869:     PetscInfo(ctx,"Using combined (merged) vector scatter begin and end\n");
870:   }
871:   ctx->packtogether = PETSC_FALSE;
872:   PetscOptionsGetBool(PETSC_NULL,"-vecscatter_packtogether",&ctx->packtogether,PETSC_NULL);
873:   if (ctx->packtogether) {
874:     PetscInfo(ctx,"Pack all messages before sending\n");
875:   }

877:   VecGetLocalSize(xin,&ctx->from_n);
878:   VecGetLocalSize(yin,&ctx->to_n);

880:   /*
881:       if ix or iy is not included; assume just grabbing entire vector
882:   */
883:   if (!ix && xin_type == VEC_SEQ_ID) {
884:     ISCreateStride(comm,ctx->from_n,0,1,&ix);
885:     tix  = ix;
886:   } else if (!ix && xin_type == VEC_MPI_ID) {
887:     if (yin_type == VEC_MPI_ID) {
888:       PetscInt ntmp, low;
889:       VecGetLocalSize(xin,&ntmp);
890:       VecGetOwnershipRange(xin,&low,PETSC_NULL);
891:       ISCreateStride(comm,ntmp,low,1,&ix);
892:     } else{
893:       PetscInt Ntmp;
894:       VecGetSize(xin,&Ntmp);
895:       ISCreateStride(comm,Ntmp,0,1,&ix);
896:     }
897:     tix  = ix;
898:   } else if (!ix) {
899:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ix not given, but not Seq or MPI vector");
900:   }

902:   if (!iy && yin_type == VEC_SEQ_ID) {
903:     ISCreateStride(comm,ctx->to_n,0,1,&iy);
904:     tiy  = iy;
905:   } else if (!iy && yin_type == VEC_MPI_ID) {
906:     if (xin_type == VEC_MPI_ID) {
907:       PetscInt ntmp, low;
908:       VecGetLocalSize(yin,&ntmp);
909:       VecGetOwnershipRange(yin,&low,PETSC_NULL);
910:       ISCreateStride(comm,ntmp,low,1,&iy);
911:     } else{
912:       PetscInt Ntmp;
913:       VecGetSize(yin,&Ntmp);
914:       ISCreateStride(comm,Ntmp,0,1,&iy);
915:     }
916:     tiy  = iy;
917:   } else if (!iy) {
918:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"iy not given, but not Seq or MPI vector");
919:   }

921:   /*
922:      Determine types of index sets
923:   */
924:   PetscObjectTypeCompare((PetscObject)ix,ISBLOCK,&flag);
925:   if (flag) ix_type = IS_BLOCK_ID;
926:   PetscObjectTypeCompare((PetscObject)iy,ISBLOCK,&flag);
927:   if (flag) iy_type = IS_BLOCK_ID;
928:   PetscObjectTypeCompare((PetscObject)ix,ISSTRIDE,&flag);
929:   if (flag) ix_type = IS_STRIDE_ID;
930:   PetscObjectTypeCompare((PetscObject)iy,ISSTRIDE,&flag);
931:   if (flag) iy_type = IS_STRIDE_ID;

933:   /* ===========================================================================================================
934:         Check for special cases
935:      ==========================================================================================================*/
936:   /* ---------------------------------------------------------------------------*/
937:   if (xin_type == VEC_SEQ_ID && yin_type == VEC_SEQ_ID) {
938:     if (ix_type == IS_GENERAL_ID && iy_type == IS_GENERAL_ID){
939:       PetscInt               nx,ny;
940:       const PetscInt         *idx,*idy;
941:       VecScatter_Seq_General *to = PETSC_NULL,*from = PETSC_NULL;

943:       ISGetLocalSize(ix,&nx);
944:       ISGetLocalSize(iy,&ny);
945:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
946:       ISGetIndices(ix,&idx);
947:       ISGetIndices(iy,&idy);
948:       PetscMalloc2(1,VecScatter_Seq_General,&to,1,VecScatter_Seq_General,&from);
949:       PetscMalloc2(nx,PetscInt,&to->vslots,nx,PetscInt,&from->vslots);
950:       to->n             = nx;
951: #if defined(PETSC_USE_DEBUG)
952:       VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
953: #endif
954:       PetscMemcpy(to->vslots,idy,nx*sizeof(PetscInt));
955:       from->n           = nx;
956: #if defined(PETSC_USE_DEBUG)
957:       VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
958: #endif
959:        PetscMemcpy(from->vslots,idx,nx*sizeof(PetscInt));
960:       to->type          = VEC_SCATTER_SEQ_GENERAL;
961:       from->type        = VEC_SCATTER_SEQ_GENERAL;
962:       ctx->todata       = (void*)to;
963:       ctx->fromdata     = (void*)from;
964:       ctx->begin        = VecScatterBegin_SGtoSG;
965:       ctx->end          = 0;
966:       ctx->destroy      = VecScatterDestroy_SGtoSG;
967:       ctx->copy         = VecScatterCopy_SGToSG;
968:       PetscInfo(xin,"Special case: sequential vector general scatter\n");
969:       goto functionend;
970:     } else if (ix_type == IS_STRIDE_ID &&  iy_type == IS_STRIDE_ID){
971:       PetscInt               nx,ny,to_first,to_step,from_first,from_step;
972:       VecScatter_Seq_Stride  *from8 = PETSC_NULL,*to8 = PETSC_NULL;

974:       ISGetLocalSize(ix,&nx);
975:       ISGetLocalSize(iy,&ny);
976:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
977:       ISStrideGetInfo(iy,&to_first,&to_step);
978:       ISStrideGetInfo(ix,&from_first,&from_step);
979:       PetscMalloc2(1,VecScatter_Seq_Stride,&to8,1,VecScatter_Seq_Stride,&from8);
980:       to8->n             = nx;
981:       to8->first         = to_first;
982:       to8->step          = to_step;
983:       from8->n           = nx;
984:       from8->first       = from_first;
985:       from8->step        = from_step;
986:       to8->type          = VEC_SCATTER_SEQ_STRIDE;
987:       from8->type        = VEC_SCATTER_SEQ_STRIDE;
988:       ctx->todata       = (void*)to8;
989:       ctx->fromdata     = (void*)from8;
990:       ctx->begin        = VecScatterBegin_SStoSS;
991:       ctx->end          = 0;
992:       ctx->destroy      = VecScatterDestroy_SStoSS;
993:       ctx->copy         = VecScatterCopy_SStoSS;
994:       PetscInfo(xin,"Special case: sequential vector stride to stride\n");
995:       goto functionend;
996:     } else if (ix_type == IS_GENERAL_ID && iy_type == IS_STRIDE_ID){
997:       PetscInt               nx,ny,first,step;
998:       const PetscInt         *idx;
999:       VecScatter_Seq_General *from9 = PETSC_NULL;
1000:       VecScatter_Seq_Stride  *to9 = PETSC_NULL;

1002:       ISGetLocalSize(ix,&nx);
1003:       ISGetIndices(ix,&idx);
1004:       ISGetLocalSize(iy,&ny);
1005:       ISStrideGetInfo(iy,&first,&step);
1006:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1007:       PetscMalloc2(1,VecScatter_Seq_Stride,&to9,1,VecScatter_Seq_General,&from9);
1008:       PetscMalloc(nx*sizeof(PetscInt),&from9->vslots);
1009:       to9->n         = nx;
1010:       to9->first     = first;
1011:       to9->step      = step;
1012:       from9->n       = nx;
1013: #if defined(PETSC_USE_DEBUG)
1014:       VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1015: #endif
1016:       PetscMemcpy(from9->vslots,idx,nx*sizeof(PetscInt));
1017:       ctx->todata    = (void*)to9; ctx->fromdata = (void*)from9;
1018:       if (step == 1)  ctx->begin = VecScatterBegin_SGtoSS_Stride1;
1019:       else            ctx->begin = VecScatterBegin_SGtoSS;
1020:       ctx->destroy   = VecScatterDestroy_SGtoSS;
1021:       ctx->end       = 0;
1022:       ctx->copy      = VecScatterCopy_SGToSS;
1023:       to9->type      = VEC_SCATTER_SEQ_STRIDE;
1024:       from9->type    = VEC_SCATTER_SEQ_GENERAL;
1025:       PetscInfo(xin,"Special case: sequential vector general to stride\n");
1026:       goto functionend;
1027:     } else if (ix_type == IS_STRIDE_ID && iy_type == IS_GENERAL_ID){
1028:       PetscInt               nx,ny,first,step;
1029:       const PetscInt         *idy;
1030:       VecScatter_Seq_General *to10 = PETSC_NULL;
1031:       VecScatter_Seq_Stride  *from10 = PETSC_NULL;

1033:       ISGetLocalSize(ix,&nx);
1034:       ISGetIndices(iy,&idy);
1035:       ISGetLocalSize(iy,&ny);
1036:       ISStrideGetInfo(ix,&first,&step);
1037:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1038:       PetscMalloc2(1,VecScatter_Seq_General,&to10,1,VecScatter_Seq_Stride,&from10);
1039:       PetscMalloc(nx*sizeof(PetscInt),&to10->vslots);
1040:       from10->n         = nx;
1041:       from10->first     = first;
1042:       from10->step      = step;
1043:       to10->n           = nx;
1044: #if defined(PETSC_USE_DEBUG)
1045:       VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1046: #endif
1047:       PetscMemcpy(to10->vslots,idy,nx*sizeof(PetscInt));
1048:       ctx->todata     = (void*)to10;
1049:       ctx->fromdata   = (void*)from10;
1050:       if (step == 1) ctx->begin = VecScatterBegin_SStoSG_Stride1;
1051:       else           ctx->begin = VecScatterBegin_SStoSG;
1052:       ctx->destroy    = VecScatterDestroy_SStoSG;
1053:       ctx->end        = 0;
1054:       ctx->copy       = 0;
1055:       to10->type      = VEC_SCATTER_SEQ_GENERAL;
1056:       from10->type    = VEC_SCATTER_SEQ_STRIDE;
1057:       PetscInfo(xin,"Special case: sequential vector stride to general\n");
1058:       goto functionend;
1059:     } else {
1060:       PetscInt               nx,ny;
1061:       const PetscInt         *idx,*idy;
1062:       VecScatter_Seq_General *to11 = PETSC_NULL,*from11 = PETSC_NULL;
1063:       PetscBool              idnx,idny;

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

1069:       ISIdentity(ix,&idnx);
1070:       ISIdentity(iy,&idny);
1071:       if (idnx && idny) {
1072:         VecScatter_Seq_Stride *to13 = PETSC_NULL,*from13 = PETSC_NULL;
1073:         PetscMalloc2(1,VecScatter_Seq_Stride,&to13,1,VecScatter_Seq_Stride,&from13);
1074:         to13->n           = nx;
1075:         to13->first       = 0;
1076:         to13->step        = 1;
1077:         from13->n         = nx;
1078:         from13->first     = 0;
1079:         from13->step      = 1;
1080:         to13->type        = VEC_SCATTER_SEQ_STRIDE;
1081:         from13->type      = VEC_SCATTER_SEQ_STRIDE;
1082:         ctx->todata       = (void*)to13;
1083:         ctx->fromdata     = (void*)from13;
1084:         ctx->begin        = VecScatterBegin_SStoSS;
1085:         ctx->end          = 0;
1086:         ctx->destroy      = VecScatterDestroy_SStoSS;
1087:         ctx->copy         = VecScatterCopy_SStoSS;
1088:         PetscInfo(xin,"Special case: sequential copy\n");
1089:         goto functionend;
1090:       }

1092:       ISGetIndices(iy,&idy);
1093:       ISGetIndices(ix,&idx);
1094:       PetscMalloc2(1,VecScatter_Seq_General,&to11,1,VecScatter_Seq_General,&from11);
1095:       PetscMalloc2(nx,PetscInt,&to11->vslots,nx,PetscInt,&from11->vslots);
1096:       to11->n           = nx;
1097: #if defined(PETSC_USE_DEBUG)
1098:       VecScatterCheckIndices_Private(ctx->to_n,ny,idy);
1099: #endif
1100:       PetscMemcpy(to11->vslots,idy,nx*sizeof(PetscInt));
1101:       from11->n         = nx;
1102: #if defined(PETSC_USE_DEBUG)
1103:       VecScatterCheckIndices_Private(ctx->from_n,nx,idx);
1104: #endif
1105:       PetscMemcpy(from11->vslots,idx,nx*sizeof(PetscInt));
1106:       to11->type        = VEC_SCATTER_SEQ_GENERAL;
1107:       from11->type      = VEC_SCATTER_SEQ_GENERAL;
1108:       ctx->todata       = (void*)to11;
1109:       ctx->fromdata     = (void*)from11;
1110:       ctx->begin        = VecScatterBegin_SGtoSG;
1111:       ctx->end          = 0;
1112:       ctx->destroy      = VecScatterDestroy_SGtoSG;
1113:       ctx->copy         = VecScatterCopy_SGToSG;
1114:       ISRestoreIndices(ix,&idx);
1115:       ISRestoreIndices(iy,&idy);
1116:       PetscInfo(xin,"Sequential vector scatter with block indices\n");
1117:       goto functionend;
1118:     }
1119:   }
1120:   /* ---------------------------------------------------------------------------*/
1121:   if (xin_type == VEC_MPI_ID && yin_type == VEC_SEQ_ID) {

1123:   /* ===========================================================================================================
1124:         Check for special cases
1125:      ==========================================================================================================*/
1126:     islocal = PETSC_FALSE;
1127:     /* special case extracting (subset of) local portion */
1128:     if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID){
1129:       PetscInt              nx,ny,to_first,to_step,from_first,from_step;
1130:       PetscInt              start,end;
1131:       VecScatter_Seq_Stride *from12 = PETSC_NULL,*to12 = PETSC_NULL;

1133:       VecGetOwnershipRange(xin,&start,&end);
1134:       ISGetLocalSize(ix,&nx);
1135:       ISStrideGetInfo(ix,&from_first,&from_step);
1136:       ISGetLocalSize(iy,&ny);
1137:       ISStrideGetInfo(iy,&to_first,&to_step);
1138:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1139:       if (ix->min >= start && ix->max < end) islocal = PETSC_TRUE; else islocal = PETSC_FALSE;
1140:       MPI_Allreduce(&islocal,&cando,1,MPI_INT,MPI_LAND,((PetscObject)xin)->comm);
1141:       if (cando) {
1142:         PetscMalloc2(1,VecScatter_Seq_Stride,&to12,1,VecScatter_Seq_Stride,&from12);
1143:         to12->n             = nx;
1144:         to12->first         = to_first;
1145:         to12->step          = to_step;
1146:         from12->n           = nx;
1147:         from12->first       = from_first-start;
1148:         from12->step        = from_step;
1149:         to12->type          = VEC_SCATTER_SEQ_STRIDE;
1150:         from12->type        = VEC_SCATTER_SEQ_STRIDE;
1151:         ctx->todata         = (void*)to12;
1152:         ctx->fromdata       = (void*)from12;
1153:         ctx->begin          = VecScatterBegin_SStoSS;
1154:         ctx->end            = 0;
1155:         ctx->destroy        = VecScatterDestroy_SStoSS;
1156:         ctx->copy           = VecScatterCopy_SStoSS;
1157:         PetscInfo(xin,"Special case: processors only getting local values\n");
1158:         goto functionend;
1159:       }
1160:     } else {
1161:       MPI_Allreduce(&islocal,&cando,1,MPI_INT,MPI_LAND,((PetscObject)xin)->comm);
1162:     }

1164:     /* test for special case of all processors getting entire vector */
1165:     /* contains check that PetscMPIInt can handle the sizes needed */
1166:     totalv = 0;
1167:     if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID){
1168:       PetscInt             i,nx,ny,to_first,to_step,from_first,from_step,N;
1169:       PetscMPIInt          *count = PETSC_NULL,*displx;
1170:       VecScatter_MPI_ToAll *sto = PETSC_NULL;

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

1185: #if defined(PETSC_USE_64BIT_INDICES)
1186:       if (cando && (yin->map->N < PETSC_MPI_INT_MAX)) {
1187: #else
1188:       if (cando) {
1189: #endif
1190:         MPI_Comm_size(((PetscObject)ctx)->comm,&size);
1191:         PetscMalloc3(1,VecScatter_MPI_ToAll,&sto,size,PetscMPIInt,&count,size,PetscMPIInt,&displx);
1192:         range = xin->map->range;
1193:         for (i=0; i<size; i++) {
1194:           count[i]  = PetscMPIIntCast(range[i+1] - range[i]);
1195:           displx[i] = PetscMPIIntCast(range[i]);
1196:         }
1197:         sto->count        = count;
1198:         sto->displx       = displx;
1199:         sto->work1        = 0;
1200:         sto->work2        = 0;
1201:         sto->type         = VEC_SCATTER_MPI_TOALL;
1202:         ctx->todata       = (void*)sto;
1203:         ctx->fromdata     = 0;
1204:         ctx->begin        = VecScatterBegin_MPI_ToAll;
1205:         ctx->end          = 0;
1206:         ctx->destroy      = VecScatterDestroy_MPI_ToAll;
1207:         ctx->copy         = VecScatterCopy_MPI_ToAll;
1208:         PetscInfo(xin,"Special case: all processors get entire parallel vector\n");
1209:         goto functionend;
1210:       }
1211:     } else {
1212:       MPI_Allreduce(&totalv,&cando,1,MPI_INT,MPI_LAND,((PetscObject)xin)->comm);
1213:     }

1215:     /* test for special case of processor 0 getting entire vector */
1216:     /* contains check that PetscMPIInt can handle the sizes needed */
1217:     totalv = 0;
1218:     if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID){
1219:       PetscInt             i,nx,ny,to_first,to_step,from_first,from_step,N;
1220:       PetscMPIInt          rank,*count = PETSC_NULL,*displx;
1221:       VecScatter_MPI_ToAll *sto = PETSC_NULL;

1223:       PetscObjectGetComm((PetscObject)xin,&comm);
1224:       MPI_Comm_rank(comm,&rank);
1225:       ISGetLocalSize(ix,&nx);
1226:       ISStrideGetInfo(ix,&from_first,&from_step);
1227:       ISGetLocalSize(iy,&ny);
1228:       ISStrideGetInfo(iy,&to_first,&to_step);
1229:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1230:       if (!rank) {
1231:         VecGetSize(xin,&N);
1232:         if (nx != N) {
1233:           totalv = 0;
1234:         } else if (from_first == 0        && from_step == 1 &&
1235:                    from_first == to_first && from_step == to_step){
1236:           totalv = 1;
1237:         } else totalv = 0;
1238:       } else {
1239:         if (!nx) totalv = 1;
1240:         else     totalv = 0;
1241:       }
1242:       MPI_Allreduce(&totalv,&cando,1,MPI_INT,MPI_LAND,((PetscObject)xin)->comm);

1244: #if defined(PETSC_USE_64BIT_INDICES)
1245:       if (cando && (yin->map->N < PETSC_MPI_INT_MAX)) {
1246: #else
1247:       if (cando) {
1248: #endif
1249:         MPI_Comm_size(((PetscObject)ctx)->comm,&size);
1250:         PetscMalloc3(1,VecScatter_MPI_ToAll,&sto,size,PetscMPIInt,&count,size,PetscMPIInt,&displx);
1251:         range = xin->map->range;
1252:         for (i=0; i<size; i++) {
1253:           count[i] = PetscMPIIntCast(range[i+1] - range[i]);
1254:           displx[i] = PetscMPIIntCast(range[i]);
1255:         }
1256:         sto->count        = count;
1257:         sto->displx       = displx;
1258:         sto->work1        = 0;
1259:         sto->work2        = 0;
1260:         sto->type         = VEC_SCATTER_MPI_TOONE;
1261:         ctx->todata       = (void*)sto;
1262:         ctx->fromdata     = 0;
1263:         ctx->begin        = VecScatterBegin_MPI_ToOne;
1264:         ctx->end          = 0;
1265:         ctx->destroy      = VecScatterDestroy_MPI_ToAll;
1266:         ctx->copy         = VecScatterCopy_MPI_ToAll;
1267:         PetscInfo(xin,"Special case: processor zero gets entire parallel vector, rest get none\n");
1268:         goto functionend;
1269:       }
1270:     } else {
1271:       MPI_Allreduce(&totalv,&cando,1,MPI_INT,MPI_LAND,((PetscObject)xin)->comm);
1272:     }

1274:     PetscObjectTypeCompare((PetscObject)ix,ISBLOCK,&ixblock);
1275:     PetscObjectTypeCompare((PetscObject)iy,ISBLOCK,&iyblock);
1276:     PetscObjectTypeCompare((PetscObject)iy,ISSTRIDE,&iystride);
1277:     if (ixblock) {
1278:       /* special case block to block */
1279:       if (iyblock) {
1280:         PetscInt       nx,ny,bsx,bsy;
1281:         const PetscInt *idx,*idy;
1282:         ISGetBlockSize(iy,&bsy);
1283:         ISGetBlockSize(ix,&bsx);
1284:         if (bsx == bsy && VecScatterOptimizedBS(bsx)) {
1285:           ISBlockGetLocalSize(ix,&nx);
1286:           ISBlockGetIndices(ix,&idx);
1287:           ISBlockGetLocalSize(iy,&ny);
1288:           ISBlockGetIndices(iy,&idy);
1289:           if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1290:           VecScatterCreate_PtoS(nx,idx,ny,idy,xin,yin,bsx,ctx);
1291:           ISBlockRestoreIndices(ix,&idx);
1292:           ISBlockRestoreIndices(iy,&idy);
1293:           PetscInfo(xin,"Special case: blocked indices\n");
1294:           goto functionend;
1295:         }
1296:       /* special case block to stride */
1297:       } else if (iystride) {
1298:         PetscInt ystart,ystride,ysize,bsx;
1299:         ISStrideGetInfo(iy,&ystart,&ystride);
1300:         ISGetLocalSize(iy,&ysize);
1301:         ISGetBlockSize(ix,&bsx);
1302:         /* see if stride index set is equivalent to block index set */
1303:         if (VecScatterOptimizedBS(bsx) && ((ystart % bsx) == 0) && (ystride == 1) && ((ysize % bsx) == 0)) {
1304:           PetscInt       nx,il,*idy;
1305:           const PetscInt *idx;
1306:           ISBlockGetLocalSize(ix,&nx);
1307:           ISBlockGetIndices(ix,&idx);
1308:           if (ysize != bsx*nx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1309:           PetscMalloc(nx*sizeof(PetscInt),&idy);
1310:           if (nx) {
1311:             idy[0] = ystart/bsx;
1312:             for (il=1; il<nx; il++) idy[il] = idy[il-1] + 1;
1313:           }
1314:           VecScatterCreate_PtoS(nx,idx,nx,idy,xin,yin,bsx,ctx);
1315:           PetscFree(idy);
1316:           ISBlockRestoreIndices(ix,&idx);
1317:           PetscInfo(xin,"Special case: blocked indices to stride\n");
1318:           goto functionend;
1319:         }
1320:       }
1321:     }
1322:     /* left over general case */
1323:     {
1324:       PetscInt       nx,ny;
1325:       const PetscInt *idx,*idy;
1326:       ISGetLocalSize(ix,&nx);
1327:       ISGetIndices(ix,&idx);
1328:       ISGetLocalSize(iy,&ny);
1329:       ISGetIndices(iy,&idy);
1330:       if (nx != ny) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match (%d %d)",nx,ny);
1331:       VecScatterCreate_PtoS(nx,idx,ny,idy,xin,yin,1,ctx);
1332:       ISRestoreIndices(ix,&idx);
1333:       ISRestoreIndices(iy,&idy);
1334:       PetscInfo(xin,"General case: MPI to Seq\n");
1335:       goto functionend;
1336:     }
1337:   }
1338:   /* ---------------------------------------------------------------------------*/
1339:   if (xin_type == VEC_SEQ_ID && yin_type == VEC_MPI_ID) {
1340:   /* ===========================================================================================================
1341:         Check for special cases
1342:      ==========================================================================================================*/
1343:     /* special case local copy portion */
1344:     islocal = PETSC_FALSE;
1345:     if (ix_type == IS_STRIDE_ID && iy_type == IS_STRIDE_ID){
1346:       PetscInt              nx,ny,to_first,to_step,from_step,start,end,from_first;
1347:       VecScatter_Seq_Stride *from = PETSC_NULL,*to = PETSC_NULL;

1349:       VecGetOwnershipRange(yin,&start,&end);
1350:       ISGetLocalSize(ix,&nx);
1351:       ISStrideGetInfo(ix,&from_first,&from_step);
1352:       ISGetLocalSize(iy,&ny);
1353:       ISStrideGetInfo(iy,&to_first,&to_step);
1354:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1355:       if (iy->min >= start && iy->max < end) islocal = PETSC_TRUE; else islocal = PETSC_FALSE;
1356:       MPI_Allreduce(&islocal,&cando,1,MPI_INT,MPI_LAND,((PetscObject)yin)->comm);
1357:       if (cando) {
1358:         PetscMalloc2(1,VecScatter_Seq_Stride,&to,1,VecScatter_Seq_Stride,&from);
1359:         to->n             = nx;
1360:         to->first         = to_first-start;
1361:         to->step          = to_step;
1362:         from->n           = nx;
1363:         from->first       = from_first;
1364:         from->step        = from_step;
1365:         to->type          = VEC_SCATTER_SEQ_STRIDE;
1366:         from->type        = VEC_SCATTER_SEQ_STRIDE;
1367:         ctx->todata       = (void*)to;
1368:         ctx->fromdata     = (void*)from;
1369:         ctx->begin        = VecScatterBegin_SStoSS;
1370:         ctx->end          = 0;
1371:         ctx->destroy      = VecScatterDestroy_SStoSS;
1372:         ctx->copy         = VecScatterCopy_SStoSS;
1373:         PetscInfo(xin,"Special case: sequential stride to MPI stride\n");
1374:         goto functionend;
1375:       }
1376:     } else {
1377:       MPI_Allreduce(&islocal,&cando,1,MPI_INT,MPI_LAND,((PetscObject)yin)->comm);
1378:     }
1379:       /* special case block to stride */
1380:     if (ix_type == IS_BLOCK_ID && iy_type == IS_STRIDE_ID){
1381:       PetscInt ystart,ystride,ysize,bsx;
1382:       ISStrideGetInfo(iy,&ystart,&ystride);
1383:       ISGetLocalSize(iy,&ysize);
1384:       ISGetBlockSize(ix,&bsx);
1385:       /* see if stride index set is equivalent to block index set */
1386:       if (VecScatterOptimizedBS(bsx) && ((ystart % bsx) == 0) && (ystride == 1) && ((ysize % bsx) == 0)) {
1387:         PetscInt       nx,il,*idy;
1388:         const PetscInt *idx;
1389:         ISBlockGetLocalSize(ix,&nx);
1390:         ISBlockGetIndices(ix,&idx);
1391:         if (ysize != bsx*nx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1392:         PetscMalloc(nx*sizeof(PetscInt),&idy);
1393:         if (nx) {
1394:           idy[0] = ystart/bsx;
1395:           for (il=1; il<nx; il++) idy[il] = idy[il-1] + 1;
1396:         }
1397:         VecScatterCreate_StoP(nx,idx,nx,idy,xin,yin,bsx,ctx);
1398:         PetscFree(idy);
1399:         ISBlockRestoreIndices(ix,&idx);
1400:         PetscInfo(xin,"Special case: Blocked indices to stride\n");
1401:         goto functionend;
1402:       }
1403:     }

1405:     /* general case */
1406:     {
1407:       PetscInt       nx,ny;
1408:       const PetscInt *idx,*idy;
1409:       ISGetLocalSize(ix,&nx);
1410:       ISGetIndices(ix,&idx);
1411:       ISGetLocalSize(iy,&ny);
1412:       ISGetIndices(iy,&idy);
1413:       if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1414:       VecScatterCreate_StoP(nx,idx,ny,idy,xin,yin,1,ctx);
1415:       ISRestoreIndices(ix,&idx);
1416:       ISRestoreIndices(iy,&idy);
1417:       PetscInfo(xin,"General case: Seq to MPI\n");
1418:       goto functionend;
1419:     }
1420:   }
1421:   /* ---------------------------------------------------------------------------*/
1422:   if (xin_type == VEC_MPI_ID && yin_type == VEC_MPI_ID) {
1423:     /* no special cases for now */
1424:     PetscInt       nx,ny;
1425:     const PetscInt *idx,*idy;
1426:     ISGetLocalSize(ix,&nx);
1427:     ISGetIndices(ix,&idx);
1428:     ISGetLocalSize(iy,&ny);
1429:     ISGetIndices(iy,&idy);
1430:     if (nx != ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
1431:     VecScatterCreate_PtoP(nx,idx,ny,idy,xin,yin,1,ctx);
1432:     ISRestoreIndices(ix,&idx);
1433:     ISRestoreIndices(iy,&idy);
1434:     PetscInfo(xin,"General case: MPI to MPI\n");
1435:     goto functionend;
1436:   }

1438:   functionend:
1439:   *newctx = ctx;
1440:   ISDestroy(&tix);
1441:   ISDestroy(&tiy);
1442:   flag = PETSC_FALSE;
1443:   PetscOptionsGetBool(PETSC_NULL,"-vecscatter_view_info",&flag,PETSC_NULL);
1444:   if (flag) {
1445:     PetscViewer viewer;
1446:     PetscViewerASCIIGetStdout(comm,&viewer);
1447:     PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1448:     VecScatterView(ctx,viewer);
1449:     PetscViewerPopFormat(viewer);
1450:   }
1451:   flag = PETSC_FALSE;
1452:   PetscOptionsGetBool(PETSC_NULL,"-vecscatter_view",&flag,PETSC_NULL);
1453:   if (flag) {
1454:     PetscViewer viewer;
1455:     PetscViewerASCIIGetStdout(comm,&viewer);
1456:     VecScatterView(ctx,viewer);
1457:   }
1458:   return(0);
1459: }

1461: /* ------------------------------------------------------------------*/
1464: /*@
1465:    VecScatterGetMerged - Returns true if the scatter is completed in the VecScatterBegin()
1466:       and the VecScatterEnd() does nothing

1468:    Not Collective

1470:    Input Parameter:
1471: .   ctx - scatter context created with VecScatterCreate()

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

1476:    Level: developer

1478: .seealso: VecScatterCreate(), VecScatterEnd(), VecScatterBegin()
1479: @*/
1480: PetscErrorCode  VecScatterGetMerged(VecScatter ctx,PetscBool  *flg)
1481: {
1484:   *flg = ctx->beginandendtogether;
1485:   return(0);
1486: }

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

1494:    Neighbor-wise Collective on VecScatter and Vec

1496:    Input Parameters:
1497: +  inctx - scatter context generated by VecScatterCreate()
1498: .  x - the vector from which we scatter
1499: .  y - the vector to which we scatter
1500: .  addv - either ADD_VALUES or INSERT_VALUES, with INSERT_VALUES mode any location 
1501:           not scattered to retains its old value; i.e. the vector is NOT first zeroed.
1502: -  mode - the scattering mode, usually SCATTER_FORWARD.  The available modes are:
1503:     SCATTER_FORWARD or SCATTER_REVERSE


1506:    Level: intermediate

1508:    Options Database: See VecScatterCreate()

1510:    Notes:
1511:    The vectors x and y need not be the same vectors used in the call 
1512:    to VecScatterCreate(), but x must have the same parallel data layout
1513:    as that passed in as the x to VecScatterCreate(), similarly for the y.
1514:    Most likely they have been obtained from VecDuplicate().

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

1519:    If you use SCATTER_REVERSE the two arguments x and y should be reversed, from 
1520:    the SCATTER_FORWARD.
1521:    
1522:    y[iy[i]] = x[ix[i]], for i=0,...,ni-1

1524:    This scatter is far more general than the conventional
1525:    scatter, since it can be a gather or a scatter or a combination,
1526:    depending on the indices ix and iy.  If x is a parallel vector and y
1527:    is sequential, VecScatterBegin() can serve to gather values to a
1528:    single processor.  Similarly, if y is parallel and x sequential, the
1529:    routine can scatter from one processor to many processors.

1531:    Concepts: scatter^between vectors
1532:    Concepts: gather^between vectors

1534: .seealso: VecScatterCreate(), VecScatterEnd()
1535: @*/
1536: PetscErrorCode  VecScatterBegin(VecScatter inctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
1537: {
1539: #if defined(PETSC_USE_DEBUG)
1540:   PetscInt      to_n,from_n;
1541: #endif
1546:   if (inctx->inuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE," Scatter ctx already in use");
1547:   CHKMEMQ;

1549: #if defined(PETSC_USE_DEBUG)
1550:   /*
1551:      Error checking to make sure these vectors match the vectors used
1552:    to create the vector scatter context. -1 in the from_n and to_n indicate the
1553:    vector lengths are unknown (for example with mapped scatters) and thus 
1554:    no error checking is performed.
1555:   */
1556:   if (inctx->from_n >= 0 && inctx->to_n >= 0) {
1557:     VecGetLocalSize(x,&from_n);
1558:     VecGetLocalSize(y,&to_n);
1559:     if (mode & SCATTER_REVERSE) {
1560:       if (to_n != inctx->from_n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter reverse and vector to != ctx from size)",to_n,inctx->from_n);
1561:       if (from_n != inctx->to_n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter reverse and vector from != ctx to size)",from_n,inctx->to_n);
1562:     } else {
1563:       if (to_n != inctx->to_n)     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter forward and vector to != ctx to size)",to_n,inctx->to_n);
1564:       if (from_n != inctx->from_n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector wrong size %D for scatter %D (scatter forward and vector from != ctx from size)",from_n,inctx->from_n);
1565:     }
1566:   }
1567: #endif

1569:   inctx->inuse = PETSC_TRUE;
1570:   PetscLogEventBarrierBegin(VEC_ScatterBarrier,0,0,0,0,((PetscObject)inctx)->comm);
1571:   (*inctx->begin)(inctx,x,y,addv,mode);
1572:   if (inctx->beginandendtogether && inctx->end) {
1573:     inctx->inuse = PETSC_FALSE;
1574:     (*inctx->end)(inctx,x,y,addv,mode);
1575:   }
1576:   PetscLogEventBarrierEnd(VEC_ScatterBarrier,0,0,0,0,((PetscObject)inctx)->comm);
1577:   CHKMEMQ;
1578:   return(0);
1579: }

1581: /* --------------------------------------------------------------------*/
1584: /*@
1585:    VecScatterEnd - Ends a generalized scatter from one vector to another.  Call
1586:    after first calling VecScatterBegin().

1588:    Neighbor-wise Collective on VecScatter and Vec

1590:    Input Parameters:
1591: +  ctx - scatter context generated by VecScatterCreate()
1592: .  x - the vector from which we scatter
1593: .  y - the vector to which we scatter
1594: .  addv - either ADD_VALUES or INSERT_VALUES.
1595: -  mode - the scattering mode, usually SCATTER_FORWARD.  The available modes are:
1596:      SCATTER_FORWARD, SCATTER_REVERSE

1598:    Level: intermediate

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

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

1605: .seealso: VecScatterBegin(), VecScatterCreate()
1606: @*/
1607: PetscErrorCode  VecScatterEnd(VecScatter ctx,Vec x,Vec y,InsertMode addv,ScatterMode mode)
1608: {

1615:   ctx->inuse = PETSC_FALSE;
1616:   if (!ctx->end) return(0);
1617:   CHKMEMQ;
1618:   if (!ctx->beginandendtogether) {
1619:     PetscLogEventBegin(VEC_ScatterEnd,ctx,x,y,0);
1620:     (*(ctx)->end)(ctx,x,y,addv,mode);
1621:     PetscLogEventEnd(VEC_ScatterEnd,ctx,x,y,0);
1622:   }
1623:   CHKMEMQ;
1624:   return(0);
1625: }

1629: /*@C
1630:    VecScatterDestroy - Destroys a scatter context created by 
1631:    VecScatterCreate().

1633:    Collective on VecScatter

1635:    Input Parameter:
1636: .  ctx - the scatter context

1638:    Level: intermediate

1640: .seealso: VecScatterCreate(), VecScatterCopy()
1641: @*/
1642: PetscErrorCode  VecScatterDestroy(VecScatter *ctx)
1643: {

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

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

1655:   (*(*ctx)->destroy)(*ctx);
1656:   PetscHeaderDestroy(ctx);
1657:   return(0);
1658: }

1662: /*@
1663:    VecScatterCopy - Makes a copy of a scatter context.

1665:    Collective on VecScatter

1667:    Input Parameter:
1668: .  sctx - the scatter context

1670:    Output Parameter:
1671: .  ctx - the context copy

1673:    Level: advanced

1675: .seealso: VecScatterCreate(), VecScatterDestroy()
1676: @*/
1677: PetscErrorCode  VecScatterCopy(VecScatter sctx,VecScatter *ctx)
1678: {

1684:   if (!sctx->copy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot copy this type");
1685:   PetscHeaderCreate(*ctx,_p_VecScatter,int,VEC_SCATTER_CLASSID,0,"VecScatter","VecScatter","Vec",((PetscObject)sctx)->comm,VecScatterDestroy,VecScatterView);
1686:   (*ctx)->to_n   = sctx->to_n;
1687:   (*ctx)->from_n = sctx->from_n;
1688:   (*sctx->copy)(sctx,*ctx);
1689:   return(0);
1690: }


1693: /* ------------------------------------------------------------------*/
1696: /*@
1697:    VecScatterView - Views a vector scatter context.

1699:    Collective on VecScatter

1701:    Input Parameters:
1702: +  ctx - the scatter context
1703: -  viewer - the viewer for displaying the context

1705:    Level: intermediate

1707: @*/
1708: PetscErrorCode  VecScatterView(VecScatter ctx,PetscViewer viewer)
1709: {

1714:   if (!viewer) {
1715:     PetscViewerASCIIGetStdout(((PetscObject)ctx)->comm,&viewer);
1716:   }
1718:   if (ctx->view) {
1719:     (*ctx->view)(ctx,viewer);
1720:   }
1721:   return(0);
1722: }

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

1730:    Collective on VecScatter

1732:    Input Parameters:
1733: +  scat - vector scatter context
1734: .  from - remapping for "from" indices (may be PETSC_NULL)
1735: -  to   - remapping for "to" indices (may be PETSC_NULL)

1737:    Level: developer

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

1743:           In the sequential case the todata is the indices where the 
1744:           data is put and the fromdata is where it is taken from.
1745:           This is backwards from the paralllel case! CRY! CRY! CRY!

1747: @*/
1748: PetscErrorCode  VecScatterRemap(VecScatter scat,PetscInt *rto,PetscInt *rfrom)
1749: {
1750:   VecScatter_Seq_General *to,*from;
1751:   VecScatter_MPI_General *mto;
1752:   PetscInt               i;


1759:   from = (VecScatter_Seq_General *)scat->fromdata;
1760:   mto  = (VecScatter_MPI_General *)scat->todata;

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

1764:   if (rto) {
1765:     if (mto->type == VEC_SCATTER_MPI_GENERAL) {
1766:       /* handle off processor parts */
1767:       for (i=0; i<mto->starts[mto->n]; i++) {
1768:         mto->indices[i] = rto[mto->indices[i]];
1769:       }
1770:       /* handle local part */
1771:       to = &mto->local;
1772:       for (i=0; i<to->n; i++) {
1773:         to->vslots[i] = rto[to->vslots[i]];
1774:       }
1775:     } else if (from->type == VEC_SCATTER_SEQ_GENERAL) {
1776:       for (i=0; i<from->n; i++) {
1777:         from->vslots[i] = rto[from->vslots[i]];
1778:       }
1779:     } else if (from->type == VEC_SCATTER_SEQ_STRIDE) {
1780:       VecScatter_Seq_Stride *sto = (VecScatter_Seq_Stride*)from;
1781: 
1782:       /* if the remapping is the identity and stride is identity then skip remap */
1783:       if (sto->step == 1 && sto->first == 0) {
1784:         for (i=0; i<sto->n; i++) {
1785:           if (rto[i] != i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1786:         }
1787:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1788:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
1789:   }

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

1793:   /*
1794:      Mark then vector lengths as unknown because we do not know the 
1795:    lengths of the remapped vectors
1796:   */
1797:   scat->from_n = -1;
1798:   scat->to_n   = -1;

1800:   return(0);
1801: }