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: // Duplicate a vector with the provided array on host (if not NULL) for the new vector.
     22: // In that case, the array should be long enough to take into account ghost points if any.
     23: // If array is NULL, the routine will allocate memory itself.
     24: PetscErrorCode VecDuplicateWithArray_MPI(Vec win, const PetscScalar *array, Vec *v)
     25: {
     26:   Vec_MPI *vw, *w = (Vec_MPI *)win->data;

     28:   PetscFunctionBegin;
     29:   PetscCall(VecCreateWithLayout_Private(win->map, v));

     31:   PetscCall(VecCreate_MPI_Private(*v, PETSC_TRUE, w->nghost, array)); // array could be NULL
     32:   vw           = (Vec_MPI *)(*v)->data;
     33:   (*v)->ops[0] = win->ops[0];

     35:   /* save local representation of the parallel vector (and scatter) if it exists */
     36:   if (w->localrep) {
     37:     PetscScalar *arr = ((Vec_MPI *)(*v)->data)->array;

     39:     PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, win->map->bs, win->map->n + w->nghost, arr, &vw->localrep));
     40:     vw->localrep->ops[0] = w->localrep->ops[0];

     42:     vw->localupdate = w->localupdate;
     43:     PetscCall(PetscObjectReference((PetscObject)vw->localupdate));

     45:     vw->ghost = w->ghost;
     46:     PetscCall(PetscObjectReference((PetscObject)vw->ghost));
     47:   }

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

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

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

     61: PetscErrorCode VecDuplicate_MPI(Vec win, Vec *v)
     62: {
     63:   PetscFunctionBegin;
     64:   PetscCall(VecDuplicateWithArray_MPI(win, NULL, v));
     65:   PetscFunctionReturn(PETSC_SUCCESS);
     66: }

     68: static PetscErrorCode VecDuplicateVecs_MPI_GEMV(Vec w, PetscInt m, Vec *V[])
     69: {
     70:   Vec_MPI *wmpi = (Vec_MPI *)w->data;

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

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

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

    108: static PetscErrorCode VecSetOption_MPI(Vec V, VecOption op, PetscBool flag)
    109: {
    110:   Vec_MPI *v = (Vec_MPI *)V->data;

    112:   PetscFunctionBegin;
    113:   switch (op) {
    114:   case VEC_IGNORE_OFF_PROC_ENTRIES:
    115:     V->stash.donotstash = flag;
    116:     break;
    117:   case VEC_IGNORE_NEGATIVE_INDICES:
    118:     V->stash.ignorenegidx = flag;
    119:     break;
    120:   case VEC_SUBSET_OFF_PROC_ENTRIES:
    121:     v->assembly_subset = flag;              /* See the same logic in MatAssembly wrt MAT_SUBSET_OFF_PROC_ENTRIES */
    122:     if (!v->assembly_subset) {              /* User indicates "do not reuse the communication pattern" */
    123:       PetscCall(VecAssemblyReset_MPI(V));   /* Reset existing pattern to free memory */
    124:       v->first_assembly_done = PETSC_FALSE; /* Mark the first assembly is not done */
    125:     }
    126:     break;
    127:   }
    128:   PetscFunctionReturn(PETSC_SUCCESS);
    129: }

    131: PetscErrorCode VecResetArray_MPI(Vec vin)
    132: {
    133:   Vec_MPI *v = (Vec_MPI *)vin->data;

    135:   PetscFunctionBegin;
    136:   v->array         = v->unplacedarray;
    137:   v->unplacedarray = NULL;
    138:   if (v->localrep) PetscCall(VecResetArray(v->localrep));
    139:   PetscFunctionReturn(PETSC_SUCCESS);
    140: }

    142: static PetscErrorCode VecAssemblySend_MPI_Private(MPI_Comm comm, const PetscMPIInt tag[], PetscMPIInt rankid, PetscMPIInt rank, void *sdata, MPI_Request req[], void *ctx)
    143: {
    144:   Vec                X   = (Vec)ctx;
    145:   Vec_MPI           *x   = (Vec_MPI *)X->data;
    146:   VecAssemblyHeader *hdr = (VecAssemblyHeader *)sdata;
    147:   PetscInt           bs  = X->map->bs;

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

    165: static PetscErrorCode VecAssemblyRecv_MPI_Private(MPI_Comm comm, const PetscMPIInt tag[], PetscMPIInt rank, void *rdata, MPI_Request req[], void *ctx)
    166: {
    167:   Vec                X   = (Vec)ctx;
    168:   Vec_MPI           *x   = (Vec_MPI *)X->data;
    169:   VecAssemblyHeader *hdr = (VecAssemblyHeader *)rdata;
    170:   PetscInt           bs  = X->map->bs;
    171:   VecAssemblyFrame  *frame;

    173:   PetscFunctionBegin;
    174:   PetscCall(PetscSegBufferGet(x->segrecvframe, 1, &frame));

    176:   if (hdr->count) {
    177:     PetscCall(PetscSegBufferGet(x->segrecvint, hdr->count, &frame->ints));
    178:     PetscCallMPI(MPIU_Irecv(frame->ints, hdr->count, MPIU_INT, rank, tag[0], comm, &req[0]));
    179:     PetscCall(PetscSegBufferGet(x->segrecvscalar, hdr->count, &frame->scalars));
    180:     PetscCallMPI(MPIU_Irecv(frame->scalars, hdr->count, MPIU_SCALAR, rank, tag[1], comm, &req[1]));
    181:     frame->pendings = 2;
    182:   } else {
    183:     frame->ints     = NULL;
    184:     frame->scalars  = NULL;
    185:     frame->pendings = 0;
    186:   }

    188:   if (hdr->bcount) {
    189:     PetscCall(PetscSegBufferGet(x->segrecvint, hdr->bcount, &frame->intb));
    190:     PetscCallMPI(MPIU_Irecv(frame->intb, hdr->bcount, MPIU_INT, rank, tag[2], comm, &req[2]));
    191:     PetscCall(PetscSegBufferGet(x->segrecvscalar, hdr->bcount * bs, &frame->scalarb));
    192:     PetscCallMPI(MPIU_Irecv(frame->scalarb, hdr->bcount * bs, MPIU_SCALAR, rank, tag[3], comm, &req[3]));
    193:     frame->pendingb = 2;
    194:   } else {
    195:     frame->intb     = NULL;
    196:     frame->scalarb  = NULL;
    197:     frame->pendingb = 0;
    198:   }
    199:   PetscFunctionReturn(PETSC_SUCCESS);
    200: }

    202: static PetscErrorCode VecAssemblyBegin_MPI_BTS(Vec X)
    203: {
    204:   Vec_MPI    *x = (Vec_MPI *)X->data;
    205:   MPI_Comm    comm;
    206:   PetscMPIInt i;
    207:   PetscInt    j, jb, bs;

    209:   PetscFunctionBegin;
    210:   if (X->stash.donotstash) PetscFunctionReturn(PETSC_SUCCESS);

    212:   PetscCall(PetscObjectGetComm((PetscObject)X, &comm));
    213:   PetscCall(VecGetBlockSize(X, &bs));
    214:   if (PetscDefined(USE_DEBUG)) {
    215:     InsertMode addv;

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

    222:   PetscCall(VecStashSortCompress_Private(&X->stash));
    223:   PetscCall(VecStashSortCompress_Private(&X->bstash));

    225:   if (!x->sendranks) {
    226:     PetscMPIInt nowners, bnowners, *owners, *bowners;
    227:     PetscInt    ntmp;

    229:     PetscCall(VecStashGetOwnerList_Private(&X->stash, X->map, &nowners, &owners));
    230:     PetscCall(VecStashGetOwnerList_Private(&X->bstash, X->map, &bnowners, &bowners));
    231:     PetscCall(PetscMergeMPIIntArray(nowners, owners, bnowners, bowners, &ntmp, &x->sendranks));
    232:     PetscCall(PetscMPIIntCast(ntmp, &x->nsendranks));
    233:     PetscCall(PetscFree(owners));
    234:     PetscCall(PetscFree(bowners));
    235:     PetscCall(PetscMalloc1(x->nsendranks, &x->sendhdr));
    236:     PetscCall(PetscCalloc1(x->nsendranks, &x->sendptrs));
    237:   }
    238:   for (i = 0, j = 0, jb = 0; i < x->nsendranks; i++) {
    239:     PetscMPIInt rank = x->sendranks[i];

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

    259:   if (!x->segrecvint) PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 1000, &x->segrecvint));
    260:   if (!x->segrecvscalar) PetscCall(PetscSegBufferCreate(sizeof(PetscScalar), 1000, &x->segrecvscalar));
    261:   if (!x->segrecvframe) PetscCall(PetscSegBufferCreate(sizeof(VecAssemblyFrame), 50, &x->segrecvframe));
    262:   if (x->first_assembly_done) { /* this is not the first assembly */
    263:     PetscMPIInt tag[4];
    264:     for (i = 0; i < 4; i++) PetscCall(PetscCommGetNewTag(comm, &tag[i]));
    265:     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));
    266:     for (i = 0; i < x->nrecvranks; i++) PetscCall(VecAssemblyRecv_MPI_Private(comm, tag, x->recvranks[i], x->recvhdr + i, x->recvreqs + 4 * i, X));
    267:     x->use_status = PETSC_TRUE;
    268:   } else { /* First time assembly */
    269:     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));
    270:     x->use_status = PETSC_FALSE;
    271:   }

    273:   /* The first_assembly_done flag is only meaningful when x->assembly_subset is set.
    274:      This line says when assembly_subset is set, then we mark that the first assembly is done.
    275:    */
    276:   x->first_assembly_done = x->assembly_subset;

    278:   {
    279:     PetscInt nstash, reallocs;
    280:     PetscCall(VecStashGetInfo_Private(&X->stash, &nstash, &reallocs));
    281:     PetscCall(PetscInfo(X, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
    282:     PetscCall(VecStashGetInfo_Private(&X->bstash, &nstash, &reallocs));
    283:     PetscCall(PetscInfo(X, "Block-Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
    284:   }
    285:   PetscFunctionReturn(PETSC_SUCCESS);
    286: }

    288: static PetscErrorCode VecAssemblyEnd_MPI_BTS(Vec X)
    289: {
    290:   Vec_MPI          *x  = (Vec_MPI *)X->data;
    291:   PetscInt          bs = X->map->bs;
    292:   PetscMPIInt       npending, *some_indices, r;
    293:   MPI_Status       *some_statuses;
    294:   PetscScalar      *xarray;
    295:   VecAssemblyFrame *frame;

    297:   PetscFunctionBegin;
    298:   if (X->stash.donotstash) {
    299:     X->stash.insertmode  = NOT_SET_VALUES;
    300:     X->bstash.insertmode = NOT_SET_VALUES;
    301:     PetscFunctionReturn(PETSC_SUCCESS);
    302:   }

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

    324:       npending--;
    325:       if (!blockmsg) { /* Scalar stash */
    326:         PetscMPIInt count;

    328:         if (--frame[i].pendings > 0) continue;
    329:         if (x->use_status) PetscCallMPI(MPI_Get_count(&some_statuses[ii], intmsg ? MPIU_INT : MPIU_SCALAR, &count));
    330:         else PetscCall(PetscMPIIntCast(x->recvhdr[i].count, &count));
    331:         for (j = 0, recvint = frame[i].ints, recvscalar = frame[i].scalars; j < count; j++, recvint++) {
    332:           PetscInt loc = *recvint - X->map->rstart;

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

    381:   X->stash.insertmode  = NOT_SET_VALUES;
    382:   X->bstash.insertmode = NOT_SET_VALUES;
    383:   PetscCall(VecStashScatterEnd_Private(&X->stash));
    384:   PetscCall(VecStashScatterEnd_Private(&X->bstash));
    385:   PetscFunctionReturn(PETSC_SUCCESS);
    386: }

    388: PetscErrorCode VecAssemblyReset_MPI(Vec X)
    389: {
    390:   Vec_MPI *x = (Vec_MPI *)X->data;

    392:   PetscFunctionBegin;
    393:   PetscCall(PetscFree(x->sendreqs));
    394:   PetscCall(PetscFree(x->recvreqs));
    395:   PetscCall(PetscFree(x->sendranks));
    396:   PetscCall(PetscFree(x->recvranks));
    397:   PetscCall(PetscFree(x->sendhdr));
    398:   PetscCall(PetscFree(x->recvhdr));
    399:   PetscCall(PetscFree(x->sendptrs));
    400:   PetscCall(PetscSegBufferDestroy(&x->segrecvint));
    401:   PetscCall(PetscSegBufferDestroy(&x->segrecvscalar));
    402:   PetscCall(PetscSegBufferDestroy(&x->segrecvframe));
    403:   PetscFunctionReturn(PETSC_SUCCESS);
    404: }

    406: static PetscErrorCode VecSetFromOptions_MPI(Vec X, PetscOptionItems PetscOptionsObject)
    407: {
    408: #if !defined(PETSC_HAVE_MPIUNI)
    409:   PetscBool flg = PETSC_FALSE, set;

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

    427: static PetscErrorCode VecGetLocalToGlobalMapping_MPI_VecGhost(Vec X, ISLocalToGlobalMapping *ltg)
    428: {
    429:   PetscInt       *indices, n, nghost, rstart, i;
    430:   IS              ghostis;
    431:   const PetscInt *ghostidx;

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

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

    545: /*
    546:     VecCreate_MPI_Private - Basic create routine called by VecCreate_MPI() (i.e. VecCreateMPI()),
    547:     VecCreateMPIWithArray(), VecCreate_Shared() (i.e. VecCreateShared()), VecCreateGhost(),
    548:     VecDuplicate_MPI(), VecCreateGhostWithArray(), VecDuplicate_MPI(), and VecDuplicate_Shared()

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

    559:   PetscFunctionBegin;
    560:   PetscCall(PetscNew(&s));
    561:   v->data   = (void *)s;
    562:   v->ops[0] = DvOps;

    564:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-vec_mdot_use_gemv", &mdot_use_gemv, NULL));
    565:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-vec_maxpy_use_gemv", &maxpy_use_gemv, NULL));

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

    570:   if (mdot_use_gemv) {
    571:     v->ops[0].mdot        = VecMDot_MPI_GEMV;
    572:     v->ops[0].mdot_local  = VecMDot_Seq_GEMV;
    573:     v->ops[0].mtdot       = VecMTDot_MPI_GEMV;
    574:     v->ops[0].mtdot_local = VecMTDot_Seq_GEMV;
    575:   }
    576:   if (maxpy_use_gemv) v->ops[0].maxpy = VecMAXPY_Seq_GEMV;

    578:   s->nghost      = nghost;
    579:   v->petscnative = PETSC_TRUE;
    580:   if (array) v->offloadmask = PETSC_OFFLOAD_CPU;

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

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

    595:   /* By default parallel vectors do not have local representation */
    596:   s->localrep    = NULL;
    597:   s->localupdate = NULL;
    598:   s->ghost       = NULL;

    600:   v->stash.insertmode  = NOT_SET_VALUES;
    601:   v->bstash.insertmode = NOT_SET_VALUES;
    602:   /* create the stashes. The block-size for bstash is set later when
    603:      VecSetValuesBlocked is called.
    604:   */
    605:   PetscCall(VecStashCreate_Private(PetscObjectComm((PetscObject)v), 1, &v->stash));
    606:   PetscCall(VecStashCreate_Private(PetscObjectComm((PetscObject)v), v->map->bs, &v->bstash));

    608: #if defined(PETSC_HAVE_MATLAB)
    609:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", VecMatlabEnginePut_Default));
    610:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", VecMatlabEngineGet_Default));
    611: #endif
    612:   PetscCall(PetscObjectChangeTypeName((PetscObject)v, VECMPI));
    613:   PetscFunctionReturn(PETSC_SUCCESS);
    614: }

    616: /*
    617:   Create a VECMPI with the given layout and array

    619:   Collective

    621:   Input Parameter:
    622: + map   - the layout
    623: - array - the array on host

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

    636: /*MC
    637:    VECMPI - VECMPI = "mpi" - The basic parallel vector

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

    642:   Level: beginner

    644: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecSetType()`, `VecSetFromOptions()`, `VecCreateMPIWithArray()`, `VECSEQ`, `VecCreateMPI()`
    645: M*/

    647: PetscErrorCode VecCreate_MPI(Vec vv)
    648: {
    649:   PetscFunctionBegin;
    650:   PetscCall(VecCreate_MPI_Private(vv, PETSC_TRUE, 0, NULL));
    651:   PetscFunctionReturn(PETSC_SUCCESS);
    652: }

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

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

    660:   Level: beginner

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

    665: PETSC_EXTERN PetscErrorCode VecCreate_Standard(Vec v)
    666: {
    667:   PetscMPIInt size;

    669:   PetscFunctionBegin;
    670:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
    671:   if (size == 1) PetscCall(VecSetType(v, VECSEQ));
    672:   else PetscCall(VecSetType(v, VECMPI));
    673:   PetscFunctionReturn(PETSC_SUCCESS);
    674: }

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

    680:   Collective

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

    689:   Output Parameter:
    690: . vv - the vector

    692:   Level: intermediate

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

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

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

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

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

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

    724:   Collective

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

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

    737:   Level: advanced

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

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

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

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

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

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

    783: /*@
    784:   VecGhostGetGhostIS - Return ghosting indices of a ghost vector

    786:   Input Parameters:
    787: . X - ghost vector

    789:   Output Parameter:
    790: . ghost - ghosting indices

    792:   Level: beginner

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

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

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

    814:   Collective

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

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

    826:   Level: advanced

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

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

    834: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeq()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateMPI()`,
    835:           `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`, `VecGhostUpdateBegin()`,
    836:           `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`, `VecGhostUpdateEnd()`,
    837:           `VecCreateGhostBlock()`, `VecCreateGhostBlockWithArray()`, `VecMPISetGhost()`
    838: @*/
    839: PetscErrorCode VecCreateGhost(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[], Vec *vv)
    840: {
    841:   PetscFunctionBegin;
    842:   PetscCall(VecCreateGhostWithArray(comm, n, N, nghost, ghosts, NULL, vv));
    843:   PetscFunctionReturn(PETSC_SUCCESS);
    844: }

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

    849:   Collective

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

    856:   Level: advanced

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

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

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

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

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

    887:     PetscCall(PetscObjectGetComm((PetscObject)vv, &comm));
    888:     PetscCall(VecGetType(vv, &vtype));             // vtype might be VECMPIKOKKOS
    889:     PetscCall(PetscStrallocpy(vtype, &savedType)); // save it as *vtype could be changed by VecCreate_MPI_Private()

    891:     n = vv->map->n;
    892:     N = vv->map->N;
    893:     PetscUseTypeMethod(vv, destroy);
    894:     PetscCall(VecSetSizes(vv, n, N));
    895:     PetscCall(VecCreate_MPI_Private(vv, PETSC_TRUE, nghost, NULL));
    896:     w = (Vec_MPI *)vv->data;
    897:     /* Create local representation */
    898:     PetscCall(VecGetArray(vv, &larray));
    899:     PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, n + nghost, larray, &w->localrep));
    900:     PetscCall(VecRestoreArray(vv, &larray));

    902:     /*
    903:      Create scatter context for scattering (updating) ghost values
    904:      */
    905:     PetscCall(ISCreateGeneral(comm, nghost, ghosts, PETSC_COPY_VALUES, &from));
    906:     PetscCall(ISCreateStride(PETSC_COMM_SELF, nghost, n, 1, &to));
    907:     PetscCall(VecScatterCreate(vv, from, w->localrep, to, &w->localupdate));
    908:     PetscCall(ISDestroy(&to));

    910:     w->ghost                         = from;
    911:     vv->ops->getlocaltoglobalmapping = VecGetLocalToGlobalMapping_MPI_VecGhost;
    912:     PetscCall(VecSetType(vv, savedType));
    913:     PetscCall(PetscFree(savedType));
    914:   } else {
    915:     PetscCheck(vv->ops->create != VecCreate_MPI, PetscObjectComm((PetscObject)vv), PETSC_ERR_ARG_WRONGSTATE, "Must set local or global size before setting ghosting");
    916:     PetscCheck(((PetscObject)vv)->type_name, PetscObjectComm((PetscObject)vv), PETSC_ERR_ARG_WRONGSTATE, "Must set type to VECMPI before ghosting");
    917:   }
    918:   PetscFunctionReturn(PETSC_SUCCESS);
    919: }

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

    925:   Collective

    927:   Input Parameters:
    928: + comm   - the MPI communicator to use
    929: . bs     - block size
    930: . n      - local vector length
    931: . N      - global vector length (or `PETSC_DETERMINE` to have calculated if `n` is given)
    932: . nghost - number of local ghost blocks
    933: . 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)
    934: - array  - the space to store the vector values (as long as $n + nghost*bs$)

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

    939:   Level: advanced

    941:   Notes:
    942:   Use `VecGhostGetLocalForm()` to access the local, ghosted representation
    943:   of the vector.

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

    949: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`,
    950:           `VecCreateGhost()`, `VecCreateSeqWithArray()`, `VecCreateMPIWithArray()`,
    951:           `VecCreateGhostWithArray()`, `VecCreateGhostBlock()`, `VecGhostUpdateBegin()`, `VecGhostUpdateEnd()`
    952: @*/
    953: PetscErrorCode VecCreateGhostBlockWithArray(MPI_Comm comm, PetscInt bs, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[], const PetscScalar array[], Vec *vv)
    954: {
    955:   Vec_MPI               *w;
    956:   PetscScalar           *larray;
    957:   IS                     from, to;
    958:   ISLocalToGlobalMapping ltog;
    959:   PetscInt               rstart, i, nb, *indices;

    961:   PetscFunctionBegin;
    962:   *vv = NULL;

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

    980:   /*
    981:        Create scatter context for scattering (updating) ghost values
    982:   */
    983:   PetscCall(ISCreateBlock(comm, bs, nghost, ghosts, PETSC_COPY_VALUES, &from));
    984:   PetscCall(ISCreateStride(PETSC_COMM_SELF, bs * nghost, n, 1, &to));
    985:   PetscCall(VecScatterCreate(*vv, from, w->localrep, to, &w->localupdate));
    986:   PetscCall(ISDestroy(&to));
    987:   PetscCall(ISDestroy(&from));

    989:   /* set local to global mapping for ghosted vector */
    990:   nb = n / bs;
    991:   PetscCall(PetscMalloc1(nb + nghost, &indices));
    992:   PetscCall(VecGetOwnershipRange(*vv, &rstart, NULL));
    993:   rstart = rstart / bs;

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

    998:   PetscCall(ISLocalToGlobalMappingCreate(comm, bs, nb + nghost, indices, PETSC_OWN_POINTER, &ltog));
    999:   PetscCall(VecSetLocalToGlobalMapping(*vv, ltog));
    1000:   PetscCall(ISLocalToGlobalMappingDestroy(&ltog));
    1001:   PetscFunctionReturn(PETSC_SUCCESS);
    1002: }

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

    1008:   Collective

    1010:   Input Parameters:
    1011: + comm   - the MPI communicator to use
    1012: . bs     - the block size
    1013: . n      - local vector length
    1014: . N      - global vector length (or `PETSC_DETERMINE` to have calculated if `n` is given)
    1015: . nghost - number of local ghost blocks
    1016: - ghosts - global indices of ghost blocks, counts are by block, not by individual index, these do not need to be in increasing order (sorted)

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

    1021:   Level: advanced

    1023:   Notes:
    1024:   Use `VecGhostGetLocalForm()` to access the local, ghosted representation
    1025:   of the vector.

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

    1031: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreateSeq()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateMPI()`,
    1032:           `VecGhostGetLocalForm()`, `VecGhostRestoreLocalForm()`, `VecGhostUpdateBegin()`, `VecGhostUpdateEnd()`,
    1033:           `VecCreateGhostWithArray()`, `VecCreateMPIWithArray()`, `VecCreateGhostBlockWithArray()`
    1034: @*/
    1035: PetscErrorCode VecCreateGhostBlock(MPI_Comm comm, PetscInt bs, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[], Vec *vv)
    1036: {
    1037:   PetscFunctionBegin;
    1038:   PetscCall(VecCreateGhostBlockWithArray(comm, bs, n, N, nghost, ghosts, NULL, vv));
    1039:   PetscFunctionReturn(PETSC_SUCCESS);
    1040: }