Actual source code: vscat.c

      1: #define PETSCVEC_DLL
      2: /*
      3:      Code for creating scatters between vectors. This file
      4:   includes the code for scattering between sequential vectors and
      5:   some special cases for parallel scatters.
      6: */

     8:  #include src/vec/is/isimpl.h
     9:  #include vecimpl.h

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

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

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

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

     37:    This code was written by Cameron Cooper, Occidental College, Fall 1995,
     38:    will working at ANL as a SERS student.
     39: */
     42: PetscErrorCode VecScatterBegin_MPI_ToAll(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
     43: {
     44: #if defined(PETSC_USE_64BIT_INT)
     46:   SETERRQ(PETSC_ERR_SUP,"Not implemented due to limited MPI support for extremely long gathers");
     47: #else
     49:   PetscInt       yy_n,xx_n,*range;
     50:   PetscScalar    *xv,*yv;
     51:   PetscMap       map;

     54:   VecGetLocalSize(y,&yy_n);
     55:   VecGetLocalSize(x,&xx_n);
     56:   VecGetArray(y,&yv);
     57:   if (x != y) {VecGetArray(x,&xv);}

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

     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:         PetscMalloc(xx_n*sizeof(PetscScalar),&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:           PetscMalloc(xx_n*sizeof(PetscScalar),& xvt2);
     87:           scat->work2 = xvt2;
     88:         }
     89:         VecGetPetscMap(y,&map);
     90:         PetscMapGetGlobalRange(map,&range);
     91:         MPI_Gatherv(yv,yy_n,MPIU_SCALAR,xvt2,scat->count,range,MPIU_SCALAR,0,ctx->comm);
     92: #if defined(PETSC_USE_COMPLEX)
     93:         MPI_Reduce(xv,xvt,2*xx_n,MPIU_REAL,MPI_SUM,0,ctx->comm);
     94: #else
     95:         MPI_Reduce(xv,xvt,xx_n,MPIU_SCALAR,MPI_SUM,0,ctx->comm);
     96: #endif
     97:         if (addv == ADD_VALUES) {
     98:           for (i=0; i<xx_n; i++) {
     99:             xvt[i] += xvt2[i];
    100:           }
    101: #if !defined(PETSC_USE_COMPLEX)
    102:         } else if (addv == MAX_VALUES) {
    103:           for (i=0; i<xx_n; i++) {
    104:             xvt[i] = PetscMax(xvt[i],xvt2[i]);
    105:           }
    106: #endif
    107:         } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
    108:         MPI_Scatterv(xvt,scat->count,map->range,MPIU_SCALAR,yv,yy_n,MPIU_SCALAR,0,ctx->comm);
    109:       } else {
    110:         VecGetPetscMap(y,&map);
    111:         PetscMapGetGlobalRange(map,&range);
    112:         MPI_Gatherv(yv,yy_n,MPIU_SCALAR,0, 0,0,MPIU_SCALAR,0,ctx->comm);
    113: #if defined(PETSC_USE_COMPLEX)
    114:         MPI_Reduce(xv,xvt,2*xx_n,MPIU_REAL,MPI_SUM,0,ctx->comm);
    115: #else
    116:         MPI_Reduce(xv,xvt,xx_n,MPIU_SCALAR,MPI_SUM,0,ctx->comm);
    117: #endif
    118:         MPI_Scatterv(0,scat->count,range,MPIU_SCALAR,yv,yy_n,MPIU_SCALAR,0,ctx->comm);
    119:       }
    120:     }
    121:   } else {
    122:     PetscScalar          *yvt;
    123:     VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
    124:     PetscInt             i;

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

    156: /*
    157:       This is special scatter code for when the entire parallel vector is
    158:    copied to processor 0.

    160: */
    163: PetscErrorCode VecScatterBegin_MPI_ToOne(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
    164: {
    165: #if defined(PETSC_USE_64BIT_INT)
    167:   SETERRQ(PETSC_ERR_SUP,"Not implemented due to limited MPI support for extremely long gathers");
    168: #else
    170:   PetscMPIInt    rank;
    171:   PetscInt       yy_n,xx_n,*range;
    172:   PetscScalar    *xv,*yv;
    173:   MPI_Comm       comm;
    174:   PetscMap       map;

    177:   VecGetLocalSize(y,&yy_n);
    178:   VecGetLocalSize(x,&xx_n);
    179:   VecGetArray(x,&xv);
    180:   VecGetArray(y,&yv);

    182:   PetscObjectGetComm((PetscObject)x,&comm);
    183:   MPI_Comm_rank(comm,&rank);

    185:   /* --------  Reverse scatter; spread from processor 0 to other processors */
    186:   if (mode & SCATTER_REVERSE) {
    187:     PetscScalar          *yvt;
    188:     VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
    189:     PetscInt             i;

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

    220:     VecGetPetscMap(x,&map);
    221:     PetscMapGetGlobalRange(map,&range);
    222:     if (addv == INSERT_VALUES) {
    223:       MPI_Gatherv(xv,xx_n,MPIU_SCALAR,yv,scat->count,range,MPIU_SCALAR,0,ctx->comm);
    224:     } else {
    225:       if (!rank) {
    226:         if (scat->work1) yvt = scat->work1;
    227:         else {
    228:           PetscMalloc(yy_n*sizeof(PetscScalar),&yvt);
    229:           scat->work1 = yvt;
    230:         }
    231:       }
    232:       MPI_Gatherv(xv,xx_n,MPIU_SCALAR,yvt,scat->count,range,MPIU_SCALAR,0,ctx->comm);
    233:       if (!rank) {
    234:         if (addv == ADD_VALUES) {
    235:           for (i=0; i<yy_n; i++) {
    236:             yv[i] += yvt[i];
    237:           }
    238: #if !defined(PETSC_USE_COMPLEX)
    239:         } else if (addv == MAX_VALUES) {
    240:           for (i=0; i<yy_n; i++) {
    241:             yv[i] = PetscMax(yv[i],yvt[i]);
    242:           }
    243: #endif
    244:         }  else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
    245:       }
    246:     }
    247:   }
    248:   VecRestoreArray(x,&xv);
    249:   VecRestoreArray(y,&yv);
    250: #endif
    251:   return(0);
    252: }

    254: /*
    255:        The follow to are used for both VecScatterBegin_MPI_ToAll() and VecScatterBegin_MPI_ToOne()
    256: */
    259: PetscErrorCode VecScatterDestroy_MPI_ToAll(VecScatter ctx)
    260: {
    261:   VecScatter_MPI_ToAll *scat = (VecScatter_MPI_ToAll*)ctx->todata;
    262:   PetscErrorCode       ierr;

    265:   if (scat->work1) {PetscFree(scat->work1);}
    266:   if (scat->work2) {PetscFree(scat->work2);}
    267:   PetscFree2(ctx->todata,scat->count);
    268:   PetscHeaderDestroy(ctx);
    269:   return(0);
    270: }

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

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

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

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

    298: PetscErrorCode VecScatterDestroy_SStoSG(VecScatter ctx)
    299: {
    300:   PetscErrorCode         ierr;

    303:   PetscFree3(ctx->todata,((VecScatter_Seq_General*)ctx->todata)->vslots,ctx->fromdata);
    304:   PetscHeaderDestroy(ctx);
    305:   return(0);
    306: }

    310: PetscErrorCode VecScatterDestroy_SStoSS(VecScatter ctx)
    311: {

    315:   PetscFree2(ctx->todata,ctx->fromdata);
    316:   PetscHeaderDestroy(ctx);
    317:   return(0);
    318: }

    320: /* -------------------------------------------------------------------------*/
    323: PetscErrorCode VecScatterCopy_MPI_ToAll(VecScatter in,VecScatter out)
    324: {
    325:   VecScatter_MPI_ToAll *in_to = (VecScatter_MPI_ToAll*)in->todata,*sto;
    326:   PetscErrorCode       ierr;
    327:   PetscMPIInt          size;
    328:   PetscInt             i;

    331:   out->postrecvs      = 0;
    332:   out->begin          = in->begin;
    333:   out->end            = in->end;
    334:   out->copy           = in->copy;
    335:   out->destroy        = in->destroy;
    336:   out->view           = in->view;

    338:   MPI_Comm_size(out->comm,&size);
    339:   PetscMalloc2(1,VecScatter_MPI_ToAll,&sto,size,PetscMPIInt,&sto->count);
    340:   sto->type           = in_to->type;
    341:   for (i=0; i<size; i++) {
    342:     sto->count[i] = in_to->count[i];
    343:   }
    344:   sto->work1         = 0;
    345:   sto->work2         = 0;
    346:   out->todata        = (void*)sto;
    347:   out->fromdata      = (void*)0;
    348:   return(0);
    349: }

    351: /* --------------------------------------------------------------------------------------*/
    352: /*
    353:         Scatter: sequential general to sequential general
    354: */
    357: PetscErrorCode VecScatterBegin_SGtoSG(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
    358: {
    359:   VecScatter_Seq_General *gen_to = (VecScatter_Seq_General*)ctx->todata;
    360:   VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
    361:   PetscErrorCode         ierr;
    362:   PetscInt               i,n = gen_from->n,*fslots,*tslots;
    363:   PetscScalar            *xv,*yv;
    364:
    366:   VecGetArray(x,&xv);
    367:   if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}
    368:   if (mode & SCATTER_REVERSE){
    369:     gen_to   = (VecScatter_Seq_General*)ctx->fromdata;
    370:     gen_from = (VecScatter_Seq_General*)ctx->todata;
    371:   }
    372:   fslots = gen_from->vslots;
    373:   tslots = gen_to->vslots;

    375:   if (addv == INSERT_VALUES) {
    376:     for (i=0; i<n; i++) {yv[tslots[i]] = xv[fslots[i]];}
    377:   } else if (addv == ADD_VALUES) {
    378:     for (i=0; i<n; i++) {yv[tslots[i]] += xv[fslots[i]];}
    379: #if !defined(PETSC_USE_COMPLEX)
    380:   } else  if (addv == MAX_VALUES) {
    381:     for (i=0; i<n; i++) {yv[tslots[i]] = PetscMax(yv[tslots[i]],xv[fslots[i]]);}
    382: #endif
    383:   } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
    384:   VecRestoreArray(x,&xv);
    385:   if (x != y) {VecRestoreArray(y,&yv);}
    386:   return(0);
    387: }

    389: /*
    390:     Scatter: sequential general to sequential stride 1
    391: */
    394: PetscErrorCode VecScatterBegin_SGtoSS_Stride1(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
    395: {
    396:   VecScatter_Seq_Stride  *gen_to   = (VecScatter_Seq_Stride*)ctx->todata;
    397:   VecScatter_Seq_General *gen_from = (VecScatter_Seq_General*)ctx->fromdata;
    398:   PetscInt               i,n = gen_from->n,*fslots = gen_from->vslots;
    399:   PetscErrorCode         ierr;
    400:   PetscInt               first = gen_to->first;
    401:   PetscScalar            *xv,*yv;
    402:
    404:   VecGetArray(x,&xv);
    405:   if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}
    406:   if (mode & SCATTER_REVERSE){
    407:     xv += first;
    408:     if (addv == INSERT_VALUES) {
    409:       for (i=0; i<n; i++) {yv[fslots[i]] = xv[i];}
    410:     } else  if (addv == ADD_VALUES) {
    411:       for (i=0; i<n; i++) {yv[fslots[i]] += xv[i];}
    412: #if !defined(PETSC_USE_COMPLEX)
    413:     } else  if (addv == MAX_VALUES) {
    414:       for (i=0; i<n; i++) {yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[i]);}
    415: #endif
    416:     } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
    417:   } else {
    418:     yv += first;
    419:     if (addv == INSERT_VALUES) {
    420:       for (i=0; i<n; i++) {yv[i] = xv[fslots[i]];}
    421:     } else  if (addv == ADD_VALUES) {
    422:       for (i=0; i<n; i++) {yv[i] += xv[fslots[i]];}
    423: #if !defined(PETSC_USE_COMPLEX)
    424:     } else if (addv == MAX_VALUES) {
    425:       for (i=0; i<n; i++) {yv[i] = PetscMax(yv[i],xv[fslots[i]]);}
    426: #endif
    427:     } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
    428:   }
    429:   VecRestoreArray(x,&xv);
    430:   if (x != y) {VecRestoreArray(y,&yv);}
    431:   return(0);
    432: }

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

    452:   if (mode & SCATTER_REVERSE){
    453:     if (addv == INSERT_VALUES) {
    454:       for (i=0; i<n; i++) {yv[fslots[i]] = xv[first + i*step];}
    455:     } else if (addv == ADD_VALUES) {
    456:       for (i=0; i<n; i++) {yv[fslots[i]] += xv[first + i*step];}
    457: #if !defined(PETSC_USE_COMPLEX)
    458:     } else if (addv == MAX_VALUES) {
    459:       for (i=0; i<n; i++) {yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[first + i*step]);}
    460: #endif
    461:     } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
    462:   } else {
    463:     if (addv == INSERT_VALUES) {
    464:       for (i=0; i<n; i++) {yv[first + i*step] = xv[fslots[i]];}
    465:     } else if (addv == ADD_VALUES) {
    466:       for (i=0; i<n; i++) {yv[first + i*step] += xv[fslots[i]];}
    467: #if !defined(PETSC_USE_COMPLEX)
    468:     } else if (addv == MAX_VALUES) {
    469:       for (i=0; i<n; i++) {yv[first + i*step] = PetscMax(yv[first + i*step],xv[fslots[i]]);}
    470: #endif
    471:     } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
    472:   }
    473:   VecRestoreArray(x,&xv);
    474:   if (x != y) {VecRestoreArray(y,&yv);}
    475:   return(0);
    476: }

    478: /*
    479:     Scatter: sequential stride 1 to sequential general
    480: */
    483: PetscErrorCode VecScatterBegin_SStoSG_Stride1(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
    484: {
    485:   VecScatter_Seq_Stride  *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
    486:   VecScatter_Seq_General *gen_to   = (VecScatter_Seq_General*)ctx->todata;
    487:   PetscInt               i,n = gen_from->n,*fslots = gen_to->vslots;
    488:   PetscErrorCode         ierr;
    489:   PetscInt               first = gen_from->first;
    490:   PetscScalar            *xv,*yv;
    491:
    493:   VecGetArray(x,&xv);
    494:   if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}

    496:   if (mode & SCATTER_REVERSE){
    497:     yv += first;
    498:     if (addv == INSERT_VALUES) {
    499:       for (i=0; i<n; i++) {yv[i] = xv[fslots[i]];}
    500:     } else  if (addv == ADD_VALUES) {
    501:       for (i=0; i<n; i++) {yv[i] += xv[fslots[i]];}
    502: #if !defined(PETSC_USE_COMPLEX)
    503:     } else  if (addv == MAX_VALUES) {
    504:       for (i=0; i<n; i++) {yv[i] = PetscMax(yv[i],xv[fslots[i]]);}
    505: #endif
    506:     } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
    507:   } else {
    508:     xv += first;
    509:     if (addv == INSERT_VALUES) {
    510:       for (i=0; i<n; i++) {yv[fslots[i]] = xv[i];}
    511:     } else  if (addv == ADD_VALUES) {
    512:       for (i=0; i<n; i++) {yv[fslots[i]] += xv[i];}
    513: #if !defined(PETSC_USE_COMPLEX)
    514:     } else  if (addv == MAX_VALUES) {
    515:       for (i=0; i<n; i++) {yv[fslots[i]] = PetscMax(yv[fslots[i]],xv[i]);}
    516: #endif
    517:     } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
    518:   }
    519:   VecRestoreArray(x,&xv);
    520:   if (x != y) {VecRestoreArray(y,&yv);}
    521:   return(0);
    522: }

    526: /*
    527:    Scatter: sequential stride to sequential general
    528: */
    529: PetscErrorCode VecScatterBegin_SStoSG(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
    530: {
    531:   VecScatter_Seq_Stride  *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
    532:   VecScatter_Seq_General *gen_to   = (VecScatter_Seq_General*)ctx->todata;
    533:   PetscInt               i,n = gen_from->n,*fslots = gen_to->vslots;
    534:   PetscErrorCode         ierr;
    535:   PetscInt               first = gen_from->first,step = gen_from->step;
    536:   PetscScalar            *xv,*yv;
    537:
    539:   VecGetArray(x,&xv);
    540:   if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}

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

    568: /*
    569:      Scatter: sequential stride to sequential stride
    570: */
    573: PetscErrorCode VecScatterBegin_SStoSS(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
    574: {
    575:   VecScatter_Seq_Stride *gen_to   = (VecScatter_Seq_Stride*)ctx->todata;
    576:   VecScatter_Seq_Stride *gen_from = (VecScatter_Seq_Stride*)ctx->fromdata;
    577:   PetscInt              i,n = gen_from->n,to_first = gen_to->first,to_step = gen_to->step;
    578:   PetscErrorCode        ierr;
    579:   PetscInt              from_first = gen_from->first,from_step = gen_from->step;
    580:   PetscScalar           *xv,*yv;
    581:
    583:   VecGetArray(x,&xv);
    584:   if (x != y) {VecGetArray(y,&yv);} else {yv = xv;}

    586:   if (mode & SCATTER_REVERSE){
    587:     from_first = gen_to->first;
    588:     to_first   = gen_from->first;
    589:     from_step  = gen_to->step;
    590:     to_step    = gen_from->step;
    591:   }

    593:   if (addv == INSERT_VALUES) {
    594:     if (to_step == 1 && from_step == 1) {
    595:       PetscMemcpy(yv+to_first,xv+from_first,n*sizeof(PetscScalar));
    596:     } else  {
    597:       for (i=0; i<n; i++) {
    598:         yv[to_first + i*to_step] = xv[from_first+i*from_step];
    599:       }
    600:     }
    601:   } else if (addv == ADD_VALUES) {
    602:     if (to_step == 1 && from_step == 1) {
    603:       yv += to_first; xv += from_first;
    604:       for (i=0; i<n; i++) {
    605:         yv[i] += xv[i];
    606:       }
    607:     } else {
    608:       for (i=0; i<n; i++) {
    609:         yv[to_first + i*to_step] += xv[from_first+i*from_step];
    610:       }
    611:     }
    612: #if !defined(PETSC_USE_COMPLEX)
    613:   } else if (addv == MAX_VALUES) {
    614:     if (to_step == 1 && from_step == 1) {
    615:       yv += to_first; xv += from_first;
    616:       for (i=0; i<n; i++) {
    617:         yv[i] = PetscMax(yv[i],xv[i]);
    618:       }
    619:     } else {
    620:       for (i=0; i<n; i++) {
    621:         yv[to_first + i*to_step] = PetscMax(yv[to_first + i*to_step],xv[from_first+i*from_step]);
    622:       }
    623:     }
    624: #endif
    625:   } else {SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Wrong insert option");}
    626:   VecRestoreArray(x,&xv);
    627:   if (x != y) {VecRestoreArray(y,&yv);}
    628:   return(0);
    629: }

    631: /* --------------------------------------------------------------------------*/

    636: PetscErrorCode VecScatterCopy_SGToSG(VecScatter in,VecScatter out)
    637: {
    638:   PetscErrorCode         ierr;
    639:   VecScatter_Seq_General *in_to   = (VecScatter_Seq_General*)in->todata,*out_to;
    640:   VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from;
    641:
    643:   out->postrecvs     = 0;
    644:   out->begin         = in->begin;
    645:   out->end           = in->end;
    646:   out->copy          = in->copy;
    647:   out->destroy       = in->destroy;
    648:   out->view          = in->view;

    650:   PetscMalloc4(1,VecScatter_Seq_General,&out_to,in_to->n,PetscInt,&out_to->vslots,1,VecScatter_Seq_General,&out_from,in_from->n,PetscInt,&out_from->vslots);
    651:   out_to->n                      = in_to->n;
    652:   out_to->type                   = in_to->type;
    653:   out_to->nonmatching_computed   = PETSC_FALSE;
    654:   out_to->slots_nonmatching      = 0;
    655:   out_to->is_copy                = PETSC_FALSE;
    656:   PetscMemcpy(out_to->vslots,in_to->vslots,(out_to->n)*sizeof(PetscInt));

    658:   out_from->n                    = in_from->n;
    659:   out_from->type                 = in_from->type;
    660:   out_from->nonmatching_computed = PETSC_FALSE;
    661:   out_from->slots_nonmatching    = 0;
    662:   out_from->is_copy              = PETSC_FALSE;
    663:   PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));

    665:   out->todata     = (void*)out_to;
    666:   out->fromdata   = (void*)out_from;
    667:   return(0);
    668: }

    673: PetscErrorCode VecScatterCopy_SGToSS(VecScatter in,VecScatter out)
    674: {
    675:   PetscErrorCode         ierr;
    676:   VecScatter_Seq_Stride  *in_to   = (VecScatter_Seq_Stride*)in->todata,*out_to;
    677:   VecScatter_Seq_General *in_from = (VecScatter_Seq_General*)in->fromdata,*out_from;
    678:
    680:   out->postrecvs     = 0;
    681:   out->begin         = in->begin;
    682:   out->end           = in->end;
    683:   out->copy          = in->copy;
    684:   out->destroy       = in->destroy;
    685:   out->view          = in->view;

    687:   PetscMalloc3(1,VecScatter_Seq_Stride,&out_to,1,VecScatter_Seq_General,&out_from,in_from->n,PetscInt,&out_from->vslots);
    688:   out_to->n       = in_to->n;
    689:   out_to->type    = in_to->type;
    690:   out_to->first   = in_to->first;
    691:   out_to->step    = in_to->step;
    692:   out_to->type    = in_to->type;

    694:   out_from->n                    = in_from->n;
    695:   out_from->type                 = in_from->type;
    696:   out_from->nonmatching_computed = PETSC_FALSE;
    697:   out_from->slots_nonmatching    = 0;
    698:   out_from->is_copy              = PETSC_FALSE;
    699:   PetscMemcpy(out_from->vslots,in_from->vslots,(out_from->n)*sizeof(PetscInt));

    701:   out->todata     = (void*)out_to;
    702:   out->fromdata   = (void*)out_from;
    703:   return(0);
    704: }

    706: /* --------------------------------------------------------------------------*/
    707: /*
    708:     Scatter: parallel to sequential vector, sequential strides for both.
    709: */
    712: PetscErrorCode VecScatterCopy_SStoSS(VecScatter in,VecScatter out)
    713: {
    714:   VecScatter_Seq_Stride *in_to   = (VecScatter_Seq_Stride*)in->todata,*out_to;
    715:   VecScatter_Seq_Stride *in_from = (VecScatter_Seq_Stride*)in->fromdata,*out_from;
    716:   PetscErrorCode        ierr;

    719:   out->postrecvs  = 0;
    720:   out->begin      = in->begin;
    721:   out->end        = in->end;
    722:   out->copy       = in->copy;
    723:   out->destroy    = in->destroy;
    724:   out->view       = in->view;

    726:   PetscMalloc2(1,VecScatter_Seq_Stride,&out_to,1,VecScatter_Seq_Stride,&out_from);
    727:   out_to->n       = in_to->n;
    728:   out_to->type    = in_to->type;
    729:   out_to->first   = in_to->first;
    730:   out_to->step    = in_to->step;
    731:   out_to->type    = in_to->type;
    732:   out_from->n     = in_from->n;
    733:   out_from->type  = in_from->type;
    734:   out_from->first = in_from->first;
    735:   out_from->step  = in_from->step;
    736:   out_from->type  = in_from->type;
    737:   out->todata     = (void*)out_to;
    738:   out->fromdata   = (void*)out_from;
    739:   return(0);
    740: }

    742: EXTERN PetscErrorCode VecScatterCreate_PtoS(PetscInt,PetscInt *,PetscInt,PetscInt *,Vec,Vec,PetscInt,VecScatter);
    743: EXTERN PetscErrorCode VecScatterCreate_PtoP(PetscInt,PetscInt *,PetscInt,PetscInt *,Vec,Vec,VecScatter);
    744: EXTERN PetscErrorCode VecScatterCreate_StoP(PetscInt,PetscInt *,PetscInt,PetscInt *,Vec,PetscInt,VecScatter);

    746: /* =======================================================================*/
    747: #define VEC_SEQ_ID 0
    748: #define VEC_MPI_ID 1

    750: /*
    751:    Blocksizes we have optimized scatters for
    752: */

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

    758: /*@C
    759:    VecScatterCreate - Creates a vector scatter context.

    761:    Collective on Vec

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

    771:    Output Parameter:
    772: .  newctx - location to store the new scatter context

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

    783:     Level: intermediate

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

    794:    Concepts: scatter^between vectors
    795:    Concepts: gather^between vectors

    797: .seealso: VecScatterDestroy(), VecScatterCreateToAll(), VecScatterCreateToZero()
    798: @*/
    799: PetscErrorCode PETSCVEC_DLLEXPORT VecScatterCreate(Vec xin,IS ix,Vec yin,IS iy,VecScatter *newctx)
    800: {
    801:   VecScatter     ctx;
    803:   PetscMPIInt    size;
    804:   PetscInt       totalv,*range,xin_type = VEC_SEQ_ID,yin_type = VEC_SEQ_ID;
    805:   MPI_Comm       comm,ycomm;
    806:   PetscTruth     ixblock,iyblock,iystride,islocal,cando,flag;
    807:   IS             tix = 0,tiy = 0;

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

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

    830:   ctx->beginandendtogether = PETSC_FALSE;
    831:   PetscOptionsHasName(PETSC_NULL,"-vecscatter_merge",&ctx->beginandendtogether);
    832:   if (ctx->beginandendtogether) {
    833:     PetscLogInfo((ctx,"VecScatterCreate:Using combined (merged) vector scatter begin and end\n"));
    834:   }
    835:   PetscOptionsHasName(PETSC_NULL,"-vecscatter_packtogether",&ctx->packtogether);
    836:   if (ctx->packtogether) {
    837:     PetscLogInfo((ctx,"VecScatterCreate:Pack all messages before sending\n"));
    838:   }

    840:   VecGetLocalSize(xin,&ctx->from_n);
    841:   VecGetLocalSize(yin,&ctx->to_n);

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

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

    884:   /* ===========================================================================================================
    885:         Check for special cases
    886:      ==========================================================================================================*/
    887:   /* ---------------------------------------------------------------------------*/
    888:   if (xin_type == VEC_SEQ_ID && yin_type == VEC_SEQ_ID) {
    889:     if (ix->type == IS_GENERAL && iy->type == IS_GENERAL){
    890:       PetscInt               nx,ny,*idx,*idy;
    891:       VecScatter_Seq_General *to,*from;

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

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

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

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

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

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

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

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

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

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

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

    1135:       if (cando) {
    1136:         PetscMap map;

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

    1163:     /* test for special case of processor 0 getting entire vector */
    1164:     totalv = 0;
    1165:     if (ix->type == IS_STRIDE && iy->type == IS_STRIDE){
    1166:       PetscInt             i,nx,ny,to_first,to_step,from_first,from_step,N;
    1167:       PetscMPIInt          rank,*count;
    1168:       VecScatter_MPI_ToAll *sto;

    1170:       PetscObjectGetComm((PetscObject)xin,&comm);
    1171:       MPI_Comm_rank(comm,&rank);
    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_ERR_ARG_SIZ,"Local scatter sizes don't match");
    1177:       if (!rank) {
    1178:         VecGetSize(xin,&N);
    1179:         if (nx != N) {
    1180:           totalv = 0;
    1181:         } else if (from_first == 0        && from_step == 1 &&
    1182:                    from_first == to_first && from_step == to_step){
    1183:           totalv = 1;
    1184:         } else totalv = 0;
    1185:       } else {
    1186:         if (!nx) totalv = 1;
    1187:         else     totalv = 0;
    1188:       }
    1189:       MPI_Allreduce(&totalv,&cando,1,MPI_INT,MPI_LAND,xin->comm);

    1191:       if (cando) {
    1192:         PetscMap map;

    1194:         MPI_Comm_size(ctx->comm,&size);
    1195:         PetscMalloc2(1,VecScatter_MPI_ToAll,&sto,size,PetscMPIInt,&count);
    1196:         VecGetPetscMap(xin,&map);
    1197:         PetscMapGetGlobalRange(map,&range);
    1198:         for (i=0; i<size; i++) {
    1199:           count[i] = range[i+1] - range[i];
    1200:         }
    1201:         sto->count        = count;
    1202:         sto->work1        = 0;
    1203:         sto->work2        = 0;
    1204:         sto->type         = VEC_SCATTER_MPI_TOONE;
    1205:         ctx->todata       = (void*)sto;
    1206:         ctx->fromdata     = 0;
    1207:         ctx->postrecvs    = 0;
    1208:         ctx->begin        = VecScatterBegin_MPI_ToOne;
    1209:         ctx->end          = 0;
    1210:         ctx->destroy      = VecScatterDestroy_MPI_ToAll;
    1211:         ctx->copy         = VecScatterCopy_MPI_ToAll;
    1212:         PetscLogInfo((xin,"VecScatterCreate:Special case: processor zero gets entire parallel vector, rest get none\n"));
    1213:         goto functionend;
    1214:       }
    1215:     } else {
    1216:       MPI_Allreduce(&totalv,&cando,1,MPI_INT,MPI_LAND,xin->comm);
    1217:     }

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

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

    1344:     /* general case */
    1345:     {
    1346:       PetscInt nx,ny,*idx,*idy;
    1347:       ISGetLocalSize(ix,&nx);
    1348:       ISGetIndices(ix,&idx);
    1349:       ISGetLocalSize(iy,&ny);
    1350:       ISGetIndices(iy,&idy);
    1351:       if (nx != ny) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
    1352:       VecScatterCreate_StoP(nx,idx,ny,idy,yin,1,ctx);
    1353:       ISRestoreIndices(ix,&idx);
    1354:       ISRestoreIndices(iy,&idy);
    1355:       PetscLogInfo((xin,"VecScatterCreate:General case: Seq to MPI\n"));
    1356:       goto functionend;
    1357:     }
    1358:   }
    1359:   /* ---------------------------------------------------------------------------*/
    1360:   if (xin_type == VEC_MPI_ID && yin_type == VEC_MPI_ID) {
    1361:     /* no special cases for now */
    1362:     PetscInt nx,ny,*idx,*idy;
    1363:     ISGetLocalSize(ix,&nx);
    1364:     ISGetIndices(ix,&idx);
    1365:     ISGetLocalSize(iy,&ny);
    1366:     ISGetIndices(iy,&idy);
    1367:     if (nx != ny) SETERRQ(PETSC_ERR_ARG_SIZ,"Local scatter sizes don't match");
    1368:     VecScatterCreate_PtoP(nx,idx,ny,idy,xin,yin,ctx);
    1369:     ISRestoreIndices(ix,&idx);
    1370:     ISRestoreIndices(iy,&idy);
    1371:     PetscLogInfo((xin,"VecScatterCreate:General case: MPI to MPI\n"));
    1372:     goto functionend;
    1373:   }

    1375:   functionend:
    1376:   *newctx = ctx;
    1377:   if (tix) {ISDestroy(tix);}
    1378:   if (tiy) {ISDestroy(tiy);}
    1379:   PetscOptionsHasName(PETSC_NULL,"-vecscatter_view_info",&flag);
    1380:   if (flag) {
    1381:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(comm),PETSC_VIEWER_ASCII_INFO);
    1382:     VecScatterView(ctx,PETSC_VIEWER_STDOUT_(comm));
    1383:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(comm));
    1384:   }
    1385:   PetscOptionsHasName(PETSC_NULL,"-vecscatter_view",&flag);
    1386:   if (flag) {
    1387:     VecScatterView(ctx,PETSC_VIEWER_STDOUT_(comm));
    1388:   }
    1389:   return(0);
    1390: }

    1392: /* ------------------------------------------------------------------*/
    1395: /*@
    1396:    VecScatterPostRecvs - Posts the receives required for the ready-receiver
    1397:    version of the VecScatter routines.

    1399:    Collective on VecScatter

    1401:    Input Parameters:
    1402: +  x - the vector from which we scatter (not needed,can be null)
    1403: .  y - the vector to which we scatter
    1404: .  addv - either ADD_VALUES or INSERT_VALUES
    1405: .  mode - the scattering mode, usually SCATTER_FORWARD.  The available modes are:
    1406:     SCATTER_FORWARD, SCATTER_REVERSE
    1407: -  inctx - scatter context generated by VecScatterCreate()

    1409:    Output Parameter:
    1410: .  y - the vector to which we scatter

    1412:    Level: advanced

    1414:    Notes:
    1415:    If you use SCATTER_REVERSE the first two arguments should be reversed, from
    1416:    the SCATTER_FORWARD.
    1417:    The vectors x and y cannot be the same. y[iy[i]] = x[ix[i]], for i=0,...,ni-1

    1419:    This scatter is far more general than the conventional
    1420:    scatter, since it can be a gather or a scatter or a combination,
    1421:    depending on the indices ix and iy.  If x is a parallel vector and y
    1422:    is sequential, VecScatterBegin() can serve to gather values to a
    1423:    single processor.  Similarly, if y is parallel and x sequential, the
    1424:    routine can scatter from one processor to many processors.

    1426: .seealso: VecScatterCreate(), VecScatterEnd(), VecScatterBegin()
    1427: @*/
    1428: PetscErrorCode PETSCVEC_DLLEXPORT VecScatterPostRecvs(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter inctx)
    1429: {

    1436:   if (inctx->postrecvs) {
    1437:     (*inctx->postrecvs)(x,y,addv,mode,inctx);
    1438:   }
    1439:   return(0);
    1440: }

    1442: /* ------------------------------------------------------------------*/
    1445: /*@
    1446:    VecScatterGetMerged - Returns true if the scatter is completed in the VecScatterBegin()
    1447:       and the VecScatterEnd() does nothing

    1449:    Not Collective

    1451:    Input Parameter:
    1452: .   ctx - scatter context created with VecScatterCreate()

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

    1457:    Level: developer

    1459: .seealso: VecScatterCreate(), VecScatterEnd(), VecScatterBegin()
    1460: @*/
    1461: PetscErrorCode PETSCVEC_DLLEXPORT VecScatterGetMerged(VecScatter ctx,PetscTruth *flg)
    1462: {
    1465:   *flg = ctx->beginandendtogether;
    1466:   return(0);
    1467: }

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

    1475:    Collective on VecScatter and Vec

    1477:    Input Parameters:
    1478: +  x - the vector from which we scatter
    1479: .  y - the vector to which we scatter
    1480: .  addv - either ADD_VALUES or INSERT_VALUES, with INSERT_VALUES mode any location
    1481:           not scattered to retains its old value; i.e. the vector is NOT first zeroed.
    1482: .  mode - the scattering mode, usually SCATTER_FORWARD.  The available modes are:
    1483:     SCATTER_FORWARD or SCATTER_REVERSE
    1484: -  inctx - scatter context generated by VecScatterCreate()

    1486:    Level: intermediate

    1488:    Notes:
    1489:    The vectors x and y need not be the same vectors used in the call
    1490:    to VecScatterCreate(), but x must have the same parallel data layout
    1491:    as that passed in as the x to VecScatterCreate(), similarly for the y.
    1492:    Most likely they have been obtained from VecDuplicate().

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

    1497:    If you use SCATTER_REVERSE the first two arguments should be reversed, from
    1498:    the SCATTER_FORWARD.
    1499:
    1500:    y[iy[i]] = x[ix[i]], for i=0,...,ni-1

    1502:    This scatter is far more general than the conventional
    1503:    scatter, since it can be a gather or a scatter or a combination,
    1504:    depending on the indices ix and iy.  If x is a parallel vector and y
    1505:    is sequential, VecScatterBegin() can serve to gather values to a
    1506:    single processor.  Similarly, if y is parallel and x sequential, the
    1507:    routine can scatter from one processor to many processors.

    1509:    Concepts: scatter^between vectors
    1510:    Concepts: gather^between vectors

    1512: .seealso: VecScatterCreate(), VecScatterEnd()
    1513: @*/
    1514: PetscErrorCode PETSCVEC_DLLEXPORT VecScatterBegin(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter inctx)
    1515: {
    1517: #if defined(PETSC_USE_DEBUG)
    1518:   PetscInt      to_n,from_n;
    1519: #endif

    1525:   if (inctx->inuse) SETERRQ(PETSC_ERR_ARG_WRONGSTATE," Scatter ctx already in use");
    1526: #if defined(PETSC_USE_DEBUG)
    1527:   /*
    1528:      Error checking to make sure these vectors match the vectors used
    1529:    to create the vector scatter context. -1 in the from_n and to_n indicate the
    1530:    vector lengths are unknown (for example with mapped scatters) and thus
    1531:    no error checking is performed.
    1532:   */
    1533:   if (inctx->from_n >= 0 && inctx->to_n >= 0) {
    1534:     VecGetLocalSize(x,&from_n);
    1535:     VecGetLocalSize(y,&to_n);
    1536:     if (mode & SCATTER_REVERSE) {
    1537:       if (to_n != inctx->from_n) SETERRQ(PETSC_ERR_ARG_SIZ,"Vector wrong size for scatter (scatter reverse and vector to != ctx from size)");
    1538:       if (from_n != inctx->to_n) SETERRQ(PETSC_ERR_ARG_SIZ,"Vector wrong size for scatter (scatter reverse and vector from != ctx to size)");
    1539:     } else {
    1540:       if (to_n != inctx->to_n)     SETERRQ(PETSC_ERR_ARG_SIZ,"Vector wrong size for scatter (scatter forward and vector to != ctx to size)");
    1541:       if (from_n != inctx->from_n) SETERRQ(PETSC_ERR_ARG_SIZ,"Vector wrong size for scatter (scatter forward and vector from != ctx from size)");
    1542:     }
    1543:   }
    1544: #endif

    1546:   inctx->inuse = PETSC_TRUE;
    1547:   PetscLogEventBarrierBegin(VEC_ScatterBarrier,0,0,0,0,inctx->comm);
    1548:   (*inctx->begin)(x,y,addv,mode,inctx);
    1549:   if (inctx->beginandendtogether && inctx->end) {
    1550:     inctx->inuse = PETSC_FALSE;
    1551:     (*inctx->end)(x,y,addv,mode,inctx);
    1552:   }
    1553:   PetscLogEventBarrierEnd(VEC_ScatterBarrier,0,0,0,0,inctx->comm);
    1554:   return(0);
    1555: }

    1557: /* --------------------------------------------------------------------*/
    1560: /*@
    1561:    VecScatterEnd - Ends a generalized scatter from one vector to another.  Call
    1562:    after first calling VecScatterBegin().

    1564:    Collective on VecScatter and Vec

    1566:    Input Parameters:
    1567: +  x - the vector from which we scatter
    1568: .  y - the vector to which we scatter
    1569: .  addv - either ADD_VALUES or INSERT_VALUES.
    1570: .  mode - the scattering mode, usually SCATTER_FORWARD.  The available modes are:
    1571:      SCATTER_FORWARD, SCATTER_REVERSE
    1572: -  ctx - scatter context generated by VecScatterCreate()

    1574:    Level: intermediate

    1576:    Notes:
    1577:    If you use SCATTER_REVERSE the first two arguments should be reversed, from
    1578:    the SCATTER_FORWARD.
    1579:    y[iy[i]] = x[ix[i]], for i=0,...,ni-1

    1581: .seealso: VecScatterBegin(), VecScatterCreate()
    1582: @*/
    1583: PetscErrorCode PETSCVEC_DLLEXPORT VecScatterEnd(Vec x,Vec y,InsertMode addv,ScatterMode mode,VecScatter ctx)
    1584: {

    1591:   ctx->inuse = PETSC_FALSE;
    1592:   if (!ctx->end) return(0);
    1593:   if (!ctx->beginandendtogether) {
    1594:     PetscLogEventBegin(VEC_ScatterEnd,ctx,x,y,0);
    1595:     (*(ctx)->end)(x,y,addv,mode,ctx);
    1596:     PetscLogEventEnd(VEC_ScatterEnd,ctx,x,y,0);
    1597:   }
    1598:   return(0);
    1599: }

    1603: /*@C
    1604:    VecScatterDestroy - Destroys a scatter context created by
    1605:    VecScatterCreate().

    1607:    Collective on VecScatter

    1609:    Input Parameter:
    1610: .  ctx - the scatter context

    1612:    Level: intermediate

    1614: .seealso: VecScatterCreate(), VecScatterCopy()
    1615: @*/
    1616: PetscErrorCode PETSCVEC_DLLEXPORT VecScatterDestroy(VecScatter ctx)
    1617: {

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

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

    1627:   (*ctx->destroy)(ctx);
    1628:   return(0);
    1629: }

    1633: /*@C
    1634:    VecScatterCopy - Makes a copy of a scatter context.

    1636:    Collective on VecScatter

    1638:    Input Parameter:
    1639: .  sctx - the scatter context

    1641:    Output Parameter:
    1642: .  ctx - the context copy

    1644:    Level: advanced

    1646: .seealso: VecScatterCreate(), VecScatterDestroy()
    1647: @*/
    1648: PetscErrorCode PETSCVEC_DLLEXPORT VecScatterCopy(VecScatter sctx,VecScatter *ctx)
    1649: {

    1655:   if (!sctx->copy) SETERRQ(PETSC_ERR_SUP,"Cannot copy this type");
    1656:   PetscHeaderCreate(*ctx,_p_VecScatter,int,VEC_SCATTER_COOKIE,0,"VecScatter",sctx->comm,VecScatterDestroy,VecScatterView);
    1657:   (*ctx)->to_n   = sctx->to_n;
    1658:   (*ctx)->from_n = sctx->from_n;
    1659:   (*sctx->copy)(sctx,*ctx);
    1660:   return(0);
    1661: }

    1664: /* ------------------------------------------------------------------*/
    1667: /*@
    1668:    VecScatterView - Views a vector scatter context.

    1670:    Collective on VecScatter

    1672:    Input Parameters:
    1673: +  ctx - the scatter context
    1674: -  viewer - the viewer for displaying the context

    1676:    Level: intermediate

    1678: @*/
    1679: PetscErrorCode PETSCVEC_DLLEXPORT VecScatterView(VecScatter ctx,PetscViewer viewer)
    1680: {

    1685:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(ctx->comm);
    1687:   if (!ctx->view) SETERRQ(PETSC_ERR_SUP,"Cannot view this type of scatter context yet");

    1689:   (*ctx->view)(ctx,viewer);
    1690:   return(0);
    1691: }

    1695: /*@
    1696:    VecScatterRemap - Remaps the "from" and "to" indices in a
    1697:    vector scatter context. FOR EXPERTS ONLY!

    1699:    Collective on VecScatter

    1701:    Input Parameters:
    1702: +  scat - vector scatter context
    1703: .  from - remapping for "from" indices (may be PETSC_NULL)
    1704: -  to   - remapping for "to" indices (may be PETSC_NULL)

    1706:    Level: developer

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

    1712:           In the sequential case the todata is the indices where the
    1713:           data is put and the fromdata is where it is taken from.
    1714:           This is backwards from the paralllel case! CRY! CRY! CRY!

    1716: @*/
    1717: PetscErrorCode PETSCVEC_DLLEXPORT VecScatterRemap(VecScatter scat,PetscInt *rto,PetscInt *rfrom)
    1718: {
    1719:   VecScatter_Seq_General *to,*from;
    1720:   VecScatter_MPI_General *mto;
    1721:   PetscInt               i;

    1728:   from = (VecScatter_Seq_General *)scat->fromdata;
    1729:   mto  = (VecScatter_MPI_General *)scat->todata;

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

    1733:   if (rto) {
    1734:     if (mto->type == VEC_SCATTER_MPI_GENERAL) {
    1735:       /* handle off processor parts */
    1736:       for (i=0; i<mto->starts[mto->n]; i++) {
    1737:         mto->indices[i] = rto[mto->indices[i]];
    1738:       }
    1739:       /* handle local part */
    1740:       to = &mto->local;
    1741:       for (i=0; i<to->n; i++) {
    1742:         to->vslots[i] = rto[to->vslots[i]];
    1743:       }
    1744:     } else if (from->type == VEC_SCATTER_SEQ_GENERAL) {
    1745:       for (i=0; i<from->n; i++) {
    1746:         from->vslots[i] = rto[from->vslots[i]];
    1747:       }
    1748:     } else if (from->type == VEC_SCATTER_SEQ_STRIDE) {
    1749:       VecScatter_Seq_Stride *sto = (VecScatter_Seq_Stride*)from;
    1750:
    1751:       /* if the remapping is the identity and stride is identity then skip remap */
    1752:       if (sto->step == 1 && sto->first == 0) {
    1753:         for (i=0; i<sto->n; i++) {
    1754:           if (rto[i] != i) {
    1755:             SETERRQ(PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
    1756:           }
    1757:         }
    1758:       } else SETERRQ(PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
    1759:     } else SETERRQ(PETSC_ERR_ARG_SIZ,"Unable to remap such scatters");
    1760:   }

    1762:   if (rfrom) {
    1763:     SETERRQ(PETSC_ERR_SUP,"Unable to remap the FROM in scatters yet");
    1764:   }

    1766:   /*
    1767:      Mark then vector lengths as unknown because we do not know the
    1768:    lengths of the remapped vectors
    1769:   */
    1770:   scat->from_n = -1;
    1771:   scat->to_n   = -1;

    1773:   return(0);
    1774: }