Actual source code: pbvec.c

      1: /*
      2:    This file contains routines for Parallel vector operations.
      3:  */
      4: #include <petscsys.h>
      5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>

      7: PETSC_INTERN PetscErrorCode VecView_MPI_Draw(Vec, PetscViewer);

      9: PetscErrorCode VecPlaceArray_MPI(Vec vin, const PetscScalar *a)
     10: {
     11:   Vec_MPI *v = (Vec_MPI *)vin->data;

     13:   PetscFunctionBegin;
     14:   PetscCheck(!v->unplacedarray, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "VecPlaceArray() was already called on this vector, without a call to VecResetArray()");
     15:   v->unplacedarray = v->array; /* save previous array so reset can bring it back */
     16:   v->array         = (PetscScalar *)a;
     17:   if (v->localrep) PetscCall(VecPlaceArray(v->localrep, a));
     18:   PetscFunctionReturn(PETSC_SUCCESS);
     19: }

     21: PetscErrorCode VecDuplicate_MPI(Vec win, Vec *v)
     22: {
     23:   Vec_MPI     *vw, *w = (Vec_MPI *)win->data;
     24:   PetscScalar *array;

     26:   PetscFunctionBegin;
     27:   PetscCall(VecCreateWithLayout_Private(win->map, v));

     29:   PetscCall(VecCreate_MPI_Private(*v, PETSC_TRUE, w->nghost, NULL));
     30:   vw           = (Vec_MPI *)(*v)->data;
     31:   (*v)->ops[0] = win->ops[0];

     33:   /* save local representation of the parallel vector (and scatter) if it exists */
     34:   if (w->localrep) {
     35:     PetscCall(VecGetArray(*v, &array));
     36:     PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, PetscAbs(win->map->bs), win->map->n + w->nghost, array, &vw->localrep));
     37:     vw->localrep->ops[0] = w->localrep->ops[0];
     38:     PetscCall(VecRestoreArray(*v, &array));

     40:     vw->localupdate = w->localupdate;
     41:     if (vw->localupdate) PetscCall(PetscObjectReference((PetscObject)vw->localupdate));

     43:     vw->ghost = w->ghost;
     44:     if (vw->ghost) PetscCall(PetscObjectReference((PetscObject)vw->ghost));
     45:   }

     47:   /* New vector should inherit stashing property of parent */
     48:   (*v)->stash.donotstash   = win->stash.donotstash;
     49:   (*v)->stash.ignorenegidx = win->stash.ignorenegidx;

     51:   PetscCall(PetscObjectListDuplicate(((PetscObject)win)->olist, &((PetscObject)*v)->olist));
     52:   PetscCall(PetscFunctionListDuplicate(((PetscObject)win)->qlist, &((PetscObject)*v)->qlist));

     54:   (*v)->stash.bs  = win->stash.bs;
     55:   (*v)->bstash.bs = win->bstash.bs;
     56:   PetscFunctionReturn(PETSC_SUCCESS);
     57: }

     59: static PetscErrorCode VecDuplicateVecs_MPI_GEMV(Vec w, PetscInt m, Vec *V[])
     60: {
     61:   Vec_MPI *wmpi = (Vec_MPI *)w->data;

     63:   PetscFunctionBegin;
     64:   // Currently only do GEMV for vectors without ghosts. Note w might be a VECMPI subclass object.
     65:   // This routine relies on the duplicate operation being VecDuplicate_MPI. If not, bail out to the default.
     66:   if (wmpi->localrep || w->ops->duplicate != VecDuplicate_MPI) {
     67:     w->ops->duplicatevecs = VecDuplicateVecs_Default;
     68:     PetscCall(VecDuplicateVecs(w, m, V));
     69:   } else {
     70:     PetscScalar *array;
     71:     PetscInt64   lda; // use 64-bit as we will do "m * lda"

     73:     PetscCall(PetscMalloc1(m, V));
     74:     VecGetLocalSizeAligned(w, 64, &lda); // get in lda the 64-bytes aligned local size

     76:     PetscCall(PetscCalloc1(m * lda, &array));
     77:     for (PetscInt i = 0; i < m; i++) {
     78:       Vec v;
     79:       PetscCall(VecCreateMPIWithLayoutAndArray_Private(w->map, PetscSafePointerPlusOffset(array, i * lda), &v));
     80:       PetscCall(VecSetType(v, ((PetscObject)w)->type_name));
     81:       PetscCall(PetscObjectListDuplicate(((PetscObject)w)->olist, &((PetscObject)v)->olist));
     82:       PetscCall(PetscFunctionListDuplicate(((PetscObject)w)->qlist, &((PetscObject)v)->qlist));
     83:       v->ops[0]             = w->ops[0];
     84:       v->stash.donotstash   = w->stash.donotstash;
     85:       v->stash.ignorenegidx = w->stash.ignorenegidx;
     86:       v->stash.bs           = w->stash.bs;
     87:       v->bstash.bs          = w->bstash.bs;
     88:       (*V)[i]               = v;
     89:     }
     90:     // So when the first vector is destroyed it will destroy the array
     91:     if (m) ((Vec_MPI *)(*V)[0]->data)->array_allocated = array;
     92:     // disable replacearray of the first vector, as freeing its memory also frees others in the group.
     93:     // But replacearray of others is ok, as they don't own their array.
     94:     if (m > 1) (*V)[0]->ops->replacearray = VecReplaceArray_Default_GEMV_Error;
     95:   }
     96:   PetscFunctionReturn(PETSC_SUCCESS);
     97: }

     99: static PetscErrorCode VecSetOption_MPI(Vec V, VecOption op, PetscBool flag)
    100: {
    101:   Vec_MPI *v = (Vec_MPI *)V->data;

    103:   PetscFunctionBegin;
    104:   switch (op) {
    105:   case VEC_IGNORE_OFF_PROC_ENTRIES:
    106:     V->stash.donotstash = flag;
    107:     break;
    108:   case VEC_IGNORE_NEGATIVE_INDICES:
    109:     V->stash.ignorenegidx = flag;
    110:     break;
    111:   case VEC_SUBSET_OFF_PROC_ENTRIES:
    112:     v->assembly_subset = flag;              /* See the same logic in MatAssembly wrt MAT_SUBSET_OFF_PROC_ENTRIES */
    113:     if (!v->assembly_subset) {              /* User indicates "do not reuse the communication pattern" */
    114:       PetscCall(VecAssemblyReset_MPI(V));   /* Reset existing pattern to free memory */
    115:       v->first_assembly_done = PETSC_FALSE; /* Mark the first assembly is not done */
    116:     }
    117:     break;
    118:   }
    119:   PetscFunctionReturn(PETSC_SUCCESS);
    120: }

    122: PetscErrorCode VecResetArray_MPI(Vec vin)
    123: {
    124:   Vec_MPI *v = (Vec_MPI *)vin->data;

    126:   PetscFunctionBegin;
    127:   v->array         = v->unplacedarray;
    128:   v->unplacedarray = NULL;
    129:   if (v->localrep) PetscCall(VecResetArray(v->localrep));
    130:   PetscFunctionReturn(PETSC_SUCCESS);
    131: }

    133: static PetscErrorCode VecAssemblySend_MPI_Private(MPI_Comm comm, const PetscMPIInt tag[], PetscMPIInt rankid, PetscMPIInt rank, void *sdata, MPI_Request req[], void *ctx)
    134: {
    135:   Vec                X   = (Vec)ctx;
    136:   Vec_MPI           *x   = (Vec_MPI *)X->data;
    137:   VecAssemblyHeader *hdr = (VecAssemblyHeader *)sdata;
    138:   PetscInt           bs  = X->map->bs;

    140:   PetscFunctionBegin;
    141:   /* x->first_assembly_done indicates we are reusing a communication network. In that case, some
    142:      messages can be empty, but we have to send them this time if we sent them before because the
    143:      receiver is expecting them.
    144:    */
    145:   if (hdr->count || (x->first_assembly_done && x->sendptrs[rankid].ints)) {
    146:     PetscCallMPI(MPIU_Isend(x->sendptrs[rankid].ints, hdr->count, MPIU_INT, rank, tag[0], comm, &req[0]));
    147:     PetscCallMPI(MPIU_Isend(x->sendptrs[rankid].scalars, hdr->count, MPIU_SCALAR, rank, tag[1], comm, &req[1]));
    148:   }
    149:   if (hdr->bcount || (x->first_assembly_done && x->sendptrs[rankid].intb)) {
    150:     PetscCallMPI(MPIU_Isend(x->sendptrs[rankid].intb, hdr->bcount, MPIU_INT, rank, tag[2], comm, &req[2]));
    151:     PetscCallMPI(MPIU_Isend(x->sendptrs[rankid].scalarb, hdr->bcount * bs, MPIU_SCALAR, rank, tag[3], comm, &req[3]));
    152:   }
    153:   PetscFunctionReturn(PETSC_SUCCESS);
    154: }

    156: static PetscErrorCode VecAssemblyRecv_MPI_Private(MPI_Comm comm, const PetscMPIInt tag[], PetscMPIInt rank, void *rdata, MPI_Request req[], void *ctx)
    157: {
    158:   Vec                X   = (Vec)ctx;
    159:   Vec_MPI           *x   = (Vec_MPI *)X->data;
    160:   VecAssemblyHeader *hdr = (VecAssemblyHeader *)rdata;
    161:   PetscInt           bs  = X->map->bs;
    162:   VecAssemblyFrame  *frame;

    164:   PetscFunctionBegin;
    165:   PetscCall(PetscSegBufferGet(x->segrecvframe, 1, &frame));

    167:   if (hdr->count) {
    168:     PetscCall(PetscSegBufferGet(x->segrecvint, hdr->count, &frame->ints));
    169:     PetscCallMPI(MPIU_Irecv(frame->ints, hdr->count, MPIU_INT, rank, tag[0], comm, &req[0]));
    170:     PetscCall(PetscSegBufferGet(x->segrecvscalar, hdr->count, &frame->scalars));
    171:     PetscCallMPI(MPIU_Irecv(frame->scalars, hdr->count, MPIU_SCALAR, rank, tag[1], comm, &req[1]));
    172:     frame->pendings = 2;
    173:   } else {
    174:     frame->ints     = NULL;
    175:     frame->scalars  = NULL;
    176:     frame->pendings = 0;
    177:   }

    179:   if (hdr->bcount) {
    180:     PetscCall(PetscSegBufferGet(x->segrecvint, hdr->bcount, &frame->intb));
    181:     PetscCallMPI(MPIU_Irecv(frame->intb, hdr->bcount, MPIU_INT, rank, tag[2], comm, &req[2]));
    182:     PetscCall(PetscSegBufferGet(x->segrecvscalar, hdr->bcount * bs, &frame->scalarb));
    183:     PetscCallMPI(MPIU_Irecv(frame->scalarb, hdr->bcount * bs, MPIU_SCALAR, rank, tag[3], comm, &req[3]));
    184:     frame->pendingb = 2;
    185:   } else {
    186:     frame->intb     = NULL;
    187:     frame->scalarb  = NULL;
    188:     frame->pendingb = 0;
    189:   }
    190:   PetscFunctionReturn(PETSC_SUCCESS);
    191: }

    193: static PetscErrorCode VecAssemblyBegin_MPI_BTS(Vec X)
    194: {
    195:   Vec_MPI    *x = (Vec_MPI *)X->data;
    196:   MPI_Comm    comm;
    197:   PetscMPIInt i;
    198:   PetscInt    j, jb, bs;

    200:   PetscFunctionBegin;
    201:   if (X->stash.donotstash) PetscFunctionReturn(PETSC_SUCCESS);

    203:   PetscCall(PetscObjectGetComm((PetscObject)X, &comm));
    204:   PetscCall(VecGetBlockSize(X, &bs));
    205:   if (PetscDefined(USE_DEBUG)) {
    206:     InsertMode addv;

    208:     PetscCallMPI(MPIU_Allreduce((PetscEnum *)&X->stash.insertmode, (PetscEnum *)&addv, 1, MPIU_ENUM, MPI_BOR, comm));
    209:     PetscCheck(addv != (ADD_VALUES | INSERT_VALUES), comm, PETSC_ERR_ARG_NOTSAMETYPE, "Some processors inserted values while others added");
    210:   }
    211:   X->bstash.insertmode = X->stash.insertmode; /* Block stash implicitly tracks InsertMode of scalar stash */

    213:   PetscCall(VecStashSortCompress_Private(&X->stash));
    214:   PetscCall(VecStashSortCompress_Private(&X->bstash));

    216:   if (!x->sendranks) {
    217:     PetscMPIInt nowners, bnowners, *owners, *bowners;
    218:     PetscInt    ntmp;

    220:     PetscCall(VecStashGetOwnerList_Private(&X->stash, X->map, &nowners, &owners));
    221:     PetscCall(VecStashGetOwnerList_Private(&X->bstash, X->map, &bnowners, &bowners));
    222:     PetscCall(PetscMergeMPIIntArray(nowners, owners, bnowners, bowners, &ntmp, &x->sendranks));
    223:     PetscCall(PetscMPIIntCast(ntmp, &x->nsendranks));
    224:     PetscCall(PetscFree(owners));
    225:     PetscCall(PetscFree(bowners));
    226:     PetscCall(PetscMalloc1(x->nsendranks, &x->sendhdr));
    227:     PetscCall(PetscCalloc1(x->nsendranks, &x->sendptrs));
    228:   }
    229:   for (i = 0, j = 0, jb = 0; i < x->nsendranks; i++) {
    230:     PetscMPIInt rank = x->sendranks[i];

    232:     x->sendhdr[i].insertmode = X->stash.insertmode;
    233:     /* Initialize pointers for non-empty stashes the first time around.  Subsequent assemblies with
    234:      * VEC_SUBSET_OFF_PROC_ENTRIES will leave the old pointers (dangling because the stash has been collected) when
    235:      * there is nothing new to send, so that size-zero messages get sent instead. */
    236:     x->sendhdr[i].count = 0;
    237:     if (X->stash.n) {
    238:       x->sendptrs[i].ints    = &X->stash.idx[j];
    239:       x->sendptrs[i].scalars = &X->stash.array[j];
    240:       for (; j < X->stash.n && X->stash.idx[j] < X->map->range[rank + 1]; j++) x->sendhdr[i].count++;
    241:     }
    242:     x->sendhdr[i].bcount = 0;
    243:     if (X->bstash.n) {
    244:       x->sendptrs[i].intb    = &X->bstash.idx[jb];
    245:       x->sendptrs[i].scalarb = &X->bstash.array[jb * bs];
    246:       for (; jb < X->bstash.n && X->bstash.idx[jb] * bs < X->map->range[rank + 1]; jb++) x->sendhdr[i].bcount++;
    247:     }
    248:   }

    250:   if (!x->segrecvint) PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 1000, &x->segrecvint));
    251:   if (!x->segrecvscalar) PetscCall(PetscSegBufferCreate(sizeof(PetscScalar), 1000, &x->segrecvscalar));
    252:   if (!x->segrecvframe) PetscCall(PetscSegBufferCreate(sizeof(VecAssemblyFrame), 50, &x->segrecvframe));
    253:   if (x->first_assembly_done) { /* this is not the first assembly */
    254:     PetscMPIInt tag[4];
    255:     for (i = 0; i < 4; i++) PetscCall(PetscCommGetNewTag(comm, &tag[i]));
    256:     for (i = 0; i < x->nsendranks; i++) PetscCall(VecAssemblySend_MPI_Private(comm, tag, i, x->sendranks[i], x->sendhdr + i, x->sendreqs + 4 * i, X));
    257:     for (i = 0; i < x->nrecvranks; i++) PetscCall(VecAssemblyRecv_MPI_Private(comm, tag, x->recvranks[i], x->recvhdr + i, x->recvreqs + 4 * i, X));
    258:     x->use_status = PETSC_TRUE;
    259:   } else { /* First time assembly */
    260:     PetscCall(PetscCommBuildTwoSidedFReq(comm, 3, MPIU_INT, x->nsendranks, x->sendranks, (PetscInt *)x->sendhdr, &x->nrecvranks, &x->recvranks, &x->recvhdr, 4, &x->sendreqs, &x->recvreqs, VecAssemblySend_MPI_Private, VecAssemblyRecv_MPI_Private, X));
    261:     x->use_status = PETSC_FALSE;
    262:   }

    264:   /* The first_assembly_done flag is only meaningful when x->assembly_subset is set.
    265:      This line says when assembly_subset is set, then we mark that the first assembly is done.
    266:    */
    267:   x->first_assembly_done = x->assembly_subset;

    269:   {
    270:     PetscInt nstash, reallocs;
    271:     PetscCall(VecStashGetInfo_Private(&X->stash, &nstash, &reallocs));
    272:     PetscCall(PetscInfo(X, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
    273:     PetscCall(VecStashGetInfo_Private(&X->bstash, &nstash, &reallocs));
    274:     PetscCall(PetscInfo(X, "Block-Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
    275:   }
    276:   PetscFunctionReturn(PETSC_SUCCESS);
    277: }

    279: static PetscErrorCode VecAssemblyEnd_MPI_BTS(Vec X)
    280: {
    281:   Vec_MPI          *x  = (Vec_MPI *)X->data;
    282:   PetscInt          bs = X->map->bs;
    283:   PetscMPIInt       npending, *some_indices, r;
    284:   MPI_Status       *some_statuses;
    285:   PetscScalar      *xarray;
    286:   VecAssemblyFrame *frame;

    288:   PetscFunctionBegin;
    289:   if (X->stash.donotstash) {
    290:     X->stash.insertmode  = NOT_SET_VALUES;
    291:     X->bstash.insertmode = NOT_SET_VALUES;
    292:     PetscFunctionReturn(PETSC_SUCCESS);
    293:   }

    295:   PetscCheck(x->segrecvframe, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing segrecvframe! Probably you forgot to call VecAssemblyBegin() first");
    296:   PetscCall(VecGetArray(X, &xarray));
    297:   PetscCall(PetscSegBufferExtractInPlace(x->segrecvframe, &frame));
    298:   PetscCall(PetscMalloc2(4 * x->nrecvranks, &some_indices, (x->use_status ? (size_t)4 * x->nrecvranks : 0), &some_statuses));
    299:   for (r = 0, npending = 0; r < x->nrecvranks; r++) npending += frame[r].pendings + frame[r].pendingb;
    300:   while (npending > 0) {
    301:     PetscMPIInt ndone = 0, ii;
    302:     /* Filling MPI_Status fields requires some resources from the MPI library.  We skip it on the first assembly, or
    303:      * when VEC_SUBSET_OFF_PROC_ENTRIES has not been set, because we could exchange exact sizes in the initial
    304:      * rendezvous.  When the rendezvous is elided, however, we use MPI_Status to get actual message lengths, so that
    305:      * subsequent assembly can set a proper subset of the values. */
    306:     PetscCallMPI(MPI_Waitsome(4 * x->nrecvranks, x->recvreqs, &ndone, some_indices, x->use_status ? some_statuses : MPI_STATUSES_IGNORE));
    307:     for (ii = 0; ii < ndone; ii++) {
    308:       PetscInt     i     = some_indices[ii] / 4, j, k;
    309:       InsertMode   imode = (InsertMode)x->recvhdr[i].insertmode;
    310:       PetscInt    *recvint;
    311:       PetscScalar *recvscalar;
    312:       PetscBool    intmsg   = (PetscBool)(some_indices[ii] % 2 == 0);
    313:       PetscBool    blockmsg = (PetscBool)((some_indices[ii] % 4) / 2 == 1);

    315:       npending--;
    316:       if (!blockmsg) { /* Scalar stash */
    317:         PetscMPIInt count;

    319:         if (--frame[i].pendings > 0) continue;
    320:         if (x->use_status) {
    321:           PetscCallMPI(MPI_Get_count(&some_statuses[ii], intmsg ? MPIU_INT : MPIU_SCALAR, &count));
    322:         } else {
    323:           PetscCall(PetscMPIIntCast(x->recvhdr[i].count, &count));
    324:         }
    325:         for (j = 0, recvint = frame[i].ints, recvscalar = frame[i].scalars; j < count; j++, recvint++) {
    326:           PetscInt loc = *recvint - X->map->rstart;

    328:           PetscCheck(*recvint >= X->map->rstart && X->map->rend > *recvint, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Received vector entry %" PetscInt_FMT " out of local range [%" PetscInt_FMT ",%" PetscInt_FMT ")]", *recvint, X->map->rstart, X->map->rend);
    329:           switch (imode) {
    330:           case ADD_VALUES:
    331:             xarray[loc] += *recvscalar++;
    332:             break;
    333:           case INSERT_VALUES:
    334:             xarray[loc] = *recvscalar++;
    335:             break;
    336:           default:
    337:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Insert mode not supported 0x%x", imode);
    338:           }
    339:         }
    340:       } else { /* Block stash */
    341:         PetscMPIInt count;
    342:         if (--frame[i].pendingb > 0) continue;
    343:         if (x->use_status) {
    344:           PetscCallMPI(MPI_Get_count(&some_statuses[ii], intmsg ? MPIU_INT : MPIU_SCALAR, &count));
    345:           if (!intmsg) count /= bs; /* Convert from number of scalars to number of blocks */
    346:         } else {
    347:           PetscCall(PetscMPIIntCast(x->recvhdr[i].bcount, &count));
    348:         }
    349:         for (j = 0, recvint = frame[i].intb, recvscalar = frame[i].scalarb; j < count; j++, recvint++) {
    350:           PetscInt loc = (*recvint) * bs - X->map->rstart;
    351:           switch (imode) {
    352:           case ADD_VALUES:
    353:             for (k = loc; k < loc + bs; k++) xarray[k] += *recvscalar++;
    354:             break;
    355:           case INSERT_VALUES:
    356:             for (k = loc; k < loc + bs; k++) xarray[k] = *recvscalar++;
    357:             break;
    358:           default:
    359:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Insert mode not supported 0x%x", imode);
    360:           }
    361:         }
    362:       }
    363:     }
    364:   }
    365:   PetscCall(VecRestoreArray(X, &xarray));
    366:   PetscCallMPI(MPI_Waitall(4 * x->nsendranks, x->sendreqs, MPI_STATUSES_IGNORE));
    367:   PetscCall(PetscFree2(some_indices, some_statuses));
    368:   if (x->assembly_subset) {
    369:     PetscCall(PetscSegBufferExtractInPlace(x->segrecvint, NULL));
    370:     PetscCall(PetscSegBufferExtractInPlace(x->segrecvscalar, NULL));
    371:   } else {
    372:     PetscCall(VecAssemblyReset_MPI(X));
    373:   }

    375:   X->stash.insertmode  = NOT_SET_VALUES;
    376:   X->bstash.insertmode = NOT_SET_VALUES;
    377:   PetscCall(VecStashScatterEnd_Private(&X->stash));
    378:   PetscCall(VecStashScatterEnd_Private(&X->bstash));
    379:   PetscFunctionReturn(PETSC_SUCCESS);
    380: }

    382: PetscErrorCode VecAssemblyReset_MPI(Vec X)
    383: {
    384:   Vec_MPI *x = (Vec_MPI *)X->data;

    386:   PetscFunctionBegin;
    387:   PetscCall(PetscFree(x->sendreqs));
    388:   PetscCall(PetscFree(x->recvreqs));
    389:   PetscCall(PetscFree(x->sendranks));
    390:   PetscCall(PetscFree(x->recvranks));
    391:   PetscCall(PetscFree(x->sendhdr));
    392:   PetscCall(PetscFree(x->recvhdr));
    393:   PetscCall(PetscFree(x->sendptrs));
    394:   PetscCall(PetscSegBufferDestroy(&x->segrecvint));
    395:   PetscCall(PetscSegBufferDestroy(&x->segrecvscalar));
    396:   PetscCall(PetscSegBufferDestroy(&x->segrecvframe));
    397:   PetscFunctionReturn(PETSC_SUCCESS);
    398: }

    400: static PetscErrorCode VecSetFromOptions_MPI(Vec X, PetscOptionItems *PetscOptionsObject)
    401: {
    402: #if !defined(PETSC_HAVE_MPIUNI)
    403:   PetscBool flg = PETSC_FALSE, set;

    405:   PetscFunctionBegin;
    406:   PetscOptionsHeadBegin(PetscOptionsObject, "VecMPI Options");
    407:   PetscCall(PetscOptionsBool("-vec_assembly_legacy", "Use MPI 1 version of assembly", "", flg, &flg, &set));
    408:   if (set) {
    409:     X->ops->assemblybegin = flg ? VecAssemblyBegin_MPI : VecAssemblyBegin_MPI_BTS;
    410:     X->ops->assemblyend   = flg ? VecAssemblyEnd_MPI : VecAssemblyEnd_MPI_BTS;
    411:   }
    412:   PetscOptionsHeadEnd();
    413: #else
    414:   PetscFunctionBegin;
    415:   X->ops->assemblybegin = VecAssemblyBegin_MPI;
    416:   X->ops->assemblyend   = VecAssemblyEnd_MPI;
    417: #endif
    418:   PetscFunctionReturn(PETSC_SUCCESS);
    419: }

    421: static PetscErrorCode VecGetLocalToGlobalMapping_MPI_VecGhost(Vec X, ISLocalToGlobalMapping *ltg)
    422: {
    423:   PetscInt       *indices, n, nghost, rstart, i;
    424:   IS              ghostis;
    425:   const PetscInt *ghostidx;

    427:   PetscFunctionBegin;
    428:   *ltg = X->map->mapping;
    429:   if (X->map->mapping) PetscFunctionReturn(PETSC_SUCCESS);
    430:   PetscCall(VecGhostGetGhostIS(X, &ghostis));
    431:   if (!ghostis) PetscFunctionReturn(PETSC_SUCCESS);
    432:   PetscCall(ISGetLocalSize(ghostis, &nghost));
    433:   PetscCall(VecGetLocalSize(X, &n));
    434:   PetscCall(ISGetIndices(ghostis, &ghostidx));
    435:   /* set local to global mapping for ghosted vector */
    436:   PetscCall(PetscMalloc1(n + nghost, &indices));
    437:   PetscCall(VecGetOwnershipRange(X, &rstart, NULL));
    438:   for (i = 0; i < n; i++) indices[i] = rstart + i;
    439:   PetscCall(PetscArraycpy(indices + n, ghostidx, nghost));
    440:   PetscCall(ISRestoreIndices(ghostis, &ghostidx));
    441:   PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, 1, n + nghost, indices, PETSC_OWN_POINTER, &X->map->mapping));
    442:   *ltg = X->map->mapping;
    443:   PetscFunctionReturn(PETSC_SUCCESS);
    444: }

    446: static struct _VecOps DvOps = {
    447:   PetscDesignatedInitializer(duplicate, VecDuplicate_MPI), /* 1 */
    448:   PetscDesignatedInitializer(duplicatevecs, VecDuplicateVecs_Default),
    449:   PetscDesignatedInitializer(destroyvecs, VecDestroyVecs_Default),
    450:   PetscDesignatedInitializer(dot, VecDot_MPI),
    451:   PetscDesignatedInitializer(mdot, VecMDot_MPI),
    452:   PetscDesignatedInitializer(norm, VecNorm_MPI),
    453:   PetscDesignatedInitializer(tdot, VecTDot_MPI),
    454:   PetscDesignatedInitializer(mtdot, VecMTDot_MPI),
    455:   PetscDesignatedInitializer(scale, VecScale_Seq),
    456:   PetscDesignatedInitializer(copy, VecCopy_Seq), /* 10 */
    457:   PetscDesignatedInitializer(set, VecSet_Seq),
    458:   PetscDesignatedInitializer(swap, VecSwap_Seq),
    459:   PetscDesignatedInitializer(axpy, VecAXPY_Seq),
    460:   PetscDesignatedInitializer(axpby, VecAXPBY_Seq),
    461:   PetscDesignatedInitializer(maxpy, VecMAXPY_Seq),
    462:   PetscDesignatedInitializer(aypx, VecAYPX_Seq),
    463:   PetscDesignatedInitializer(waxpy, VecWAXPY_Seq),
    464:   PetscDesignatedInitializer(axpbypcz, VecAXPBYPCZ_Seq),
    465:   PetscDesignatedInitializer(pointwisemult, VecPointwiseMult_Seq),
    466:   PetscDesignatedInitializer(pointwisedivide, VecPointwiseDivide_Seq),
    467:   PetscDesignatedInitializer(setvalues, VecSetValues_MPI), /* 20 */
    468:   PetscDesignatedInitializer(assemblybegin, VecAssemblyBegin_MPI_BTS),
    469:   PetscDesignatedInitializer(assemblyend, VecAssemblyEnd_MPI_BTS),
    470:   PetscDesignatedInitializer(getarray, NULL),
    471:   PetscDesignatedInitializer(getsize, VecGetSize_MPI),
    472:   PetscDesignatedInitializer(getlocalsize, VecGetSize_Seq),
    473:   PetscDesignatedInitializer(restorearray, NULL),
    474:   PetscDesignatedInitializer(max, VecMax_MPI),
    475:   PetscDesignatedInitializer(min, VecMin_MPI),
    476:   PetscDesignatedInitializer(setrandom, VecSetRandom_Seq),
    477:   PetscDesignatedInitializer(setoption, VecSetOption_MPI),
    478:   PetscDesignatedInitializer(setvaluesblocked, VecSetValuesBlocked_MPI),
    479:   PetscDesignatedInitializer(destroy, VecDestroy_MPI),
    480:   PetscDesignatedInitializer(view, VecView_MPI),
    481:   PetscDesignatedInitializer(placearray, VecPlaceArray_MPI),
    482:   PetscDesignatedInitializer(replacearray, VecReplaceArray_Seq),
    483:   PetscDesignatedInitializer(dot_local, VecDot_Seq),
    484:   PetscDesignatedInitializer(tdot_local, VecTDot_Seq),
    485:   PetscDesignatedInitializer(norm_local, VecNorm_Seq),
    486:   PetscDesignatedInitializer(mdot_local, VecMDot_Seq),
    487:   PetscDesignatedInitializer(mtdot_local, VecMTDot_Seq),
    488:   PetscDesignatedInitializer(load, VecLoad_Default),
    489:   PetscDesignatedInitializer(reciprocal, VecReciprocal_Default),
    490:   PetscDesignatedInitializer(conjugate, VecConjugate_Seq),
    491:   PetscDesignatedInitializer(setlocaltoglobalmapping, NULL),
    492:   PetscDesignatedInitializer(getlocaltoglobalmapping, VecGetLocalToGlobalMapping_MPI_VecGhost),
    493:   PetscDesignatedInitializer(setvalueslocal, NULL),
    494:   PetscDesignatedInitializer(resetarray, VecResetArray_MPI),
    495:   PetscDesignatedInitializer(setfromoptions, VecSetFromOptions_MPI), /*set from options */
    496:   PetscDesignatedInitializer(maxpointwisedivide, VecMaxPointwiseDivide_Seq),
    497:   PetscDesignatedInitializer(pointwisemax, VecPointwiseMax_Seq),
    498:   PetscDesignatedInitializer(pointwisemaxabs, VecPointwiseMaxAbs_Seq),
    499:   PetscDesignatedInitializer(pointwisemin, VecPointwiseMin_Seq),
    500:   PetscDesignatedInitializer(getvalues, VecGetValues_MPI),
    501:   PetscDesignatedInitializer(sqrt, NULL),
    502:   PetscDesignatedInitializer(abs, NULL),
    503:   PetscDesignatedInitializer(exp, NULL),
    504:   PetscDesignatedInitializer(log, NULL),
    505:   PetscDesignatedInitializer(shift, NULL),
    506:   PetscDesignatedInitializer(create, NULL), /* really? */
    507:   PetscDesignatedInitializer(stridegather, VecStrideGather_Default),
    508:   PetscDesignatedInitializer(stridescatter, VecStrideScatter_Default),
    509:   PetscDesignatedInitializer(dotnorm2, NULL),
    510:   PetscDesignatedInitializer(getsubvector, NULL),
    511:   PetscDesignatedInitializer(restoresubvector, NULL),
    512:   PetscDesignatedInitializer(getarrayread, NULL),
    513:   PetscDesignatedInitializer(restorearrayread, NULL),
    514:   PetscDesignatedInitializer(stridesubsetgather, VecStrideSubSetGather_Default),
    515:   PetscDesignatedInitializer(stridesubsetscatter, VecStrideSubSetScatter_Default),
    516:   PetscDesignatedInitializer(viewnative, VecView_MPI),
    517:   PetscDesignatedInitializer(loadnative, NULL),
    518:   PetscDesignatedInitializer(createlocalvector, NULL),
    519:   PetscDesignatedInitializer(getlocalvector, NULL),
    520:   PetscDesignatedInitializer(restorelocalvector, NULL),
    521:   PetscDesignatedInitializer(getlocalvectorread, NULL),
    522:   PetscDesignatedInitializer(restorelocalvectorread, NULL),
    523:   PetscDesignatedInitializer(bindtocpu, NULL),
    524:   PetscDesignatedInitializer(getarraywrite, NULL),
    525:   PetscDesignatedInitializer(restorearraywrite, NULL),
    526:   PetscDesignatedInitializer(getarrayandmemtype, NULL),
    527:   PetscDesignatedInitializer(restorearrayandmemtype, NULL),
    528:   PetscDesignatedInitializer(getarrayreadandmemtype, NULL),
    529:   PetscDesignatedInitializer(restorearrayreadandmemtype, NULL),
    530:   PetscDesignatedInitializer(getarraywriteandmemtype, NULL),
    531:   PetscDesignatedInitializer(restorearraywriteandmemtype, NULL),
    532:   PetscDesignatedInitializer(concatenate, NULL),
    533:   PetscDesignatedInitializer(sum, NULL),
    534:   PetscDesignatedInitializer(setpreallocationcoo, VecSetPreallocationCOO_MPI),
    535:   PetscDesignatedInitializer(setvaluescoo, VecSetValuesCOO_MPI),
    536:   PetscDesignatedInitializer(errorwnorm, NULL),
    537:   PetscDesignatedInitializer(maxpby, NULL),
    538: };

    540: /*
    541:     VecCreate_MPI_Private - Basic create routine called by VecCreate_MPI() (i.e. VecCreateMPI()),
    542:     VecCreateMPIWithArray(), VecCreate_Shared() (i.e. VecCreateShared()), VecCreateGhost(),
    543:     VecDuplicate_MPI(), VecCreateGhostWithArray(), VecDuplicate_MPI(), and VecDuplicate_Shared()

    545:     If alloc is true and array is NULL then this routine allocates the space, otherwise
    546:     no space is allocated.
    547: */
    548: PetscErrorCode VecCreate_MPI_Private(Vec v, PetscBool alloc, PetscInt nghost, const PetscScalar array[])
    549: {
    550:   Vec_MPI  *s;
    551:   PetscBool mdot_use_gemv  = PETSC_TRUE;
    552:   PetscBool maxpy_use_gemv = PETSC_FALSE; // default is false as we saw bad performance with vendors' GEMV with tall skinny matrices.

    554:   PetscFunctionBegin;
    555:   PetscCall(PetscNew(&s));
    556:   v->data   = (void *)s;
    557:   v->ops[0] = DvOps;

    559:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-vec_mdot_use_gemv", &mdot_use_gemv, NULL));
    560:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-vec_maxpy_use_gemv", &maxpy_use_gemv, NULL));

    562:   // allocate multiple vectors together
    563:   if (mdot_use_gemv || maxpy_use_gemv) v->ops[0].duplicatevecs = VecDuplicateVecs_MPI_GEMV;

    565:   if (mdot_use_gemv) {
    566:     v->ops[0].mdot        = VecMDot_MPI_GEMV;
    567:     v->ops[0].mdot_local  = VecMDot_Seq_GEMV;
    568:     v->ops[0].mtdot       = VecMTDot_MPI_GEMV;
    569:     v->ops[0].mtdot_local = VecMTDot_Seq_GEMV;
    570:   }
    571:   if (maxpy_use_gemv) v->ops[0].maxpy = VecMAXPY_Seq_GEMV;

    573:   s->nghost      = nghost;
    574:   v->petscnative = PETSC_TRUE;
    575:   if (array) v->offloadmask = PETSC_OFFLOAD_CPU;

    577:   PetscCall(PetscLayoutSetUp(v->map));

    579:   s->array           = (PetscScalar *)array;
    580:   s->array_allocated = NULL;
    581:   if (alloc && !array) {
    582:     PetscInt n = v->map->n + nghost;
    583:     PetscCall(PetscCalloc1(n, &s->array));
    584:     s->array_allocated = s->array;
    585:     PetscCall(PetscObjectComposedDataSetReal((PetscObject)v, NormIds[NORM_2], 0));
    586:     PetscCall(PetscObjectComposedDataSetReal((PetscObject)v, NormIds[NORM_1], 0));
    587:     PetscCall(PetscObjectComposedDataSetReal((PetscObject)v, NormIds[NORM_INFINITY], 0));
    588:   }

    590:   /* By default parallel vectors do not have local representation */
    591:   s->localrep    = NULL;
    592:   s->localupdate = NULL;
    593:   s->ghost       = NULL;

    595:   v->stash.insertmode  = NOT_SET_VALUES;
    596:   v->bstash.insertmode = NOT_SET_VALUES;
    597:   /* create the stashes. The block-size for bstash is set later when
    598:      VecSetValuesBlocked is called.
    599:   */
    600:   PetscCall(VecStashCreate_Private(PetscObjectComm((PetscObject)v), 1, &v->stash));
    601:   PetscCall(VecStashCreate_Private(PetscObjectComm((PetscObject)v), PetscAbs(v->map->bs), &v->bstash));

    603: #if defined(PETSC_HAVE_MATLAB)
    604:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", VecMatlabEnginePut_Default));
    605:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", VecMatlabEngineGet_Default));
    606: #endif
    607:   PetscCall(PetscObjectChangeTypeName((PetscObject)v, VECMPI));
    608:   PetscFunctionReturn(PETSC_SUCCESS);
    609: }

    611: /*
    612:   Create a VECMPI with the given layout and array

    614:   Collective

    616:   Input Parameter:
    617: + map   - the layout
    618: - array - the array on host

    620:   Output Parameter:
    621: . V  - The vector object
    622: */
    623: PetscErrorCode VecCreateMPIWithLayoutAndArray_Private(PetscLayout map, const PetscScalar array[], Vec *V)
    624: {
    625:   PetscFunctionBegin;
    626:   PetscCall(VecCreateWithLayout_Private(map, V));
    627:   PetscCall(VecCreate_MPI_Private(*V, PETSC_FALSE, 0, array));
    628:   PetscFunctionReturn(PETSC_SUCCESS);
    629: }

    631: /*MC
    632:    VECMPI - VECMPI = "mpi" - The basic parallel vector

    634:    Options Database Key:
    635: . -vec_type mpi - sets the vector type to `VECMPI` during a call to `VecSetFromOptions()`

    637:   Level: beginner

    639: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecSetType()`, `VecSetFromOptions()`, `VecCreateMPIWithArray()`, `VECMPI`, `VecType`, `VecCreateMPI()`, `VecCreateMPI()`
    640: M*/

    642: PetscErrorCode VecCreate_MPI(Vec vv)
    643: {
    644:   PetscFunctionBegin;
    645:   PetscCall(VecCreate_MPI_Private(vv, PETSC_TRUE, 0, NULL));
    646:   PetscFunctionReturn(PETSC_SUCCESS);
    647: }

    649: /*MC
    650:    VECSTANDARD = "standard" - A `VECSEQ` on one process and `VECMPI` on more than one process

    652:    Options Database Key:
    653: . -vec_type standard - sets a vector type to standard on calls to `VecSetFromOptions()`

    655:   Level: beginner

    657: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeq()`, `VecCreateMPI()`
    658: M*/

    660: PETSC_EXTERN PetscErrorCode VecCreate_Standard(Vec v)
    661: {
    662:   PetscMPIInt size;

    664:   PetscFunctionBegin;
    665:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
    666:   if (size == 1) {
    667:     PetscCall(VecSetType(v, VECSEQ));
    668:   } else {
    669:     PetscCall(VecSetType(v, VECMPI));
    670:   }
    671:   PetscFunctionReturn(PETSC_SUCCESS);
    672: }

    674: /*@
    675:   VecCreateMPIWithArray - Creates a parallel, array-style vector,
    676:   where the user provides the array space to store the vector values.

    678:   Collective

    680:   Input Parameters:
    681: + comm  - the MPI communicator to use
    682: . bs    - block size, same meaning as `VecSetBlockSize()`
    683: . n     - local vector length, cannot be `PETSC_DECIDE`
    684: . N     - global vector length (or `PETSC_DETERMINE` to have calculated)
    685: - array - the user provided array to store the vector values

    687:   Output Parameter:
    688: . vv - the vector

    690:   Level: intermediate

    692:   Notes:
    693:   Use `VecDuplicate()` or `VecDuplicateVecs()` to form additional vectors of the
    694:   same type as an existing vector.

    696:   If the user-provided array is `NULL`, then `VecPlaceArray()` can be used
    697:   at a later stage to SET the array for storing the vector values.

    699:   PETSc does NOT free `array` when the vector is destroyed via `VecDestroy()`.

    701:   The user should not free `array` until the vector is destroyed.

    703: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeqWithArray()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateGhost()`,
    704:           `VecCreateMPI()`, `VecCreateGhostWithArray()`, `VecPlaceArray()`
    705: @*/
    706: PetscErrorCode VecCreateMPIWithArray(MPI_Comm comm, PetscInt bs, PetscInt n, PetscInt N, const PetscScalar array[], Vec *vv)
    707: {
    708:   PetscFunctionBegin;
    709:   PetscCheck(n != PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must set local size of vector");
    710:   PetscCall(PetscSplitOwnership(comm, &n, &N));
    711:   PetscCall(VecCreate(comm, vv));
    712:   PetscCall(VecSetSizes(*vv, n, N));
    713:   PetscCall(VecSetBlockSize(*vv, bs));
    714:   PetscCall(VecCreate_MPI_Private(*vv, PETSC_FALSE, 0, array));
    715:   PetscFunctionReturn(PETSC_SUCCESS);
    716: }

    718: /*@
    719:   VecCreateGhostWithArray - Creates a parallel vector with ghost padding on each processor;
    720:   the caller allocates the array space.

    722:   Collective

    724:   Input Parameters:
    725: + comm   - the MPI communicator to use
    726: . n      - local vector length
    727: . N      - global vector length (or `PETSC_DETERMINE` to have calculated if `n` is given)
    728: . nghost - number of local ghost points
    729: . ghosts - global indices of ghost points (or `NULL` if not needed), these do not need to be in increasing order (sorted)
    730: - array  - the space to store the vector values (as long as n + nghost)

    732:   Output Parameter:
    733: . vv - the global vector representation (without ghost points as part of vector)

    735:   Level: advanced

    737:   Notes:
    738:   Use `VecGhostGetLocalForm()` to access the local, ghosted representation
    739:   of the vector.

    741:   This also automatically sets the `ISLocalToGlobalMapping()` for this vector.

    743: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`,
    744:           `VecCreateGhost()`, `VecCreateSeqWithArray()`, `VecCreateMPIWithArray()`,
    745:           `VecCreateGhostBlock()`, `VecCreateGhostBlockWithArray()`, `VecMPISetGhost()`
    746: @*/
    747: PetscErrorCode VecCreateGhostWithArray(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[], const PetscScalar array[], Vec *vv)
    748: {
    749:   Vec_MPI     *w;
    750:   PetscScalar *larray;
    751:   IS           from, to;

    753:   PetscFunctionBegin;
    754:   *vv = NULL;
    755:   PetscCheck(nghost != PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must set local ghost size");
    756:   PetscCheck(nghost >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Ghost length must be >= 0");
    757:   PetscCall(PetscSplitOwnership(comm, &n, &N));
    758:   /* Create global representation */
    759:   PetscCall(VecCreate(comm, vv));
    760:   PetscCall(VecSetSizes(*vv, n, N));
    761:   PetscCall(VecCreate_MPI_Private(*vv, PETSC_TRUE, nghost, array));
    762:   w = (Vec_MPI *)(*vv)->data;
    763:   /* Create local representation */
    764:   PetscCall(VecGetArray(*vv, &larray));
    765:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, n + nghost, larray, &w->localrep));
    766:   PetscCall(VecRestoreArray(*vv, &larray));

    768:   /*
    769:        Create scatter context for scattering (updating) ghost values
    770:   */
    771:   PetscCall(ISCreateGeneral(comm, nghost, ghosts, PETSC_COPY_VALUES, &from));
    772:   PetscCall(ISCreateStride(PETSC_COMM_SELF, nghost, n, 1, &to));
    773:   PetscCall(VecScatterCreate(*vv, from, w->localrep, to, &w->localupdate));
    774:   PetscCall(ISDestroy(&to));

    776:   w->ghost                            = from;
    777:   (*vv)->ops->getlocaltoglobalmapping = VecGetLocalToGlobalMapping_MPI_VecGhost;
    778:   PetscFunctionReturn(PETSC_SUCCESS);
    779: }

    781: /*@
    782:   VecGhostGetGhostIS - Return ghosting indices of a ghost vector

    784:   Input Parameters:
    785: . X - ghost vector

    787:   Output Parameter:
    788: . ghost - ghosting indices

    790:   Level: beginner

    792: .seealso: `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`
    793: @*/
    794: PetscErrorCode VecGhostGetGhostIS(Vec X, IS *ghost)
    795: {
    796:   Vec_MPI  *w;
    797:   PetscBool flg;

    799:   PetscFunctionBegin;
    801:   PetscAssertPointer(ghost, 2);
    802:   PetscCall(PetscObjectTypeCompare((PetscObject)X, VECMPI, &flg));
    803:   PetscCheck(flg, PetscObjectComm((PetscObject)X), PETSC_ERR_ARG_WRONGSTATE, "VecGhostGetGhostIS was not supported for vec type %s", ((PetscObject)X)->type_name);
    804:   w      = (Vec_MPI *)X->data;
    805:   *ghost = w->ghost;
    806:   PetscFunctionReturn(PETSC_SUCCESS);
    807: }

    809: /*@
    810:   VecCreateGhost - Creates a parallel vector with ghost padding on each processor.

    812:   Collective

    814:   Input Parameters:
    815: + comm   - the MPI communicator to use
    816: . n      - local vector length
    817: . N      - global vector length (or `PETSC_DETERMINE` to have calculated if `n` is given)
    818: . nghost - number of local ghost points
    819: - ghosts - global indices of ghost points, these do not need to be in increasing order (sorted)

    821:   Output Parameter:
    822: . vv - the global vector representation (without ghost points as part of vector)

    824:   Level: advanced

    826:   Notes:
    827:   Use `VecGhostGetLocalForm()` to access the local, ghosted representation
    828:   of the vector.

    830:   This also automatically sets the `ISLocalToGlobalMapping()` for this vector.

    832: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeq()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateMPI()`,
    833:           `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`, `VecGhostUpdateBegin()`,
    834:           `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`, `VecGhostUpdateEnd()`,
    835:           `VecCreateGhostBlock()`, `VecCreateGhostBlockWithArray()`, `VecMPISetGhost()`

    837: @*/
    838: PetscErrorCode VecCreateGhost(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[], Vec *vv)
    839: {
    840:   PetscFunctionBegin;
    841:   PetscCall(VecCreateGhostWithArray(comm, n, N, nghost, ghosts, NULL, vv));
    842:   PetscFunctionReturn(PETSC_SUCCESS);
    843: }

    845: /*@
    846:   VecMPISetGhost - Sets the ghost points for an MPI ghost vector

    848:   Collective

    850:   Input Parameters:
    851: + vv     - the MPI vector
    852: . nghost - number of local ghost points
    853: - ghosts - global indices of ghost points, these do not need to be in increasing order (sorted)

    855:   Level: advanced

    857:   Notes:
    858:   Use `VecGhostGetLocalForm()` to access the local, ghosted representation
    859:   of the vector.

    861:   This also automatically sets the `ISLocalToGlobalMapping()` for this vector.

    863:   You must call this AFTER you have set the type of the vector (with` VecSetType()`) and the size (with `VecSetSizes()`).

    865: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeq()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateMPI()`,
    866:           `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`, `VecGhostUpdateBegin()`,
    867:           `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`, `VecGhostUpdateEnd()`,
    868:           `VecCreateGhostBlock()`, `VecCreateGhostBlockWithArray()`
    869: @*/
    870: PetscErrorCode VecMPISetGhost(Vec vv, PetscInt nghost, const PetscInt ghosts[])
    871: {
    872:   PetscBool flg;

    874:   PetscFunctionBegin;
    875:   PetscCall(PetscObjectTypeCompare((PetscObject)vv, VECMPI, &flg));
    876:   /* if already fully existent VECMPI then basically destroy it and rebuild with ghosting */
    877:   if (flg) {
    878:     PetscInt     n, N;
    879:     Vec_MPI     *w;
    880:     PetscScalar *larray;
    881:     IS           from, to;
    882:     MPI_Comm     comm;

    884:     PetscCall(PetscObjectGetComm((PetscObject)vv, &comm));
    885:     n = vv->map->n;
    886:     N = vv->map->N;
    887:     PetscUseTypeMethod(vv, destroy);
    888:     PetscCall(VecSetSizes(vv, n, N));
    889:     PetscCall(VecCreate_MPI_Private(vv, PETSC_TRUE, nghost, NULL));
    890:     w = (Vec_MPI *)vv->data;
    891:     /* Create local representation */
    892:     PetscCall(VecGetArray(vv, &larray));
    893:     PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, n + nghost, larray, &w->localrep));
    894:     PetscCall(VecRestoreArray(vv, &larray));

    896:     /*
    897:      Create scatter context for scattering (updating) ghost values
    898:      */
    899:     PetscCall(ISCreateGeneral(comm, nghost, ghosts, PETSC_COPY_VALUES, &from));
    900:     PetscCall(ISCreateStride(PETSC_COMM_SELF, nghost, n, 1, &to));
    901:     PetscCall(VecScatterCreate(vv, from, w->localrep, to, &w->localupdate));
    902:     PetscCall(ISDestroy(&to));

    904:     w->ghost                         = from;
    905:     vv->ops->getlocaltoglobalmapping = VecGetLocalToGlobalMapping_MPI_VecGhost;
    906:   } else {
    907:     PetscCheck(vv->ops->create != VecCreate_MPI, PetscObjectComm((PetscObject)vv), PETSC_ERR_ARG_WRONGSTATE, "Must set local or global size before setting ghosting");
    908:     PetscCheck(((PetscObject)vv)->type_name, PetscObjectComm((PetscObject)vv), PETSC_ERR_ARG_WRONGSTATE, "Must set type to VECMPI before ghosting");
    909:   }
    910:   PetscFunctionReturn(PETSC_SUCCESS);
    911: }

    913: /*@
    914:   VecCreateGhostBlockWithArray - Creates a parallel vector with ghost padding on each processor;
    915:   the caller allocates the array space. Indices in the ghost region are based on blocks.

    917:   Collective

    919:   Input Parameters:
    920: + comm   - the MPI communicator to use
    921: . bs     - block size
    922: . n      - local vector length
    923: . N      - global vector length (or `PETSC_DETERMINE` to have calculated if `n` is given)
    924: . nghost - number of local ghost blocks
    925: . ghosts - global indices of ghost blocks (or `NULL` if not needed), counts are by block not by index, these do not need to be in increasing order (sorted)
    926: - array  - the space to store the vector values (as long as `n + nghost*bs`)

    928:   Output Parameter:
    929: . vv - the global vector representation (without ghost points as part of vector)

    931:   Level: advanced

    933:   Notes:
    934:   Use `VecGhostGetLocalForm()` to access the local, ghosted representation
    935:   of the vector.

    937:   n is the local vector size (total local size not the number of blocks) while nghost
    938:   is the number of blocks in the ghost portion, i.e. the number of elements in the ghost
    939:   portion is bs*nghost

    941: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`,
    942:           `VecCreateGhost()`, `VecCreateSeqWithArray()`, `VecCreateMPIWithArray()`,
    943:           `VecCreateGhostWithArray()`, `VecCreateGhostBlock()`
    944: @*/
    945: PetscErrorCode VecCreateGhostBlockWithArray(MPI_Comm comm, PetscInt bs, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[], const PetscScalar array[], Vec *vv)
    946: {
    947:   Vec_MPI               *w;
    948:   PetscScalar           *larray;
    949:   IS                     from, to;
    950:   ISLocalToGlobalMapping ltog;
    951:   PetscInt               rstart, i, nb, *indices;

    953:   PetscFunctionBegin;
    954:   *vv = NULL;

    956:   PetscCheck(n != PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must set local size");
    957:   PetscCheck(nghost != PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must set local ghost size");
    958:   PetscCheck(nghost >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Ghost length must be >= 0");
    959:   PetscCheck(n % bs == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local size must be a multiple of block size");
    960:   PetscCall(PetscSplitOwnership(comm, &n, &N));
    961:   /* Create global representation */
    962:   PetscCall(VecCreate(comm, vv));
    963:   PetscCall(VecSetSizes(*vv, n, N));
    964:   PetscCall(VecSetBlockSize(*vv, bs));
    965:   PetscCall(VecCreate_MPI_Private(*vv, PETSC_TRUE, nghost * bs, array));
    966:   w = (Vec_MPI *)(*vv)->data;
    967:   /* Create local representation */
    968:   PetscCall(VecGetArray(*vv, &larray));
    969:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, bs, n + bs * nghost, larray, &w->localrep));
    970:   PetscCall(VecRestoreArray(*vv, &larray));

    972:   /*
    973:        Create scatter context for scattering (updating) ghost values
    974:   */
    975:   PetscCall(ISCreateBlock(comm, bs, nghost, ghosts, PETSC_COPY_VALUES, &from));
    976:   PetscCall(ISCreateStride(PETSC_COMM_SELF, bs * nghost, n, 1, &to));
    977:   PetscCall(VecScatterCreate(*vv, from, w->localrep, to, &w->localupdate));
    978:   PetscCall(ISDestroy(&to));
    979:   PetscCall(ISDestroy(&from));

    981:   /* set local to global mapping for ghosted vector */
    982:   nb = n / bs;
    983:   PetscCall(PetscMalloc1(nb + nghost, &indices));
    984:   PetscCall(VecGetOwnershipRange(*vv, &rstart, NULL));
    985:   rstart = rstart / bs;

    987:   for (i = 0; i < nb; i++) indices[i] = rstart + i;
    988:   for (i = 0; i < nghost; i++) indices[nb + i] = ghosts[i];

    990:   PetscCall(ISLocalToGlobalMappingCreate(comm, bs, nb + nghost, indices, PETSC_OWN_POINTER, &ltog));
    991:   PetscCall(VecSetLocalToGlobalMapping(*vv, ltog));
    992:   PetscCall(ISLocalToGlobalMappingDestroy(&ltog));
    993:   PetscFunctionReturn(PETSC_SUCCESS);
    994: }

    996: /*@
    997:   VecCreateGhostBlock - Creates a parallel vector with ghost padding on each processor.
    998:   The indicing of the ghost points is done with blocks.

    1000:   Collective

    1002:   Input Parameters:
    1003: + comm   - the MPI communicator to use
    1004: . bs     - the block size
    1005: . n      - local vector length
    1006: . N      - global vector length (or `PETSC_DETERMINE` to have calculated if `n` is given)
    1007: . nghost - number of local ghost blocks
    1008: - ghosts - global indices of ghost blocks, counts are by block, not by individual index, these do not need to be in increasing order (sorted)

    1010:   Output Parameter:
    1011: . vv - the global vector representation (without ghost points as part of vector)

    1013:   Level: advanced

    1015:   Notes:
    1016:   Use `VecGhostGetLocalForm()` to access the local, ghosted representation
    1017:   of the vector.

    1019:   `n` is the local vector size (total local size not the number of blocks) while `nghost`
    1020:   is the number of blocks in the ghost portion, i.e. the number of elements in the ghost
    1021:   portion is `bs*nghost`

    1023: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeq()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateMPI()`,
    1024:           `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`,
    1025:           `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`, `VecCreateGhostBlockWithArray()`
    1026: @*/
    1027: PetscErrorCode VecCreateGhostBlock(MPI_Comm comm, PetscInt bs, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[], Vec *vv)
    1028: {
    1029:   PetscFunctionBegin;
    1030:   PetscCall(VecCreateGhostBlockWithArray(comm, bs, n, N, nghost, ghosts, NULL, vv));
    1031:   PetscFunctionReturn(PETSC_SUCCESS);
    1032: }