Actual source code: rvector.c

      1: /*
      2:      Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
      3:    These are the vector functions the user calls.
      4: */
      5: #include <petsc/private/vecimpl.h>

      7: PetscInt VecGetSubVectorSavedStateId = -1;

      9: #if PetscDefined(USE_DEBUG)
     10: // this is a no-op '0' macro in optimized builds
     11: PetscErrorCode VecValidValues_Internal(Vec vec, PetscInt argnum, PetscBool begin)
     12: {
     13:   PetscFunctionBegin;
     14:   if (vec->petscnative || vec->ops->getarray) {
     15:     PetscInt           n;
     16:     const PetscScalar *x;
     17:     PetscOffloadMask   mask;

     19:     PetscCall(VecGetOffloadMask(vec, &mask));
     20:     if (!PetscOffloadHost(mask)) PetscFunctionReturn(PETSC_SUCCESS);
     21:     PetscCall(VecGetLocalSize(vec, &n));
     22:     PetscCall(VecGetArrayRead(vec, &x));
     23:     for (PetscInt i = 0; i < n; i++) {
     24:       if (begin) {
     25:         PetscCheck(!PetscIsInfOrNanScalar(x[i]), PETSC_COMM_SELF, PETSC_ERR_FP, "Vec entry at local location %" PetscInt_FMT " is not-a-number or infinite at beginning of function: Parameter number %" PetscInt_FMT, i, argnum);
     26:       } else {
     27:         PetscCheck(!PetscIsInfOrNanScalar(x[i]), PETSC_COMM_SELF, PETSC_ERR_FP, "Vec entry at local location %" PetscInt_FMT " is not-a-number or infinite at end of function: Parameter number %" PetscInt_FMT, i, argnum);
     28:       }
     29:     }
     30:     PetscCall(VecRestoreArrayRead(vec, &x));
     31:   }
     32:   PetscFunctionReturn(PETSC_SUCCESS);
     33: }
     34: #endif

     36: /*@
     37:   VecMaxPointwiseDivide - Computes the maximum of the componentwise division `max = max_i abs(x[i]/y[i])`.

     39:   Logically Collective

     41:   Input Parameters:
     42: + x - the numerators
     43: - y - the denominators

     45:   Output Parameter:
     46: . max - the result

     48:   Level: advanced

     50:   Notes:
     51:   `x` and `y` may be the same vector

     53:   if a particular `y[i]` is zero, it is treated as 1 in the above formula

     55: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`
     56: @*/
     57: PetscErrorCode VecMaxPointwiseDivide(Vec x, Vec y, PetscReal *max)
     58: {
     59:   PetscFunctionBegin;
     62:   PetscAssertPointer(max, 3);
     65:   PetscCheckSameTypeAndComm(x, 1, y, 2);
     66:   VecCheckSameSize(x, 1, y, 2);
     67:   VecCheckAssembled(x);
     68:   VecCheckAssembled(y);
     69:   PetscCall(VecLockReadPush(x));
     70:   PetscCall(VecLockReadPush(y));
     71:   PetscUseTypeMethod(x, maxpointwisedivide, y, max);
     72:   PetscCall(VecLockReadPop(x));
     73:   PetscCall(VecLockReadPop(y));
     74:   PetscFunctionReturn(PETSC_SUCCESS);
     75: }

     77: /*@
     78:   VecDot - Computes the vector dot product.

     80:   Collective

     82:   Input Parameters:
     83: + x - first vector
     84: - y - second vector

     86:   Output Parameter:
     87: . val - the dot product

     89:   Level: intermediate

     91:   Note:
     92:   For complex vectors, `VecDot()` computes
     93: .vb
     94:   val = (x,y) = y^H x,
     95: .ve
     96:   where $y^H$ denotes the conjugate transpose of `y`. Note that this corresponds to the usual "mathematicians" complex
     97:   inner product where the SECOND argument gets the complex conjugate. Since the `BLASdot()` complex conjugates the first
     98:   first argument we call the `BLASdot()` with the arguments reversed.

    100:   Use `VecTDot()` for the indefinite form
    101: .vb
    102:   val = (x,y) = y^T x,
    103: .ve
    104:   where $y^T$ denotes the transpose of `y`.

    106: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDotRealPart()`
    107: @*/
    108: PetscErrorCode VecDot(Vec x, Vec y, PetscScalar *val)
    109: {
    110:   PetscFunctionBegin;
    113:   PetscAssertPointer(val, 3);
    116:   PetscCheckSameTypeAndComm(x, 1, y, 2);
    117:   VecCheckSameSize(x, 1, y, 2);
    118:   VecCheckAssembled(x);
    119:   VecCheckAssembled(y);

    121:   PetscCall(VecLockReadPush(x));
    122:   PetscCall(VecLockReadPush(y));
    123:   PetscCall(PetscLogEventBegin(VEC_Dot, x, y, 0, 0));
    124:   PetscUseTypeMethod(x, dot, y, val);
    125:   PetscCall(PetscLogEventEnd(VEC_Dot, x, y, 0, 0));
    126:   PetscCall(VecLockReadPop(x));
    127:   PetscCall(VecLockReadPop(y));
    128:   PetscFunctionReturn(PETSC_SUCCESS);
    129: }

    131: /*@
    132:   VecDotRealPart - Computes the real part of the vector dot product.

    134:   Collective

    136:   Input Parameters:
    137: + x - first vector
    138: - y - second vector

    140:   Output Parameter:
    141: . val - the real part of the dot product;

    143:   Level: intermediate

    145:   Notes for Users of Complex Numbers:
    146:   See `VecDot()` for more details on the definition of the dot product for complex numbers

    148:   For real numbers this returns the same value as `VecDot()`

    150:   For complex numbers in C^n (that is a vector of n components with a complex number for each component) this is equal to the usual real dot product on the
    151:   the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)

    153:   Developer Notes:
    154:   This is not currently optimized to compute only the real part of the dot product.

    156: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDot()`, `VecDotNorm2()`
    157: @*/
    158: PetscErrorCode VecDotRealPart(Vec x, Vec y, PetscReal *val)
    159: {
    160:   PetscScalar fdot;

    162:   PetscFunctionBegin;
    163:   PetscCall(VecDot(x, y, &fdot));
    164:   *val = PetscRealPart(fdot);
    165:   PetscFunctionReturn(PETSC_SUCCESS);
    166: }

    168: /*@
    169:   VecNorm  - Computes the vector norm.

    171:   Collective

    173:   Input Parameters:
    174: + x    - the vector
    175: - type - the type of the norm requested

    177:   Output Parameter:
    178: . val - the norm

    180:   Level: intermediate

    182:   Notes:
    183:   See `NormType` for descriptions of each norm.

    185:   For complex numbers `NORM_1` will return the traditional 1 norm of the 2 norm of the complex
    186:   numbers; that is the 1 norm of the absolute values of the complex entries. In PETSc 3.6 and
    187:   earlier releases it returned the 1 norm of the 1 norm of the complex entries (what is
    188:   returned by the BLAS routine `asum()`). Both are valid norms but most people expect the former.

    190:   This routine stashes the computed norm value, repeated calls before the vector entries are
    191:   changed are then rapid since the precomputed value is immediately available. Certain vector
    192:   operations such as `VecSet()` store the norms so the value is immediately available and does
    193:   not need to be explicitly computed. `VecScale()` updates any stashed norm values, thus calls
    194:   after `VecScale()` do not need to explicitly recompute the norm.

    196: .seealso: [](ch_vectors), `Vec`, `NormType`, `VecDot()`, `VecTDot()`, `VecDotBegin()`, `VecDotEnd()`, `VecNormAvailable()`,
    197:           `VecNormBegin()`, `VecNormEnd()`, `NormType()`
    198: @*/
    199: PetscErrorCode VecNorm(Vec x, NormType type, PetscReal *val)
    200: {
    201:   PetscBool flg = PETSC_TRUE;

    203:   PetscFunctionBegin;
    204:   PetscCall(VecLockReadPush(x));
    207:   VecCheckAssembled(x);
    209:   PetscAssertPointer(val, 3);

    211:   PetscCall(VecNormAvailable(x, type, &flg, val));
    212:   // check that all MPI processes call this routine together and have same availability
    213:   if (PetscDefined(USE_DEBUG)) {
    214:     PetscMPIInt b0 = (PetscMPIInt)flg, b1[2], b2[2];
    215:     b1[0]          = -b0;
    216:     b1[1]          = b0;
    217:     PetscCallMPI(MPIU_Allreduce(b1, b2, 2, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
    218:     PetscCheck(-b2[0] == b2[1], PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_WRONGSTATE, "Some MPI processes have cached %s norm, others do not. This may happen when some MPI processes call VecGetArray() and some others do not.", NormTypes[type]);
    219:     if (flg) {
    220:       PetscReal b1[2], b2[2];
    221:       b1[0] = -(*val);
    222:       b1[1] = *val;
    223:       PetscCallMPI(MPIU_Allreduce(b1, b2, 2, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)x)));
    224:       PetscCheck((PetscIsNanReal(b2[0]) && PetscIsNanReal(b2[1])) || (-b2[0] == b2[1]), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Difference in cached %s norms: local %g", NormTypes[type], (double)*val);
    225:     }
    226:   }
    227:   if (!flg) {
    228:     PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
    229:     PetscUseTypeMethod(x, norm, type, val);
    230:     PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));

    232:     if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
    233:   }
    234:   PetscCall(VecLockReadPop(x));
    235:   PetscFunctionReturn(PETSC_SUCCESS);
    236: }

    238: /*@
    239:   VecNormAvailable  - Returns the vector norm if it is already known. That is, it has been previously computed and cached in the vector

    241:   Not Collective

    243:   Input Parameters:
    244: + x    - the vector
    245: - type - one of `NORM_1` (sum_i |x[i]|), `NORM_2` sqrt(sum_i (x[i])^2), `NORM_INFINITY` max_i |x[i]|.  Also available
    246:           `NORM_1_AND_2`, which computes both norms and stores them
    247:           in a two element array.

    249:   Output Parameters:
    250: + available - `PETSC_TRUE` if the val returned is valid
    251: - val       - the norm

    253:   Level: intermediate

    255: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
    256:           `VecNormBegin()`, `VecNormEnd()`
    257: @*/
    258: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
    259: {
    260:   PetscFunctionBegin;
    263:   PetscAssertPointer(available, 3);
    264:   PetscAssertPointer(val, 4);

    266:   if (type == NORM_1_AND_2) {
    267:     *available = PETSC_FALSE;
    268:   } else {
    269:     PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
    270:   }
    271:   PetscFunctionReturn(PETSC_SUCCESS);
    272: }

    274: /*@
    275:   VecNormalize - Normalizes a vector by its 2-norm.

    277:   Collective

    279:   Input Parameter:
    280: . x - the vector

    282:   Output Parameter:
    283: . val - the vector norm before normalization. May be `NULL` if the value is not needed.

    285:   Level: intermediate

    287: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
    288: @*/
    289: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
    290: {
    291:   PetscReal norm;

    293:   PetscFunctionBegin;
    296:   PetscCall(VecSetErrorIfLocked(x, 1));
    297:   if (val) PetscAssertPointer(val, 2);
    298:   PetscCall(PetscLogEventBegin(VEC_Normalize, x, 0, 0, 0));
    299:   PetscCall(VecNorm(x, NORM_2, &norm));
    300:   if (norm == 0.0) PetscCall(PetscInfo(x, "Vector of zero norm can not be normalized; Returning only the zero norm\n"));
    301:   else if (PetscIsInfOrNanReal(norm)) PetscCall(PetscInfo(x, "Vector with infinity or NaN norm can not be normalized; Returning only the norm\n"));
    302:   else {
    303:     PetscScalar s = 1.0 / norm;
    304:     PetscCall(VecScale(x, s));
    305:   }
    306:   PetscCall(PetscLogEventEnd(VEC_Normalize, x, 0, 0, 0));
    307:   if (val) *val = norm;
    308:   PetscFunctionReturn(PETSC_SUCCESS);
    309: }

    311: /*@
    312:   VecMax - Determines the vector component with maximum real part and its location.

    314:   Collective

    316:   Input Parameter:
    317: . x - the vector

    319:   Output Parameters:
    320: + p   - the index of `val` (pass `NULL` if you don't want this) in the vector
    321: - val - the maximum component

    323:   Level: intermediate

    325:   Notes:
    326:   Returns the value `PETSC_MIN_REAL` and negative `p` if the vector is of length 0.

    328:   Returns the smallest index with the maximum value

    330:   Developer Note:
    331:   The Nag Fortran compiler does not like the symbol name VecMax

    333: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
    334: @*/
    335: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
    336: {
    337:   PetscFunctionBegin;
    340:   VecCheckAssembled(x);
    341:   if (p) PetscAssertPointer(p, 2);
    342:   PetscAssertPointer(val, 3);
    343:   PetscCall(VecLockReadPush(x));
    344:   PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
    345:   PetscUseTypeMethod(x, max, p, val);
    346:   PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
    347:   PetscCall(VecLockReadPop(x));
    348:   PetscFunctionReturn(PETSC_SUCCESS);
    349: }

    351: /*@
    352:   VecMin - Determines the vector component with minimum real part and its location.

    354:   Collective

    356:   Input Parameter:
    357: . x - the vector

    359:   Output Parameters:
    360: + p   - the index of `val` (pass `NULL` if you don't want this location) in the vector
    361: - val - the minimum component

    363:   Level: intermediate

    365:   Notes:
    366:   Returns the value `PETSC_MAX_REAL` and negative `p` if the vector is of length 0.

    368:   This returns the smallest index with the minimum value

    370:   Developer Note:
    371:   The Nag Fortran compiler does not like the symbol name VecMin

    373: .seealso: [](ch_vectors), `Vec`, `VecMax()`
    374: @*/
    375: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
    376: {
    377:   PetscFunctionBegin;
    380:   VecCheckAssembled(x);
    381:   if (p) PetscAssertPointer(p, 2);
    382:   PetscAssertPointer(val, 3);
    383:   PetscCall(VecLockReadPush(x));
    384:   PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
    385:   PetscUseTypeMethod(x, min, p, val);
    386:   PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
    387:   PetscCall(VecLockReadPop(x));
    388:   PetscFunctionReturn(PETSC_SUCCESS);
    389: }

    391: /*@
    392:   VecTDot - Computes an indefinite vector dot product. That is, this
    393:   routine does NOT use the complex conjugate.

    395:   Collective

    397:   Input Parameters:
    398: + x - first vector
    399: - y - second vector

    401:   Output Parameter:
    402: . val - the dot product

    404:   Level: intermediate

    406:   Notes for Users of Complex Numbers:
    407:   For complex vectors, `VecTDot()` computes the indefinite form
    408: .vb
    409:   val = (x,y) = y^T x,
    410: .ve
    411:   where y^T denotes the transpose of y.

    413:   Use `VecDot()` for the inner product
    414: .vb
    415:   val = (x,y) = y^H x,
    416: .ve
    417:   where y^H denotes the conjugate transpose of y.

    419: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
    420: @*/
    421: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
    422: {
    423:   PetscFunctionBegin;
    426:   PetscAssertPointer(val, 3);
    429:   PetscCheckSameTypeAndComm(x, 1, y, 2);
    430:   VecCheckSameSize(x, 1, y, 2);
    431:   VecCheckAssembled(x);
    432:   VecCheckAssembled(y);

    434:   PetscCall(VecLockReadPush(x));
    435:   PetscCall(VecLockReadPush(y));
    436:   PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
    437:   PetscUseTypeMethod(x, tdot, y, val);
    438:   PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
    439:   PetscCall(VecLockReadPop(x));
    440:   PetscCall(VecLockReadPop(y));
    441:   PetscFunctionReturn(PETSC_SUCCESS);
    442: }

    444: PetscErrorCode VecScaleAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
    445: {
    446:   PetscReal   norms[4];
    447:   PetscBool   flgs[4];
    448:   PetscScalar one = 1.0;

    450:   PetscFunctionBegin;
    453:   VecCheckAssembled(x);
    454:   PetscCall(VecSetErrorIfLocked(x, 1));
    456:   if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);

    458:   /* get current stashed norms */
    459:   for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));

    461:   PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
    462:   VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
    463:   PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));

    465:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
    466:   /* put the scaled stashed norms back into the Vec */
    467:   for (PetscInt i = 0; i < 4; i++) {
    468:     PetscReal ar = PetscAbsScalar(alpha);
    469:     if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
    470:   }
    471:   PetscFunctionReturn(PETSC_SUCCESS);
    472: }

    474: /*@
    475:   VecScale - Scales a vector.

    477:   Logically Collective

    479:   Input Parameters:
    480: + x     - the vector
    481: - alpha - the scalar

    483:   Level: intermediate

    485:   Note:
    486:   For a vector with n components, `VecScale()` computes  x[i] = alpha * x[i], for i=1,...,n.

    488: .seealso: [](ch_vectors), `Vec`, `VecSet()`
    489: @*/
    490: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
    491: {
    492:   PetscFunctionBegin;
    493:   PetscCall(VecScaleAsync_Private(x, alpha, NULL));
    494:   PetscFunctionReturn(PETSC_SUCCESS);
    495: }

    497: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
    498: {
    499:   PetscFunctionBegin;
    502:   VecCheckAssembled(x);
    504:   PetscCall(VecSetErrorIfLocked(x, 1));

    506:   if (alpha == 0) {
    507:     PetscReal norm;
    508:     PetscBool set;

    510:     PetscCall(VecNormAvailable(x, NORM_2, &set, &norm));
    511:     if (set == PETSC_TRUE && norm == 0) PetscFunctionReturn(PETSC_SUCCESS);
    512:   }
    513:   PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
    514:   VecMethodDispatch(x, dctx, VecAsyncFnName(Set), set, (Vec, PetscScalar, PetscDeviceContext), alpha);
    515:   PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
    516:   PetscCall(PetscObjectStateIncrease((PetscObject)x));

    518:   /*  norms can be simply set (if |alpha|*N not too large) */
    519:   {
    520:     PetscReal      val = PetscAbsScalar(alpha);
    521:     const PetscInt N   = x->map->N;

    523:     if (N == 0) {
    524:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0));
    525:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
    526:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
    527:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
    528:     } else if (val > PETSC_MAX_REAL / N) {
    529:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
    530:     } else {
    531:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
    532:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
    533:       val *= PetscSqrtReal((PetscReal)N);
    534:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
    535:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
    536:     }
    537:   }
    538:   PetscFunctionReturn(PETSC_SUCCESS);
    539: }

    541: /*@
    542:   VecSet - Sets all components of a vector to a single scalar value.

    544:   Logically Collective

    546:   Input Parameters:
    547: + x     - the vector
    548: - alpha - the scalar

    550:   Level: beginner

    552:   Notes:
    553:   For a vector of dimension n, `VecSet()` sets x[i] = alpha, for i=1,...,n,
    554:   so that all vector entries then equal the identical
    555:   scalar value, `alpha`.  Use the more general routine
    556:   `VecSetValues()` to set different vector entries.

    558:   You CANNOT call this after you have called `VecSetValues()` but before you call
    559:   `VecAssemblyBegin()`

    561:   If `alpha` is zero and the norm of the vector is known to be zero then this skips the unneeded zeroing process

    563: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
    564: @*/
    565: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
    566: {
    567:   PetscFunctionBegin;
    568:   PetscCall(VecSetAsync_Private(x, alpha, NULL));
    569:   PetscFunctionReturn(PETSC_SUCCESS);
    570: }

    572: PetscErrorCode VecAXPYAsync_Private(Vec y, PetscScalar alpha, Vec x, PetscDeviceContext dctx)
    573: {
    574:   PetscFunctionBegin;
    579:   PetscCheckSameTypeAndComm(x, 3, y, 1);
    580:   VecCheckSameSize(x, 3, y, 1);
    581:   VecCheckAssembled(x);
    582:   VecCheckAssembled(y);
    584:   if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
    585:   PetscCall(VecSetErrorIfLocked(y, 1));
    586:   if (x == y) {
    587:     PetscCall(VecScale(y, alpha + 1.0));
    588:     PetscFunctionReturn(PETSC_SUCCESS);
    589:   }
    590:   PetscCall(VecLockReadPush(x));
    591:   PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
    592:   VecMethodDispatch(y, dctx, VecAsyncFnName(AXPY), axpy, (Vec, PetscScalar, Vec, PetscDeviceContext), alpha, x);
    593:   PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
    594:   PetscCall(VecLockReadPop(x));
    595:   PetscCall(PetscObjectStateIncrease((PetscObject)y));
    596:   PetscFunctionReturn(PETSC_SUCCESS);
    597: }
    598: /*@
    599:   VecAXPY - Computes `y = alpha x + y`.

    601:   Logically Collective

    603:   Input Parameters:
    604: + alpha - the scalar
    605: . x     - vector scale by `alpha`
    606: - y     - vector accumulated into

    608:   Output Parameter:
    609: . y - output vector

    611:   Level: intermediate

    613:   Notes:
    614:   This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
    615: .vb
    616:     VecAXPY(y,alpha,x)                   y = alpha x           +      y
    617:     VecAYPX(y,beta,x)                    y =       x           + beta y
    618:     VecAXPBY(y,alpha,beta,x)             y = alpha x           + beta y
    619:     VecWAXPY(w,alpha,x,y)                w = alpha x           +      y
    620:     VecAXPBYPCZ(z,alpha,beta,gamma,x,y)  z = alpha x           + beta y + gamma z
    621:     VecMAXPY(y,nv,alpha[],x[])           y = sum alpha[i] x[i] +      y
    622: .ve

    624: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
    625: @*/
    626: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
    627: {
    628:   PetscFunctionBegin;
    629:   PetscCall(VecAXPYAsync_Private(y, alpha, x, NULL));
    630:   PetscFunctionReturn(PETSC_SUCCESS);
    631: }

    633: PetscErrorCode VecAYPXAsync_Private(Vec y, PetscScalar beta, Vec x, PetscDeviceContext dctx)
    634: {
    635:   PetscFunctionBegin;
    640:   PetscCheckSameTypeAndComm(x, 3, y, 1);
    641:   VecCheckSameSize(x, 1, y, 3);
    642:   VecCheckAssembled(x);
    643:   VecCheckAssembled(y);
    645:   PetscCall(VecSetErrorIfLocked(y, 1));
    646:   if (x == y) {
    647:     PetscCall(VecScale(y, beta + 1.0));
    648:     PetscFunctionReturn(PETSC_SUCCESS);
    649:   }
    650:   PetscCall(VecLockReadPush(x));
    651:   if (beta == (PetscScalar)0.0) {
    652:     PetscCall(VecCopy(x, y));
    653:   } else {
    654:     PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
    655:     VecMethodDispatch(y, dctx, VecAsyncFnName(AYPX), aypx, (Vec, PetscScalar, Vec, PetscDeviceContext), beta, x);
    656:     PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
    657:     PetscCall(PetscObjectStateIncrease((PetscObject)y));
    658:   }
    659:   PetscCall(VecLockReadPop(x));
    660:   PetscFunctionReturn(PETSC_SUCCESS);
    661: }

    663: /*@
    664:   VecAYPX - Computes `y = x + beta y`.

    666:   Logically Collective

    668:   Input Parameters:
    669: + beta - the scalar
    670: . x    - the unscaled vector
    671: - y    - the vector to be scaled

    673:   Output Parameter:
    674: . y - output vector

    676:   Level: intermediate

    678:   Developer Notes:
    679:   The implementation is optimized for `beta` of -1.0, 0.0, and 1.0

    681: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
    682: @*/
    683: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
    684: {
    685:   PetscFunctionBegin;
    686:   PetscCall(VecAYPXAsync_Private(y, beta, x, NULL));
    687:   PetscFunctionReturn(PETSC_SUCCESS);
    688: }

    690: PetscErrorCode VecAXPBYAsync_Private(Vec y, PetscScalar alpha, PetscScalar beta, Vec x, PetscDeviceContext dctx)
    691: {
    692:   PetscFunctionBegin;
    697:   PetscCheckSameTypeAndComm(x, 4, y, 1);
    698:   VecCheckSameSize(y, 1, x, 4);
    699:   VecCheckAssembled(x);
    700:   VecCheckAssembled(y);
    703:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
    704:   if (x == y) {
    705:     PetscCall(VecScale(y, alpha + beta));
    706:     PetscFunctionReturn(PETSC_SUCCESS);
    707:   }

    709:   PetscCall(VecSetErrorIfLocked(y, 1));
    710:   PetscCall(VecLockReadPush(x));
    711:   PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
    712:   VecMethodDispatch(y, dctx, VecAsyncFnName(AXPBY), axpby, (Vec, PetscScalar, PetscScalar, Vec, PetscDeviceContext), alpha, beta, x);
    713:   PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
    714:   PetscCall(PetscObjectStateIncrease((PetscObject)y));
    715:   PetscCall(VecLockReadPop(x));
    716:   PetscFunctionReturn(PETSC_SUCCESS);
    717: }

    719: /*@
    720:   VecAXPBY - Computes `y = alpha x + beta y`.

    722:   Logically Collective

    724:   Input Parameters:
    725: + alpha - first scalar
    726: . beta  - second scalar
    727: . x     - the first scaled vector
    728: - y     - the second scaled vector

    730:   Output Parameter:
    731: . y - output vector

    733:   Level: intermediate

    735:   Developer Notes:
    736:   The implementation is optimized for `alpha` and/or `beta` values of 0.0 and 1.0

    738: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
    739: @*/
    740: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
    741: {
    742:   PetscFunctionBegin;
    743:   PetscCall(VecAXPBYAsync_Private(y, alpha, beta, x, NULL));
    744:   PetscFunctionReturn(PETSC_SUCCESS);
    745: }

    747: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
    748: {
    749:   PetscFunctionBegin;
    756:   PetscCheckSameTypeAndComm(x, 5, y, 6);
    757:   PetscCheckSameTypeAndComm(x, 5, z, 1);
    758:   VecCheckSameSize(x, 5, y, 6);
    759:   VecCheckSameSize(x, 5, z, 1);
    760:   PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
    761:   PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
    762:   VecCheckAssembled(x);
    763:   VecCheckAssembled(y);
    764:   VecCheckAssembled(z);
    768:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);

    770:   PetscCall(VecSetErrorIfLocked(z, 1));
    771:   PetscCall(VecLockReadPush(x));
    772:   PetscCall(VecLockReadPush(y));
    773:   PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
    774:   VecMethodDispatch(z, dctx, VecAsyncFnName(AXPBYPCZ), axpbypcz, (Vec, PetscScalar, PetscScalar, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, beta, gamma, x, y);
    775:   PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
    776:   PetscCall(PetscObjectStateIncrease((PetscObject)z));
    777:   PetscCall(VecLockReadPop(x));
    778:   PetscCall(VecLockReadPop(y));
    779:   PetscFunctionReturn(PETSC_SUCCESS);
    780: }
    781: /*@
    782:   VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`

    784:   Logically Collective

    786:   Input Parameters:
    787: + alpha - first scalar
    788: . beta  - second scalar
    789: . gamma - third scalar
    790: . x     - first vector
    791: . y     - second vector
    792: - z     - third vector

    794:   Output Parameter:
    795: . z - output vector

    797:   Level: intermediate

    799:   Note:
    800:   `x`, `y` and `z` must be different vectors

    802:   Developer Notes:
    803:   The implementation is optimized for `alpha` of 1.0 and `gamma` of 1.0 or 0.0

    805: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
    806: @*/
    807: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
    808: {
    809:   PetscFunctionBegin;
    810:   PetscCall(VecAXPBYPCZAsync_Private(z, alpha, beta, gamma, x, y, NULL));
    811:   PetscFunctionReturn(PETSC_SUCCESS);
    812: }

    814: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
    815: {
    816:   PetscFunctionBegin;
    823:   PetscCheckSameTypeAndComm(x, 3, y, 4);
    824:   PetscCheckSameTypeAndComm(y, 4, w, 1);
    825:   VecCheckSameSize(x, 3, y, 4);
    826:   VecCheckSameSize(x, 3, w, 1);
    827:   PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
    828:   PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
    829:   VecCheckAssembled(x);
    830:   VecCheckAssembled(y);
    832:   PetscCall(VecSetErrorIfLocked(w, 1));

    834:   PetscCall(VecLockReadPush(x));
    835:   PetscCall(VecLockReadPush(y));
    836:   if (alpha == (PetscScalar)0.0) {
    837:     PetscCall(VecCopyAsync_Private(y, w, dctx));
    838:   } else {
    839:     PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
    840:     VecMethodDispatch(w, dctx, VecAsyncFnName(WAXPY), waxpy, (Vec, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, x, y);
    841:     PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
    842:     PetscCall(PetscObjectStateIncrease((PetscObject)w));
    843:   }
    844:   PetscCall(VecLockReadPop(x));
    845:   PetscCall(VecLockReadPop(y));
    846:   PetscFunctionReturn(PETSC_SUCCESS);
    847: }

    849: /*@
    850:   VecWAXPY - Computes `w = alpha x + y`.

    852:   Logically Collective

    854:   Input Parameters:
    855: + alpha - the scalar
    856: . x     - first vector, multiplied by `alpha`
    857: - y     - second vector

    859:   Output Parameter:
    860: . w - the result

    862:   Level: intermediate

    864:   Note:
    865:   `w` cannot be either `x` or `y`, but `x` and `y` can be the same

    867:   Developer Notes:
    868:   The implementation is optimized for alpha of -1.0, 0.0, and 1.0

    870: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
    871: @*/
    872: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
    873: {
    874:   PetscFunctionBegin;
    875:   PetscCall(VecWAXPYAsync_Private(w, alpha, x, y, NULL));
    876:   PetscFunctionReturn(PETSC_SUCCESS);
    877: }

    879: /*@
    880:   VecSetValues - Inserts or adds values into certain locations of a vector.

    882:   Not Collective

    884:   Input Parameters:
    885: + x    - vector to insert in
    886: . ni   - number of elements to add
    887: . ix   - indices where to add
    888: . y    - array of values. Pass `NULL` to set all zeroes.
    889: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries

    891:   Level: beginner

    893:   Notes:
    894: .vb
    895:    `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
    896: .ve

    898:   Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
    899:   options cannot be mixed without intervening calls to the assembly
    900:   routines.

    902:   These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
    903:   MUST be called after all calls to `VecSetValues()` have been completed.

    905:   VecSetValues() uses 0-based indices in Fortran as well as in C.

    907:   If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
    908:   negative indices may be passed in ix. These rows are
    909:   simply ignored. This allows easily inserting element load matrices
    910:   with homogeneous Dirichlet boundary conditions that you don't want represented
    911:   in the vector.

    913:   Fortran Note:
    914:   If any of `ix` and `y` are scalars pass them using, for example,
    915: .vb
    916:   call VecSetValues(mat, one, [ix], [y], INSERT_VALUES, ierr)
    917: .ve

    919: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
    920:           `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
    921: @*/
    922: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
    923: {
    924:   PetscFunctionBeginHot;
    926:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
    927:   PetscAssertPointer(ix, 3);
    928:   if (y) PetscAssertPointer(y, 4);

    931:   PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
    932:   PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
    933:   PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
    934:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
    935:   PetscFunctionReturn(PETSC_SUCCESS);
    936: }

    938: /*@
    939:   VecGetValues - Gets values from certain locations of a vector. Currently
    940:   can only get values on the same processor on which they are owned

    942:   Not Collective

    944:   Input Parameters:
    945: + x  - vector to get values from
    946: . ni - number of elements to get
    947: - ix - indices where to get them from (in global 1d numbering)

    949:   Output Parameter:
    950: . y - array of values, must be passed in with a length of `ni`

    952:   Level: beginner

    954:   Notes:
    955:   The user provides the allocated array y; it is NOT allocated in this routine

    957:   `VecGetValues()` gets y[i] = x[ix[i]], for i=0,...,ni-1.

    959:   `VecAssemblyBegin()` and `VecAssemblyEnd()`  MUST be called before calling this if `VecSetValues()` or related routine has been called

    961:   VecGetValues() uses 0-based indices in Fortran as well as in C.

    963:   If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
    964:   negative indices may be passed in ix. These rows are
    965:   simply ignored.

    967: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
    968: @*/
    969: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
    970: {
    971:   PetscFunctionBegin;
    973:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
    974:   PetscAssertPointer(ix, 3);
    975:   PetscAssertPointer(y, 4);
    977:   VecCheckAssembled(x);
    978:   PetscUseTypeMethod(x, getvalues, ni, ix, y);
    979:   PetscFunctionReturn(PETSC_SUCCESS);
    980: }

    982: /*@
    983:   VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.

    985:   Not Collective

    987:   Input Parameters:
    988: + x    - vector to insert in
    989: . ni   - number of blocks to add
    990: . ix   - indices where to add in block count, rather than element count
    991: . y    - array of values. Pass `NULL` to set all zeroes.
    992: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries

    994:   Level: intermediate

    996:   Notes:
    997:   `VecSetValuesBlocked()` sets x[bs*ix[i]+j] = y[bs*i+j],
    998:   for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().

    1000:   Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
    1001:   options cannot be mixed without intervening calls to the assembly
    1002:   routines.

    1004:   These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
    1005:   MUST be called after all calls to `VecSetValuesBlocked()` have been completed.

    1007:   `VecSetValuesBlocked()` uses 0-based indices in Fortran as well as in C.

    1009:   Negative indices may be passed in ix, these rows are
    1010:   simply ignored. This allows easily inserting element load matrices
    1011:   with homogeneous Dirichlet boundary conditions that you don't want represented
    1012:   in the vector.

    1014:   Fortran Note:
    1015:   If any of `ix` and `y` are scalars pass them using, for example,
    1016: .vb
    1017:   call VecSetValuesBlocked(mat, one, [ix], [y], INSERT_VALUES, ierr)
    1018: .ve

    1020: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
    1021:           `VecSetValues()`
    1022: @*/
    1023: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
    1024: {
    1025:   PetscFunctionBeginHot;
    1027:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
    1028:   PetscAssertPointer(ix, 3);
    1029:   if (y) PetscAssertPointer(y, 4);

    1032:   PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
    1033:   PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
    1034:   PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
    1035:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
    1036:   PetscFunctionReturn(PETSC_SUCCESS);
    1037: }

    1039: /*@
    1040:   VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
    1041:   using a local ordering of the nodes.

    1043:   Not Collective

    1045:   Input Parameters:
    1046: + x    - vector to insert in
    1047: . ni   - number of elements to add
    1048: . ix   - indices where to add
    1049: . y    - array of values. Pass `NULL` to set all zeroes.
    1050: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries

    1052:   Level: intermediate

    1054:   Notes:
    1055:   `VecSetValuesLocal()` sets x[ix[i]] = y[i], for i=0,...,ni-1.

    1057:   Calls to `VecSetValuesLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
    1058:   options cannot be mixed without intervening calls to the assembly
    1059:   routines.

    1061:   These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
    1062:   MUST be called after all calls to `VecSetValuesLocal()` have been completed.

    1064:   `VecSetValuesLocal()` uses 0-based indices in Fortran as well as in C.

    1066:   Fortran Note:
    1067:   If any of `ix` and `y` are scalars pass them using, for example,
    1068: .vb
    1069:   call VecSetValuesLocal(mat, one, [ix], [y], INSERT_VALUES, ierr)
    1070: .ve

    1072: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
    1073:           `VecSetValuesBlockedLocal()`
    1074: @*/
    1075: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
    1076: {
    1077:   PetscInt lixp[128], *lix = lixp;

    1079:   PetscFunctionBeginHot;
    1081:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
    1082:   PetscAssertPointer(ix, 3);
    1083:   if (y) PetscAssertPointer(y, 4);

    1086:   PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
    1087:   if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
    1088:   if (x->map->mapping) {
    1089:     if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
    1090:     PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
    1091:     PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
    1092:     if (ni > 128) PetscCall(PetscFree(lix));
    1093:   } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
    1094:   PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
    1095:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
    1096:   PetscFunctionReturn(PETSC_SUCCESS);
    1097: }

    1099: /*@
    1100:   VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
    1101:   using a local ordering of the nodes.

    1103:   Not Collective

    1105:   Input Parameters:
    1106: + x    - vector to insert in
    1107: . ni   - number of blocks to add
    1108: . ix   - indices where to add in block count, not element count
    1109: . y    - array of values. Pass `NULL` to set all zeroes.
    1110: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries

    1112:   Level: intermediate

    1114:   Notes:
    1115:   `VecSetValuesBlockedLocal()` sets x[bs*ix[i]+j] = y[bs*i+j],
    1116:   for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with `VecSetBlockSize()`.

    1118:   Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
    1119:   options cannot be mixed without intervening calls to the assembly
    1120:   routines.

    1122:   These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
    1123:   MUST be called after all calls to `VecSetValuesBlockedLocal()` have been completed.

    1125:   `VecSetValuesBlockedLocal()` uses 0-based indices in Fortran as well as in C.

    1127:   Fortran Note:
    1128:   If any of `ix` and `y` are scalars pass them using, for example,
    1129: .vb
    1130:   call VecSetValuesBlockedLocal(mat, one, [ix], [y], INSERT_VALUES, ierr)
    1131: .ve

    1133: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
    1134:           `VecSetLocalToGlobalMapping()`
    1135: @*/
    1136: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
    1137: {
    1138:   PetscInt lixp[128], *lix = lixp;

    1140:   PetscFunctionBeginHot;
    1142:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
    1143:   PetscAssertPointer(ix, 3);
    1144:   if (y) PetscAssertPointer(y, 4);
    1146:   PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
    1147:   if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
    1148:   if (x->map->mapping) {
    1149:     if (ni > (PetscInt)PETSC_STATIC_ARRAY_LENGTH(lixp)) PetscCall(PetscMalloc1(ni, &lix));
    1150:     PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
    1151:     PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
    1152:     if (ni > (PetscInt)PETSC_STATIC_ARRAY_LENGTH(lixp)) PetscCall(PetscFree(lix));
    1153:   } else {
    1154:     PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
    1155:   }
    1156:   PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
    1157:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
    1158:   PetscFunctionReturn(PETSC_SUCCESS);
    1159: }

    1161: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
    1162: {
    1163:   PetscFunctionBegin;
    1166:   VecCheckAssembled(x);
    1168:   if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
    1169:   PetscAssertPointer(y, 3);
    1170:   for (PetscInt i = 0; i < nv; ++i) {
    1173:     PetscCheckSameTypeAndComm(x, 1, y[i], 3);
    1174:     VecCheckSameSize(x, 1, y[i], 3);
    1175:     VecCheckAssembled(y[i]);
    1176:     PetscCall(VecLockReadPush(y[i]));
    1177:   }
    1178:   PetscAssertPointer(result, 4);

    1181:   PetscCall(VecLockReadPush(x));
    1182:   PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
    1183:   PetscCall((*mxdot)(x, nv, y, result));
    1184:   PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
    1185:   PetscCall(VecLockReadPop(x));
    1186:   for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
    1187:   PetscFunctionReturn(PETSC_SUCCESS);
    1188: }

    1190: /*@
    1191:   VecMTDot - Computes indefinite vector multiple dot products.
    1192:   That is, it does NOT use the complex conjugate.

    1194:   Collective

    1196:   Input Parameters:
    1197: + x  - one vector
    1198: . nv - number of vectors
    1199: - y  - array of vectors.  Note that vectors are pointers

    1201:   Output Parameter:
    1202: . val - array of the dot products

    1204:   Level: intermediate

    1206:   Notes for Users of Complex Numbers:
    1207:   For complex vectors, `VecMTDot()` computes the indefinite form
    1208: .vb
    1209:   val = (x,y) = y^T x,
    1210: .ve
    1211:   where y^T denotes the transpose of y.

    1213:   Use `VecMDot()` for the inner product
    1214: .vb
    1215:   val = (x,y) = y^H x,
    1216: .ve
    1217:   where y^H denotes the conjugate transpose of y.

    1219: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
    1220: @*/
    1221: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
    1222: {
    1223:   PetscFunctionBegin;
    1225:   PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
    1226:   PetscFunctionReturn(PETSC_SUCCESS);
    1227: }

    1229: /*@
    1230:   VecMDot - Computes multiple vector dot products.

    1232:   Collective

    1234:   Input Parameters:
    1235: + x  - one vector
    1236: . nv - number of vectors
    1237: - y  - array of vectors.

    1239:   Output Parameter:
    1240: . val - array of the dot products (does not allocate the array)

    1242:   Level: intermediate

    1244:   Notes for Users of Complex Numbers:
    1245:   For complex vectors, `VecMDot()` computes
    1246: .vb
    1247:   val = (x,y) = y^H x,
    1248: .ve
    1249:   where y^H denotes the conjugate transpose of y.

    1251:   Use `VecMTDot()` for the indefinite form
    1252: .vb
    1253:   val = (x,y) = y^T x,
    1254: .ve
    1255:   where y^T denotes the transpose of y.

    1257:   Note:
    1258:   The implementation may use BLAS 2 operations when the vectors `y` have been obtained with `VecDuplicateVecs()`

    1260: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`, `VecDuplicateVecs()`
    1261: @*/
    1262: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
    1263: {
    1264:   PetscFunctionBegin;
    1266:   PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
    1267:   PetscFunctionReturn(PETSC_SUCCESS);
    1268: }

    1270: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
    1271: {
    1272:   PetscFunctionBegin;
    1274:   VecCheckAssembled(y);
    1276:   PetscCall(VecSetErrorIfLocked(y, 1));
    1277:   PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
    1278:   if (nv) {
    1279:     PetscInt zeros = 0;

    1281:     PetscAssertPointer(alpha, 3);
    1282:     PetscAssertPointer(x, 4);
    1283:     for (PetscInt i = 0; i < nv; ++i) {
    1287:       PetscCheckSameTypeAndComm(y, 1, x[i], 4);
    1288:       VecCheckSameSize(y, 1, x[i], 4);
    1289:       PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
    1290:       VecCheckAssembled(x[i]);
    1291:       PetscCall(VecLockReadPush(x[i]));
    1292:       zeros += alpha[i] == (PetscScalar)0.0;
    1293:     }

    1295:     if (zeros < nv) {
    1296:       PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
    1297:       VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
    1298:       PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
    1299:       PetscCall(PetscObjectStateIncrease((PetscObject)y));
    1300:     }

    1302:     for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
    1303:   }
    1304:   PetscFunctionReturn(PETSC_SUCCESS);
    1305: }

    1307: /*@
    1308:   VecMAXPY - Computes `y = y + sum alpha[i] x[i]`

    1310:   Logically Collective

    1312:   Input Parameters:
    1313: + nv    - number of scalars and `x` vectors
    1314: . alpha - array of scalars
    1315: . y     - one vector
    1316: - x     - array of vectors

    1318:   Level: intermediate

    1320:   Notes:
    1321:   `y` cannot be any of the `x` vectors

    1323:   The implementation may use BLAS 2 operations when the vectors `y` have been obtained with `VecDuplicateVecs()`

    1325: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`, `VecDuplicateVecs()`
    1326: @*/
    1327: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
    1328: {
    1329:   PetscFunctionBegin;
    1330:   PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
    1331:   PetscFunctionReturn(PETSC_SUCCESS);
    1332: }

    1334: /*@
    1335:   VecMAXPBY - Computes `y = beta y + sum alpha[i] x[i]`

    1337:   Logically Collective

    1339:   Input Parameters:
    1340: + nv    - number of scalars and `x` vectors
    1341: . alpha - array of scalars
    1342: . beta  - scalar
    1343: . y     - one vector
    1344: - x     - array of vectors

    1346:   Level: intermediate

    1348:   Note:
    1349:   `y` cannot be any of the `x` vectors.

    1351:   Developer Notes:
    1352:   This is a convenience routine, but implementations might be able to optimize it, for example, when `beta` is zero.

    1354: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
    1355: @*/
    1356: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
    1357: {
    1358:   PetscFunctionBegin;
    1360:   VecCheckAssembled(y);
    1362:   PetscCall(VecSetErrorIfLocked(y, 1));
    1363:   PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);

    1366:   if (y->ops->maxpby) {
    1367:     PetscInt zeros = 0;

    1369:     if (nv) {
    1370:       PetscAssertPointer(alpha, 3);
    1371:       PetscAssertPointer(x, 5);
    1372:     }

    1374:     for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
    1378:       PetscCheckSameTypeAndComm(y, 1, x[i], 5);
    1379:       VecCheckSameSize(y, 1, x[i], 5);
    1380:       PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
    1381:       VecCheckAssembled(x[i]);
    1382:       PetscCall(VecLockReadPush(x[i]));
    1383:       zeros += alpha[i] == (PetscScalar)0.0;
    1384:     }

    1386:     if (zeros < nv) { // has nonzero alpha
    1387:       PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
    1388:       PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
    1389:       PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
    1390:       PetscCall(PetscObjectStateIncrease((PetscObject)y));
    1391:     } else {
    1392:       PetscCall(VecScale(y, beta));
    1393:     }

    1395:     for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
    1396:   } else { // no maxpby
    1397:     if (beta == 0.0) PetscCall(VecSet(y, 0.0));
    1398:     else PetscCall(VecScale(y, beta));
    1399:     PetscCall(VecMAXPY(y, nv, alpha, x));
    1400:   }
    1401:   PetscFunctionReturn(PETSC_SUCCESS);
    1402: }

    1404: /*@
    1405:   VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
    1406:   in the order they appear in the array. The concatenated vector resides on the same
    1407:   communicator and is the same type as the source vectors.

    1409:   Collective

    1411:   Input Parameters:
    1412: + nx - number of vectors to be concatenated
    1413: - X  - array containing the vectors to be concatenated in the order of concatenation

    1415:   Output Parameters:
    1416: + Y    - concatenated vector
    1417: - x_is - array of index sets corresponding to the concatenated components of `Y` (pass `NULL` if not needed)

    1419:   Level: advanced

    1421:   Notes:
    1422:   Concatenation is similar to the functionality of a `VECNEST` object; they both represent combination of
    1423:   different vector spaces. However, concatenated vectors do not store any information about their
    1424:   sub-vectors and own their own data. Consequently, this function provides index sets to enable the
    1425:   manipulation of data in the concatenated vector that corresponds to the original components at creation.

    1427:   This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
    1428:   has to operate on combined vector spaces and cannot utilize `VECNEST` objects due to incompatibility with
    1429:   bound projections.

    1431: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
    1432: @*/
    1433: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
    1434: {
    1435:   MPI_Comm comm;
    1436:   VecType  vec_type;
    1437:   Vec      Ytmp, Xtmp;
    1438:   IS      *is_tmp;
    1439:   PetscInt i, shift = 0, Xnl, Xng, Xbegin;

    1441:   PetscFunctionBegin;
    1445:   PetscAssertPointer(Y, 3);

    1447:   if ((*X)->ops->concatenate) {
    1448:     /* use the dedicated concatenation function if available */
    1449:     PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
    1450:   } else {
    1451:     /* loop over vectors and start creating IS */
    1452:     comm = PetscObjectComm((PetscObject)*X);
    1453:     PetscCall(VecGetType(*X, &vec_type));
    1454:     PetscCall(PetscMalloc1(nx, &is_tmp));
    1455:     for (i = 0; i < nx; i++) {
    1456:       PetscCall(VecGetSize(X[i], &Xng));
    1457:       PetscCall(VecGetLocalSize(X[i], &Xnl));
    1458:       PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
    1459:       PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
    1460:       shift += Xng;
    1461:     }
    1462:     /* create the concatenated vector */
    1463:     PetscCall(VecCreate(comm, &Ytmp));
    1464:     PetscCall(VecSetType(Ytmp, vec_type));
    1465:     PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
    1466:     PetscCall(VecSetUp(Ytmp));
    1467:     /* copy data from X array to Y and return */
    1468:     for (i = 0; i < nx; i++) {
    1469:       PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
    1470:       PetscCall(VecCopy(X[i], Xtmp));
    1471:       PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
    1472:     }
    1473:     *Y = Ytmp;
    1474:     if (x_is) {
    1475:       *x_is = is_tmp;
    1476:     } else {
    1477:       for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
    1478:       PetscCall(PetscFree(is_tmp));
    1479:     }
    1480:   }
    1481:   PetscFunctionReturn(PETSC_SUCCESS);
    1482: }

    1484: /* A helper function for VecGetSubVector to check if we can implement it with no-copy (i.e. the subvector shares
    1485:    memory with the original vector), and the block size of the subvector.

    1487:     Input Parameters:
    1488: +   X - the original vector
    1489: -   is - the index set of the subvector

    1491:     Output Parameters:
    1492: +   contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
    1493: .   start  - start of contiguous block, as an offset from the start of the ownership range of the original vector
    1494: -   blocksize - the block size of the subvector

    1496: */
    1497: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
    1498: {
    1499:   PetscInt  gstart, gend, lstart;
    1500:   PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
    1501:   PetscInt  n, N, ibs, vbs, bs = 1;

    1503:   PetscFunctionBegin;
    1504:   PetscCall(ISGetLocalSize(is, &n));
    1505:   PetscCall(ISGetSize(is, &N));
    1506:   PetscCall(ISGetBlockSize(is, &ibs));
    1507:   PetscCall(VecGetBlockSize(X, &vbs));
    1508:   PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
    1509:   PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
    1510:   /* block size is given by IS if ibs > 1; otherwise, check the vector */
    1511:   if (ibs > 1) {
    1512:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
    1513:     bs = ibs;
    1514:   } else {
    1515:     if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
    1516:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
    1517:     if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
    1518:   }

    1520:   *contig    = red[0];
    1521:   *start     = lstart;
    1522:   *blocksize = bs;
    1523:   PetscFunctionReturn(PETSC_SUCCESS);
    1524: }

    1526: /* A helper function for VecGetSubVector, to be used when we have to build a standalone subvector through VecScatter

    1528:     Input Parameters:
    1529: +   X - the original vector
    1530: .   is - the index set of the subvector
    1531: -   bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()

    1533:     Output Parameter:
    1534: .   Z  - the subvector, which will compose the VecScatter context on output
    1535: */
    1536: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
    1537: {
    1538:   PetscInt   n, N;
    1539:   VecScatter vscat;
    1540:   Vec        Y;

    1542:   PetscFunctionBegin;
    1543:   PetscCall(ISGetLocalSize(is, &n));
    1544:   PetscCall(ISGetSize(is, &N));
    1545:   PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
    1546:   PetscCall(VecSetSizes(Y, n, N));
    1547:   PetscCall(VecSetBlockSize(Y, bs));
    1548:   PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
    1549:   PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
    1550:   PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
    1551:   PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
    1552:   PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
    1553:   PetscCall(VecScatterDestroy(&vscat));
    1554:   *Z = Y;
    1555:   PetscFunctionReturn(PETSC_SUCCESS);
    1556: }

    1558: /*@
    1559:   VecGetSubVector - Gets a vector representing part of another vector

    1561:   Collective

    1563:   Input Parameters:
    1564: + X  - vector from which to extract a subvector
    1565: - is - index set representing portion of `X` to extract

    1567:   Output Parameter:
    1568: . Y - subvector corresponding to `is`

    1570:   Level: advanced

    1572:   Notes:
    1573:   The subvector `Y` should be returned with `VecRestoreSubVector()`.
    1574:   `X` and `is` must be defined on the same communicator

    1576:   Changes to the subvector will be reflected in the `X` vector on the call to `VecRestoreSubVector()`.

    1578:   This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
    1579:   modifying the subvector.  Other non-overlapping subvectors can still be obtained from `X` using this function.

    1581:   The resulting subvector inherits the block size from `is` if greater than one. Otherwise, the block size is guessed from the block size of the original `X`.

    1583: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
    1584: @*/
    1585: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
    1586: {
    1587:   Vec Z;

    1589:   PetscFunctionBegin;
    1592:   PetscCheckSameComm(X, 1, is, 2);
    1593:   PetscAssertPointer(Y, 3);
    1594:   if (X->ops->getsubvector) {
    1595:     PetscUseTypeMethod(X, getsubvector, is, &Z);
    1596:   } else { /* Default implementation currently does no caching */
    1597:     PetscBool contig;
    1598:     PetscInt  n, N, start, bs;

    1600:     PetscCall(ISGetLocalSize(is, &n));
    1601:     PetscCall(ISGetSize(is, &N));
    1602:     PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
    1603:     if (contig) { /* We can do a no-copy implementation */
    1604:       const PetscScalar *x;
    1605:       PetscInt           state = 0;
    1606:       PetscBool          isstd, iscuda, iship;

    1608:       PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
    1609:       PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
    1610:       PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
    1611:       if (iscuda) {
    1612: #if defined(PETSC_HAVE_CUDA)
    1613:         const PetscScalar *x_d;
    1614:         PetscMPIInt        size;
    1615:         PetscOffloadMask   flg;

    1617:         PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
    1618:         PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
    1619:         PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
    1620:         if (x) x += start;
    1621:         if (x_d) x_d += start;
    1622:         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
    1623:         if (size == 1) {
    1624:           PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
    1625:         } else {
    1626:           PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
    1627:         }
    1628:         Z->offloadmask = flg;
    1629: #endif
    1630:       } else if (iship) {
    1631: #if defined(PETSC_HAVE_HIP)
    1632:         const PetscScalar *x_d;
    1633:         PetscMPIInt        size;
    1634:         PetscOffloadMask   flg;

    1636:         PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
    1637:         PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
    1638:         PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
    1639:         if (x) x += start;
    1640:         if (x_d) x_d += start;
    1641:         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
    1642:         if (size == 1) {
    1643:           PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
    1644:         } else {
    1645:           PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
    1646:         }
    1647:         Z->offloadmask = flg;
    1648: #endif
    1649:       } else if (isstd) {
    1650:         PetscMPIInt size;

    1652:         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
    1653:         PetscCall(VecGetArrayRead(X, &x));
    1654:         if (x) x += start;
    1655:         if (size == 1) {
    1656:           PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
    1657:         } else {
    1658:           PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
    1659:         }
    1660:         PetscCall(VecRestoreArrayRead(X, &x));
    1661:       } else { /* default implementation: use place array */
    1662:         PetscCall(VecGetArrayRead(X, &x));
    1663:         PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
    1664:         PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
    1665:         PetscCall(VecSetSizes(Z, n, N));
    1666:         PetscCall(VecSetBlockSize(Z, bs));
    1667:         PetscCall(VecPlaceArray(Z, PetscSafePointerPlusOffset(x, start)));
    1668:         PetscCall(VecRestoreArrayRead(X, &x));
    1669:       }

    1671:       /* this is relevant only in debug mode */
    1672:       PetscCall(VecLockGet(X, &state));
    1673:       if (state) PetscCall(VecLockReadPush(Z));
    1674:       Z->ops->placearray   = NULL;
    1675:       Z->ops->replacearray = NULL;
    1676:     } else { /* Have to create a scatter and do a copy */
    1677:       PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
    1678:     }
    1679:   }
    1680:   /* Record the state when the subvector was gotten so we know whether its values need to be put back */
    1681:   if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
    1682:   PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
    1683:   *Y = Z;
    1684:   PetscFunctionReturn(PETSC_SUCCESS);
    1685: }

    1687: /*@
    1688:   VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`

    1690:   Collective

    1692:   Input Parameters:
    1693: + X  - vector from which subvector was obtained
    1694: . is - index set representing the subset of `X`
    1695: - Y  - subvector being restored

    1697:   Level: advanced

    1699: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
    1700: @*/
    1701: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
    1702: {
    1703:   PETSC_UNUSED PetscObjectState dummystate = 0;
    1704:   PetscBool                     unchanged;

    1706:   PetscFunctionBegin;
    1709:   PetscCheckSameComm(X, 1, is, 2);
    1710:   PetscAssertPointer(Y, 3);

    1713:   if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
    1714:   else {
    1715:     PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
    1716:     if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
    1717:       VecScatter scatter;
    1718:       PetscInt   state;

    1720:       PetscCall(VecLockGet(X, &state));
    1721:       PetscCheck(state == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec X is locked for read-only or read/write access");

    1723:       PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
    1724:       if (scatter) {
    1725:         PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
    1726:         PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
    1727:       } else {
    1728:         PetscBool iscuda, iship;
    1729:         PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
    1730:         PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));

    1732:         if (iscuda) {
    1733: #if defined(PETSC_HAVE_CUDA)
    1734:           PetscOffloadMask ymask = (*Y)->offloadmask;

    1736:           /* The offloadmask of X dictates where to move memory
    1737:               If X GPU data is valid, then move Y data on GPU if needed
    1738:               Otherwise, move back to the CPU */
    1739:           switch (X->offloadmask) {
    1740:           case PETSC_OFFLOAD_BOTH:
    1741:             if (ymask == PETSC_OFFLOAD_CPU) {
    1742:               PetscCall(VecCUDAResetArray(*Y));
    1743:             } else if (ymask == PETSC_OFFLOAD_GPU) {
    1744:               X->offloadmask = PETSC_OFFLOAD_GPU;
    1745:             }
    1746:             break;
    1747:           case PETSC_OFFLOAD_GPU:
    1748:             if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
    1749:             break;
    1750:           case PETSC_OFFLOAD_CPU:
    1751:             if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
    1752:             break;
    1753:           case PETSC_OFFLOAD_UNALLOCATED:
    1754:           case PETSC_OFFLOAD_KOKKOS:
    1755:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
    1756:           }
    1757: #endif
    1758:         } else if (iship) {
    1759: #if defined(PETSC_HAVE_HIP)
    1760:           PetscOffloadMask ymask = (*Y)->offloadmask;

    1762:           /* The offloadmask of X dictates where to move memory
    1763:               If X GPU data is valid, then move Y data on GPU if needed
    1764:               Otherwise, move back to the CPU */
    1765:           switch (X->offloadmask) {
    1766:           case PETSC_OFFLOAD_BOTH:
    1767:             if (ymask == PETSC_OFFLOAD_CPU) {
    1768:               PetscCall(VecHIPResetArray(*Y));
    1769:             } else if (ymask == PETSC_OFFLOAD_GPU) {
    1770:               X->offloadmask = PETSC_OFFLOAD_GPU;
    1771:             }
    1772:             break;
    1773:           case PETSC_OFFLOAD_GPU:
    1774:             if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
    1775:             break;
    1776:           case PETSC_OFFLOAD_CPU:
    1777:             if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
    1778:             break;
    1779:           case PETSC_OFFLOAD_UNALLOCATED:
    1780:           case PETSC_OFFLOAD_KOKKOS:
    1781:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
    1782:           }
    1783: #endif
    1784:         } else {
    1785:           /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
    1786:           PetscCall(VecResetArray(*Y));
    1787:         }
    1788:         PetscCall(PetscObjectStateIncrease((PetscObject)X));
    1789:       }
    1790:     }
    1791:   }
    1792:   PetscCall(VecDestroy(Y));
    1793:   PetscFunctionReturn(PETSC_SUCCESS);
    1794: }

    1796: /*@
    1797:   VecCreateLocalVector - Creates a vector object suitable for use with `VecGetLocalVector()` and friends. You must call `VecDestroy()` when the
    1798:   vector is no longer needed.

    1800:   Not Collective.

    1802:   Input Parameter:
    1803: . v - The vector for which the local vector is desired.

    1805:   Output Parameter:
    1806: . w - Upon exit this contains the local vector.

    1808:   Level: beginner

    1810: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
    1811: @*/
    1812: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
    1813: {
    1814:   VecType  roottype;
    1815:   PetscInt n;

    1817:   PetscFunctionBegin;
    1819:   PetscAssertPointer(w, 2);
    1820:   if (v->ops->createlocalvector) {
    1821:     PetscUseTypeMethod(v, createlocalvector, w);
    1822:     PetscFunctionReturn(PETSC_SUCCESS);
    1823:   }
    1824:   PetscCall(VecGetRootType_Private(v, &roottype));
    1825:   PetscCall(VecCreate(PETSC_COMM_SELF, w));
    1826:   PetscCall(VecGetLocalSize(v, &n));
    1827:   PetscCall(VecSetSizes(*w, n, n));
    1828:   PetscCall(VecGetBlockSize(v, &n));
    1829:   PetscCall(VecSetBlockSize(*w, n));
    1830:   PetscCall(VecSetType(*w, roottype));
    1831:   PetscFunctionReturn(PETSC_SUCCESS);
    1832: }

    1834: /*@
    1835:   VecGetLocalVectorRead - Maps the local portion of a vector into a
    1836:   vector.

    1838:   Not Collective.

    1840:   Input Parameter:
    1841: . v - The vector for which the local vector is desired.

    1843:   Output Parameter:
    1844: . w - Upon exit this contains the local vector.

    1846:   Level: beginner

    1848:   Notes:
    1849:   You must call `VecRestoreLocalVectorRead()` when the local
    1850:   vector is no longer needed.

    1852:   This function is similar to `VecGetArrayRead()` which maps the local
    1853:   portion into a raw pointer.  `VecGetLocalVectorRead()` is usually
    1854:   almost as efficient as `VecGetArrayRead()` but in certain circumstances
    1855:   `VecGetLocalVectorRead()` can be much more efficient than
    1856:   `VecGetArrayRead()`.  This is because the construction of a contiguous
    1857:   array representing the vector data required by `VecGetArrayRead()` can
    1858:   be an expensive operation for certain vector types.  For example, for
    1859:   GPU vectors `VecGetArrayRead()` requires that the data between device
    1860:   and host is synchronized.

    1862:   Unlike `VecGetLocalVector()`, this routine is not collective and
    1863:   preserves cached information.

    1865: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
    1866: @*/
    1867: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
    1868: {
    1869:   PetscFunctionBegin;
    1872:   VecCheckSameLocalSize(v, 1, w, 2);
    1873:   if (v->ops->getlocalvectorread) {
    1874:     PetscUseTypeMethod(v, getlocalvectorread, w);
    1875:   } else {
    1876:     PetscScalar *a;

    1878:     PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
    1879:     PetscCall(VecPlaceArray(w, a));
    1880:   }
    1881:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
    1882:   PetscCall(VecLockReadPush(v));
    1883:   PetscCall(VecLockReadPush(w));
    1884:   PetscFunctionReturn(PETSC_SUCCESS);
    1885: }

    1887: /*@
    1888:   VecRestoreLocalVectorRead - Unmaps the local portion of a vector
    1889:   previously mapped into a vector using `VecGetLocalVectorRead()`.

    1891:   Not Collective.

    1893:   Input Parameters:
    1894: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
    1895: - w - The vector into which the local portion of `v` was mapped.

    1897:   Level: beginner

    1899: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
    1900: @*/
    1901: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
    1902: {
    1903:   PetscFunctionBegin;
    1906:   if (v->ops->restorelocalvectorread) {
    1907:     PetscUseTypeMethod(v, restorelocalvectorread, w);
    1908:   } else {
    1909:     const PetscScalar *a;

    1911:     PetscCall(VecGetArrayRead(w, &a));
    1912:     PetscCall(VecRestoreArrayRead(v, &a));
    1913:     PetscCall(VecResetArray(w));
    1914:   }
    1915:   PetscCall(VecLockReadPop(v));
    1916:   PetscCall(VecLockReadPop(w));
    1917:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
    1918:   PetscFunctionReturn(PETSC_SUCCESS);
    1919: }

    1921: /*@
    1922:   VecGetLocalVector - Maps the local portion of a vector into a
    1923:   vector.

    1925:   Collective

    1927:   Input Parameter:
    1928: . v - The vector for which the local vector is desired.

    1930:   Output Parameter:
    1931: . w - Upon exit this contains the local vector.

    1933:   Level: beginner

    1935:   Notes:
    1936:   You must call `VecRestoreLocalVector()` when the local
    1937:   vector is no longer needed.

    1939:   This function is similar to `VecGetArray()` which maps the local
    1940:   portion into a raw pointer.  `VecGetLocalVector()` is usually about as
    1941:   efficient as `VecGetArray()` but in certain circumstances
    1942:   `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
    1943:   This is because the construction of a contiguous array representing
    1944:   the vector data required by `VecGetArray()` can be an expensive
    1945:   operation for certain vector types.  For example, for GPU vectors
    1946:   `VecGetArray()` requires that the data between device and host is
    1947:   synchronized.

    1949: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
    1950: @*/
    1951: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
    1952: {
    1953:   PetscFunctionBegin;
    1956:   VecCheckSameLocalSize(v, 1, w, 2);
    1957:   if (v->ops->getlocalvector) {
    1958:     PetscUseTypeMethod(v, getlocalvector, w);
    1959:   } else {
    1960:     PetscScalar *a;

    1962:     PetscCall(VecGetArray(v, &a));
    1963:     PetscCall(VecPlaceArray(w, a));
    1964:   }
    1965:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
    1966:   PetscFunctionReturn(PETSC_SUCCESS);
    1967: }

    1969: /*@
    1970:   VecRestoreLocalVector - Unmaps the local portion of a vector
    1971:   previously mapped into a vector using `VecGetLocalVector()`.

    1973:   Logically Collective.

    1975:   Input Parameters:
    1976: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
    1977: - w - The vector into which the local portion of `v` was mapped.

    1979:   Level: beginner

    1981: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
    1982: @*/
    1983: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
    1984: {
    1985:   PetscFunctionBegin;
    1988:   if (v->ops->restorelocalvector) {
    1989:     PetscUseTypeMethod(v, restorelocalvector, w);
    1990:   } else {
    1991:     PetscScalar *a;
    1992:     PetscCall(VecGetArray(w, &a));
    1993:     PetscCall(VecRestoreArray(v, &a));
    1994:     PetscCall(VecResetArray(w));
    1995:   }
    1996:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
    1997:   PetscCall(PetscObjectStateIncrease((PetscObject)v));
    1998:   PetscFunctionReturn(PETSC_SUCCESS);
    1999: }

    2001: /*@C
    2002:   VecGetArray - Returns a pointer to a contiguous array that contains this
    2003:   MPI processes's portion of the vector data

    2005:   Logically Collective

    2007:   Input Parameter:
    2008: . x - the vector

    2010:   Output Parameter:
    2011: . a - location to put pointer to the array

    2013:   Level: beginner

    2015:   Notes:
    2016:   For the standard PETSc vectors, `VecGetArray()` returns a pointer to the local data array and
    2017:   does not use any copies. If the underlying vector data is not stored in a contiguous array
    2018:   this routine will copy the data to a contiguous array and return a pointer to that. You MUST
    2019:   call `VecRestoreArray()` when you no longer need access to the array.

    2021:   For vectors that may also have the array data in GPU memory, for example, `VECCUDA`, this call ensures the CPU array has the
    2022:   most recent array values by copying the data from the GPU memory if needed.

    2024:   Fortran Note:
    2025: .vb
    2026:   PetscScalar, pointer :: a(:)
    2027: .ve

    2029: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecPlaceArray()`, `VecGetArray2d()`,
    2030:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`, `VecGetArrayAndMemType()`
    2031: @*/
    2032: PetscErrorCode VecGetArray(Vec x, PetscScalar *a[])
    2033: {
    2034:   PetscFunctionBegin;
    2036:   PetscCall(VecSetErrorIfLocked(x, 1));
    2037:   if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
    2038:     PetscUseTypeMethod(x, getarray, a);
    2039:   } else if (x->petscnative) { /* VECSTANDARD */
    2040:     *a = *((PetscScalar **)x->data);
    2041:   } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
    2042:   PetscFunctionReturn(PETSC_SUCCESS);
    2043: }

    2045: /*@C
    2046:   VecRestoreArray - Restores a vector after `VecGetArray()` has been called and the array is no longer needed

    2048:   Logically Collective

    2050:   Input Parameters:
    2051: + x - the vector
    2052: - a - location of pointer to array obtained from `VecGetArray()`

    2054:   Level: beginner

    2056: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
    2057:           `VecGetArrayPair()`, `VecRestoreArrayPair()`
    2058: @*/
    2059: PetscErrorCode VecRestoreArray(Vec x, PetscScalar *a[])
    2060: {
    2061:   PetscFunctionBegin;
    2063:   if (a) PetscAssertPointer(a, 2);
    2064:   if (x->ops->restorearray) {
    2065:     PetscUseTypeMethod(x, restorearray, a);
    2066:   } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
    2067:   if (a) *a = NULL;
    2068:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
    2069:   PetscFunctionReturn(PETSC_SUCCESS);
    2070: }
    2071: /*@C
    2072:   VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.

    2074:   Not Collective

    2076:   Input Parameter:
    2077: . x - the vector

    2079:   Output Parameter:
    2080: . a - the array

    2082:   Level: beginner

    2084:   Notes:
    2085:   The array must be returned using a matching call to `VecRestoreArrayRead()`.

    2087:   Unlike `VecGetArray()`, preserves cached information like vector norms.

    2089:   Standard PETSc vectors use contiguous storage so that this routine does not perform a copy.  Other vector
    2090:   implementations may require a copy, but such implementations should cache the contiguous representation so that
    2091:   only one copy is performed when this routine is called multiple times in sequence.

    2093:   For vectors that may also have the array data in GPU memory, for example, `VECCUDA`, this call ensures the CPU array has the
    2094:   most recent array values by copying the data from the GPU memory if needed.

    2096: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
    2097:           `VecGetArrayAndMemType()`
    2098: @*/
    2099: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar *a[])
    2100: {
    2101:   PetscFunctionBegin;
    2103:   PetscAssertPointer(a, 2);
    2104:   if (x->ops->getarrayread) {
    2105:     PetscUseTypeMethod(x, getarrayread, a);
    2106:   } else if (x->ops->getarray) {
    2107:     PetscObjectState state;

    2109:     /* VECNEST, VECCUDA, VECKOKKOS etc */
    2110:     // x->ops->getarray may bump the object state, but since we know this is a read-only get
    2111:     // we can just undo that
    2112:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
    2113:     PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
    2114:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
    2115:   } else if (x->petscnative) {
    2116:     /* VECSTANDARD */
    2117:     *a = *((PetscScalar **)x->data);
    2118:   } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
    2119:   PetscFunctionReturn(PETSC_SUCCESS);
    2120: }

    2122: /*@C
    2123:   VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`

    2125:   Not Collective

    2127:   Input Parameters:
    2128: + x - the vector
    2129: - a - the array

    2131:   Level: beginner

    2133: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
    2134: @*/
    2135: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar *a[])
    2136: {
    2137:   PetscFunctionBegin;
    2139:   if (a) PetscAssertPointer(a, 2);
    2140:   if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
    2141:     /* nothing */
    2142:   } else if (x->ops->restorearrayread) { /* VECNEST */
    2143:     PetscUseTypeMethod(x, restorearrayread, a);
    2144:   } else { /* No one? */
    2145:     PetscObjectState state;

    2147:     // x->ops->restorearray may bump the object state, but since we know this is a read-restore
    2148:     // we can just undo that
    2149:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
    2150:     PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
    2151:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
    2152:   }
    2153:   if (a) *a = NULL;
    2154:   PetscFunctionReturn(PETSC_SUCCESS);
    2155: }

    2157: /*@C
    2158:   VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
    2159:   MPI processes's portion of the vector data.

    2161:   Logically Collective

    2163:   Input Parameter:
    2164: . x - the vector

    2166:   Output Parameter:
    2167: . a - location to put pointer to the array

    2169:   Level: intermediate

    2171:   Note:
    2172:   The values in this array are NOT valid, the caller of this routine is responsible for putting
    2173:   values into the array; any values it does not set will be invalid.

    2175:   The array must be returned using a matching call to `VecRestoreArrayWrite()`.

    2177:   For vectors associated with GPUs, the host and device vectors are not synchronized before
    2178:   giving access. If you need correct values in the array use `VecGetArray()`

    2180: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecPlaceArray()`, `VecGetArray2d()`,
    2181:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`, `VecGetArrayAndMemType()`
    2182: @*/
    2183: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar *a[])
    2184: {
    2185:   PetscFunctionBegin;
    2187:   PetscAssertPointer(a, 2);
    2188:   PetscCall(VecSetErrorIfLocked(x, 1));
    2189:   if (x->ops->getarraywrite) {
    2190:     PetscUseTypeMethod(x, getarraywrite, a);
    2191:   } else {
    2192:     PetscCall(VecGetArray(x, a));
    2193:   }
    2194:   PetscFunctionReturn(PETSC_SUCCESS);
    2195: }

    2197: /*@C
    2198:   VecRestoreArrayWrite - Restores a vector after `VecGetArrayWrite()` has been called.

    2200:   Logically Collective

    2202:   Input Parameters:
    2203: + x - the vector
    2204: - a - location of pointer to array obtained from `VecGetArray()`

    2206:   Level: beginner

    2208: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
    2209:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
    2210: @*/
    2211: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar *a[])
    2212: {
    2213:   PetscFunctionBegin;
    2215:   if (a) PetscAssertPointer(a, 2);
    2216:   if (x->ops->restorearraywrite) {
    2217:     PetscUseTypeMethod(x, restorearraywrite, a);
    2218:   } else if (x->ops->restorearray) {
    2219:     PetscUseTypeMethod(x, restorearray, a);
    2220:   }
    2221:   if (a) *a = NULL;
    2222:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
    2223:   PetscFunctionReturn(PETSC_SUCCESS);
    2224: }

    2226: /*@C
    2227:   VecGetArrays - Returns a pointer to the arrays in a set of vectors
    2228:   that were created by a call to `VecDuplicateVecs()`.

    2230:   Logically Collective; No Fortran Support

    2232:   Input Parameters:
    2233: + x - the vectors
    2234: - n - the number of vectors

    2236:   Output Parameter:
    2237: . a - location to put pointer to the array

    2239:   Level: intermediate

    2241:   Note:
    2242:   You MUST call `VecRestoreArrays()` when you no longer need access to the arrays.

    2244: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
    2245: @*/
    2246: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
    2247: {
    2248:   PetscInt      i;
    2249:   PetscScalar **q;

    2251:   PetscFunctionBegin;
    2252:   PetscAssertPointer(x, 1);
    2254:   PetscAssertPointer(a, 3);
    2255:   PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
    2256:   PetscCall(PetscMalloc1(n, &q));
    2257:   for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
    2258:   *a = q;
    2259:   PetscFunctionReturn(PETSC_SUCCESS);
    2260: }

    2262: /*@C
    2263:   VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
    2264:   has been called.

    2266:   Logically Collective; No Fortran Support

    2268:   Input Parameters:
    2269: + x - the vector
    2270: . n - the number of vectors
    2271: - a - location of pointer to arrays obtained from `VecGetArrays()`

    2273:   Notes:
    2274:   For regular PETSc vectors this routine does not involve any copies. For
    2275:   any special vectors that do not store local vector data in a contiguous
    2276:   array, this routine will copy the data back into the underlying
    2277:   vector data structure from the arrays obtained with `VecGetArrays()`.

    2279:   Level: intermediate

    2281: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
    2282: @*/
    2283: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
    2284: {
    2285:   PetscInt      i;
    2286:   PetscScalar **q = *a;

    2288:   PetscFunctionBegin;
    2289:   PetscAssertPointer(x, 1);
    2291:   PetscAssertPointer(a, 3);

    2293:   for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
    2294:   PetscCall(PetscFree(q));
    2295:   PetscFunctionReturn(PETSC_SUCCESS);
    2296: }

    2298: /*@C
    2299:   VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g.,
    2300:   `VECCUDA`), the returned pointer will be a device pointer to the device memory that contains
    2301:   this MPI processes's portion of the vector data.

    2303:   Logically Collective; No Fortran Support

    2305:   Input Parameter:
    2306: . x - the vector

    2308:   Output Parameters:
    2309: + a     - location to put pointer to the array
    2310: - mtype - memory type of the array

    2312:   Level: beginner

    2314:   Note:
    2315:   Device data is guaranteed to have the latest value. Otherwise, when this is a host vector
    2316:   (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host
    2317:   pointer.

    2319:   For `VECKOKKOS`, if Kokkos is configured without device (e.g., use serial or openmp), per
    2320:   this function, the vector works like `VECSEQ`/`VECMPI`; otherwise, it works like `VECCUDA` or
    2321:   `VECHIP` etc.

    2323:   Use `VecRestoreArrayAndMemType()` when the array access is no longer needed.

    2325: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`,
    2326:           `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
    2327: @*/
    2328: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar *a[], PetscMemType *mtype)
    2329: {
    2330:   PetscFunctionBegin;
    2333:   if (a) PetscAssertPointer(a, 2);
    2334:   if (mtype) PetscAssertPointer(mtype, 3);
    2335:   PetscCall(VecSetErrorIfLocked(x, 1));
    2336:   if (x->ops->getarrayandmemtype) {
    2337:     /* VECCUDA, VECKOKKOS etc */
    2338:     PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
    2339:   } else {
    2340:     /* VECSTANDARD, VECNEST, VECVIENNACL */
    2341:     PetscCall(VecGetArray(x, a));
    2342:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
    2343:   }
    2344:   PetscFunctionReturn(PETSC_SUCCESS);
    2345: }

    2347: /*@C
    2348:   VecRestoreArrayAndMemType - Restores a vector after `VecGetArrayAndMemType()` has been called.

    2350:   Logically Collective; No Fortran Support

    2352:   Input Parameters:
    2353: + x - the vector
    2354: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`

    2356:   Level: beginner

    2358: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`,
    2359:           `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
    2360: @*/
    2361: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar *a[])
    2362: {
    2363:   PetscFunctionBegin;
    2366:   if (a) PetscAssertPointer(a, 2);
    2367:   if (x->ops->restorearrayandmemtype) {
    2368:     /* VECCUDA, VECKOKKOS etc */
    2369:     PetscUseTypeMethod(x, restorearrayandmemtype, a);
    2370:   } else {
    2371:     /* VECNEST, VECVIENNACL */
    2372:     PetscCall(VecRestoreArray(x, a));
    2373:   } /* VECSTANDARD does nothing */
    2374:   if (a) *a = NULL;
    2375:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
    2376:   PetscFunctionReturn(PETSC_SUCCESS);
    2377: }

    2379: /*@C
    2380:   VecGetArrayReadAndMemType - Like `VecGetArrayRead()`, but if the input vector is a device vector, it will return a read-only device pointer.
    2381:   The returned pointer is guaranteed to point to up-to-date data. For host vectors, it functions as `VecGetArrayRead()`.

    2383:   Not Collective; No Fortran Support

    2385:   Input Parameter:
    2386: . x - the vector

    2388:   Output Parameters:
    2389: + a     - the array
    2390: - mtype - memory type of the array

    2392:   Level: beginner

    2394:   Notes:
    2395:   The array must be returned using a matching call to `VecRestoreArrayReadAndMemType()`.

    2397: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
    2398: @*/
    2399: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar *a[], PetscMemType *mtype)
    2400: {
    2401:   PetscFunctionBegin;
    2404:   PetscAssertPointer(a, 2);
    2405:   if (mtype) PetscAssertPointer(mtype, 3);
    2406:   if (x->ops->getarrayreadandmemtype) {
    2407:     /* VECCUDA/VECHIP though they are also petscnative */
    2408:     PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
    2409:   } else if (x->ops->getarrayandmemtype) {
    2410:     /* VECKOKKOS */
    2411:     PetscObjectState state;

    2413:     // see VecGetArrayRead() for why
    2414:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
    2415:     PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
    2416:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
    2417:   } else {
    2418:     PetscCall(VecGetArrayRead(x, a));
    2419:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
    2420:   }
    2421:   PetscFunctionReturn(PETSC_SUCCESS);
    2422: }

    2424: /*@C
    2425:   VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`

    2427:   Not Collective; No Fortran Support

    2429:   Input Parameters:
    2430: + x - the vector
    2431: - a - the array

    2433:   Level: beginner

    2435: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
    2436: @*/
    2437: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar *a[])
    2438: {
    2439:   PetscFunctionBegin;
    2442:   if (a) PetscAssertPointer(a, 2);
    2443:   if (x->ops->restorearrayreadandmemtype) {
    2444:     /* VECCUDA/VECHIP */
    2445:     PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
    2446:   } else if (!x->petscnative) {
    2447:     /* VECNEST */
    2448:     PetscCall(VecRestoreArrayRead(x, a));
    2449:   }
    2450:   if (a) *a = NULL;
    2451:   PetscFunctionReturn(PETSC_SUCCESS);
    2452: }

    2454: /*@C
    2455:   VecGetArrayWriteAndMemType - Like `VecGetArrayWrite()`, but if this is a device vector it will always return
    2456:   a device pointer to the device memory that contains this processor's portion of the vector data.

    2458:   Logically Collective; No Fortran Support

    2460:   Input Parameter:
    2461: . x - the vector

    2463:   Output Parameters:
    2464: + a     - the array
    2465: - mtype - memory type of the array

    2467:   Level: beginner

    2469:   Note:
    2470:   The array must be returned using a matching call to `VecRestoreArrayWriteAndMemType()`, where it will label the device memory as most recent.

    2472: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
    2473: @*/
    2474: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar *a[], PetscMemType *mtype)
    2475: {
    2476:   PetscFunctionBegin;
    2479:   PetscCall(VecSetErrorIfLocked(x, 1));
    2480:   PetscAssertPointer(a, 2);
    2481:   if (mtype) PetscAssertPointer(mtype, 3);
    2482:   if (x->ops->getarraywriteandmemtype) {
    2483:     /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
    2484:     PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
    2485:   } else if (x->ops->getarrayandmemtype) {
    2486:     PetscCall(VecGetArrayAndMemType(x, a, mtype));
    2487:   } else {
    2488:     /* VECNEST, VECVIENNACL */
    2489:     PetscCall(VecGetArrayWrite(x, a));
    2490:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
    2491:   }
    2492:   PetscFunctionReturn(PETSC_SUCCESS);
    2493: }

    2495: /*@C
    2496:   VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`

    2498:   Logically Collective; No Fortran Support

    2500:   Input Parameters:
    2501: + x - the vector
    2502: - a - the array

    2504:   Level: beginner

    2506: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
    2507: @*/
    2508: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar *a[])
    2509: {
    2510:   PetscFunctionBegin;
    2513:   PetscCall(VecSetErrorIfLocked(x, 1));
    2514:   if (a) PetscAssertPointer(a, 2);
    2515:   if (x->ops->restorearraywriteandmemtype) {
    2516:     /* VECCUDA/VECHIP */
    2517:     PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
    2518:     PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
    2519:   } else if (x->ops->restorearrayandmemtype) {
    2520:     PetscCall(VecRestoreArrayAndMemType(x, a));
    2521:   } else {
    2522:     PetscCall(VecRestoreArray(x, a));
    2523:   }
    2524:   if (a) *a = NULL;
    2525:   PetscFunctionReturn(PETSC_SUCCESS);
    2526: }

    2528: /*@
    2529:   VecPlaceArray - Allows one to replace the array in a vector with an
    2530:   array provided by the user. This is useful to avoid copying an array
    2531:   into a vector.

    2533:   Logically Collective

    2535:   Input Parameters:
    2536: + vec   - the vector
    2537: - array - the array

    2539:   Level: developer

    2541:   Notes:
    2542:   Adding `const` to `array` was an oversight, as subsequent operations on `vec` would
    2543:   likely modify the data in `array`. However, we have kept it to avoid breaking APIs.

    2545:   Use `VecReplaceArray()` instead to permanently replace the array

    2547:   You can return to the original array with a call to `VecResetArray()`. `vec` does not take
    2548:   ownership of `array` in any way.

    2550:   The user must free `array` themselves but be careful not to
    2551:   do so before the vector has either been destroyed, had its original array restored with
    2552:   `VecResetArray()` or permanently replaced with `VecReplaceArray()`.

    2554: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
    2555: @*/
    2556: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
    2557: {
    2558:   PetscFunctionBegin;
    2561:   if (array) PetscAssertPointer(array, 2);
    2562:   PetscUseTypeMethod(vec, placearray, array);
    2563:   PetscCall(PetscObjectStateIncrease((PetscObject)vec));
    2564:   PetscFunctionReturn(PETSC_SUCCESS);
    2565: }

    2567: /*@C
    2568:   VecReplaceArray - Allows one to replace the array in a vector with an
    2569:   array provided by the user. This is useful to avoid copying an array
    2570:   into a vector.

    2572:   Logically Collective; No Fortran Support

    2574:   Input Parameters:
    2575: + vec   - the vector
    2576: - array - the array

    2578:   Level: developer

    2580:   Notes:
    2581:   Adding `const` to `array` was an oversight, as subsequent operations on `vec` would
    2582:   likely modify the data in `array`. However, we have kept it to avoid breaking APIs.

    2584:   This permanently replaces the array and frees the memory associated
    2585:   with the old array. Use `VecPlaceArray()` to temporarily replace the array.

    2587:   The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
    2588:   freed by the user. It will be freed when the vector is destroyed.

    2590: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
    2591: @*/
    2592: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
    2593: {
    2594:   PetscFunctionBegin;
    2597:   PetscUseTypeMethod(vec, replacearray, array);
    2598:   PetscCall(PetscObjectStateIncrease((PetscObject)vec));
    2599:   PetscFunctionReturn(PETSC_SUCCESS);
    2600: }

    2602: /*@C
    2603:   VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
    2604:   processor's portion of the vector data.  You MUST call `VecRestoreArray2d()`
    2605:   when you no longer need access to the array.

    2607:   Logically Collective

    2609:   Input Parameters:
    2610: + x      - the vector
    2611: . m      - first dimension of two dimensional array
    2612: . n      - second dimension of two dimensional array
    2613: . mstart - first index you will use in first coordinate direction (often 0)
    2614: - nstart - first index in the second coordinate direction (often 0)

    2616:   Output Parameter:
    2617: . a - location to put pointer to the array

    2619:   Level: developer

    2621:   Notes:
    2622:   For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
    2623:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
    2624:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
    2625:   the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.

    2627:   For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    2629: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
    2630:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
    2631:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    2632: @*/
    2633: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
    2634: {
    2635:   PetscInt     i, N;
    2636:   PetscScalar *aa;

    2638:   PetscFunctionBegin;
    2640:   PetscAssertPointer(a, 6);
    2642:   PetscCall(VecGetLocalSize(x, &N));
    2643:   PetscCheck(m * n == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 2d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n);
    2644:   PetscCall(VecGetArray(x, &aa));

    2646:   PetscCall(PetscMalloc1(m, a));
    2647:   for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
    2648:   *a -= mstart;
    2649:   PetscFunctionReturn(PETSC_SUCCESS);
    2650: }

    2652: /*@C
    2653:   VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
    2654:   processor's portion of the vector data.  You MUST call `VecRestoreArray2dWrite()`
    2655:   when you no longer need access to the array.

    2657:   Logically Collective

    2659:   Input Parameters:
    2660: + x      - the vector
    2661: . m      - first dimension of two dimensional array
    2662: . n      - second dimension of two dimensional array
    2663: . mstart - first index you will use in first coordinate direction (often 0)
    2664: - nstart - first index in the second coordinate direction (often 0)

    2666:   Output Parameter:
    2667: . a - location to put pointer to the array

    2669:   Level: developer

    2671:   Notes:
    2672:   For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
    2673:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
    2674:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
    2675:   the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.

    2677:   For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    2679: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
    2680:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
    2681:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    2682: @*/
    2683: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
    2684: {
    2685:   PetscInt     i, N;
    2686:   PetscScalar *aa;

    2688:   PetscFunctionBegin;
    2690:   PetscAssertPointer(a, 6);
    2692:   PetscCall(VecGetLocalSize(x, &N));
    2693:   PetscCheck(m * n == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 2d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n);
    2694:   PetscCall(VecGetArrayWrite(x, &aa));

    2696:   PetscCall(PetscMalloc1(m, a));
    2697:   for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
    2698:   *a -= mstart;
    2699:   PetscFunctionReturn(PETSC_SUCCESS);
    2700: }

    2702: /*@C
    2703:   VecRestoreArray2d - Restores a vector after `VecGetArray2d()` has been called.

    2705:   Logically Collective

    2707:   Input Parameters:
    2708: + x      - the vector
    2709: . m      - first dimension of two dimensional array
    2710: . n      - second dimension of the two dimensional array
    2711: . mstart - first index you will use in first coordinate direction (often 0)
    2712: . nstart - first index in the second coordinate direction (often 0)
    2713: - a      - location of pointer to array obtained from `VecGetArray2d()`

    2715:   Level: developer

    2717:   Notes:
    2718:   For regular PETSc vectors this routine does not involve any copies. For
    2719:   any special vectors that do not store local vector data in a contiguous
    2720:   array, this routine will copy the data back into the underlying
    2721:   vector data structure from the array obtained with `VecGetArray()`.

    2723:   This routine actually zeros out the `a` pointer.

    2725: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
    2726:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
    2727:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    2728: @*/
    2729: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
    2730: {
    2731:   void *dummy;

    2733:   PetscFunctionBegin;
    2735:   PetscAssertPointer(a, 6);
    2737:   dummy = (void *)(*a + mstart);
    2738:   PetscCall(PetscFree(dummy));
    2739:   PetscCall(VecRestoreArray(x, NULL));
    2740:   *a = NULL;
    2741:   PetscFunctionReturn(PETSC_SUCCESS);
    2742: }

    2744: /*@C
    2745:   VecRestoreArray2dWrite - Restores a vector after `VecGetArray2dWrite()` has been called.

    2747:   Logically Collective

    2749:   Input Parameters:
    2750: + x      - the vector
    2751: . m      - first dimension of two dimensional array
    2752: . n      - second dimension of the two dimensional array
    2753: . mstart - first index you will use in first coordinate direction (often 0)
    2754: . nstart - first index in the second coordinate direction (often 0)
    2755: - a      - location of pointer to array obtained from `VecGetArray2d()`

    2757:   Level: developer

    2759:   Notes:
    2760:   For regular PETSc vectors this routine does not involve any copies. For
    2761:   any special vectors that do not store local vector data in a contiguous
    2762:   array, this routine will copy the data back into the underlying
    2763:   vector data structure from the array obtained with `VecGetArray()`.

    2765:   This routine actually zeros out the `a` pointer.

    2767: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
    2768:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
    2769:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    2770: @*/
    2771: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
    2772: {
    2773:   void *dummy;

    2775:   PetscFunctionBegin;
    2777:   PetscAssertPointer(a, 6);
    2779:   dummy = (void *)(*a + mstart);
    2780:   PetscCall(PetscFree(dummy));
    2781:   PetscCall(VecRestoreArrayWrite(x, NULL));
    2782:   PetscFunctionReturn(PETSC_SUCCESS);
    2783: }

    2785: /*@C
    2786:   VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
    2787:   processor's portion of the vector data.  You MUST call `VecRestoreArray1d()`
    2788:   when you no longer need access to the array.

    2790:   Logically Collective

    2792:   Input Parameters:
    2793: + x      - the vector
    2794: . m      - first dimension of two dimensional array
    2795: - mstart - first index you will use in first coordinate direction (often 0)

    2797:   Output Parameter:
    2798: . a - location to put pointer to the array

    2800:   Level: developer

    2802:   Notes:
    2803:   For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
    2804:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
    2805:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

    2807:   For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    2809: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
    2810:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
    2811:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    2812: @*/
    2813: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
    2814: {
    2815:   PetscInt N;

    2817:   PetscFunctionBegin;
    2819:   PetscAssertPointer(a, 4);
    2821:   PetscCall(VecGetLocalSize(x, &N));
    2822:   PetscCheck(m == N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local array size %" PetscInt_FMT " does not match 1d array dimensions %" PetscInt_FMT, N, m);
    2823:   PetscCall(VecGetArray(x, a));
    2824:   *a -= mstart;
    2825:   PetscFunctionReturn(PETSC_SUCCESS);
    2826: }

    2828: /*@C
    2829:   VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
    2830:   processor's portion of the vector data.  You MUST call `VecRestoreArray1dWrite()`
    2831:   when you no longer need access to the array.

    2833:   Logically Collective

    2835:   Input Parameters:
    2836: + x      - the vector
    2837: . m      - first dimension of two dimensional array
    2838: - mstart - first index you will use in first coordinate direction (often 0)

    2840:   Output Parameter:
    2841: . a - location to put pointer to the array

    2843:   Level: developer

    2845:   Notes:
    2846:   For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
    2847:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
    2848:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

    2850:   For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    2852: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
    2853:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
    2854:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    2855: @*/
    2856: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
    2857: {
    2858:   PetscInt N;

    2860:   PetscFunctionBegin;
    2862:   PetscAssertPointer(a, 4);
    2864:   PetscCall(VecGetLocalSize(x, &N));
    2865:   PetscCheck(m == N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local array size %" PetscInt_FMT " does not match 1d array dimensions %" PetscInt_FMT, N, m);
    2866:   PetscCall(VecGetArrayWrite(x, a));
    2867:   *a -= mstart;
    2868:   PetscFunctionReturn(PETSC_SUCCESS);
    2869: }

    2871: /*@C
    2872:   VecRestoreArray1d - Restores a vector after `VecGetArray1d()` has been called.

    2874:   Logically Collective

    2876:   Input Parameters:
    2877: + x      - the vector
    2878: . m      - first dimension of two dimensional array
    2879: . mstart - first index you will use in first coordinate direction (often 0)
    2880: - a      - location of pointer to array obtained from `VecGetArray1d()`

    2882:   Level: developer

    2884:   Notes:
    2885:   For regular PETSc vectors this routine does not involve any copies. For
    2886:   any special vectors that do not store local vector data in a contiguous
    2887:   array, this routine will copy the data back into the underlying
    2888:   vector data structure from the array obtained with `VecGetArray1d()`.

    2890:   This routine actually zeros out the `a` pointer.

    2892: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
    2893:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
    2894:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    2895: @*/
    2896: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
    2897: {
    2898:   PetscFunctionBegin;
    2901:   PetscCall(VecRestoreArray(x, NULL));
    2902:   *a = NULL;
    2903:   PetscFunctionReturn(PETSC_SUCCESS);
    2904: }

    2906: /*@C
    2907:   VecRestoreArray1dWrite - Restores a vector after `VecGetArray1dWrite()` has been called.

    2909:   Logically Collective

    2911:   Input Parameters:
    2912: + x      - the vector
    2913: . m      - first dimension of two dimensional array
    2914: . mstart - first index you will use in first coordinate direction (often 0)
    2915: - a      - location of pointer to array obtained from `VecGetArray1d()`

    2917:   Level: developer

    2919:   Notes:
    2920:   For regular PETSc vectors this routine does not involve any copies. For
    2921:   any special vectors that do not store local vector data in a contiguous
    2922:   array, this routine will copy the data back into the underlying
    2923:   vector data structure from the array obtained with `VecGetArray1d()`.

    2925:   This routine actually zeros out the `a` pointer.

    2927: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
    2928:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
    2929:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    2930: @*/
    2931: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
    2932: {
    2933:   PetscFunctionBegin;
    2936:   PetscCall(VecRestoreArrayWrite(x, NULL));
    2937:   *a = NULL;
    2938:   PetscFunctionReturn(PETSC_SUCCESS);
    2939: }

    2941: /*@C
    2942:   VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
    2943:   processor's portion of the vector data.  You MUST call `VecRestoreArray3d()`
    2944:   when you no longer need access to the array.

    2946:   Logically Collective

    2948:   Input Parameters:
    2949: + x      - the vector
    2950: . m      - first dimension of three dimensional array
    2951: . n      - second dimension of three dimensional array
    2952: . p      - third dimension of three dimensional array
    2953: . mstart - first index you will use in first coordinate direction (often 0)
    2954: . nstart - first index in the second coordinate direction (often 0)
    2955: - pstart - first index in the third coordinate direction (often 0)

    2957:   Output Parameter:
    2958: . a - location to put pointer to the array

    2960:   Level: developer

    2962:   Notes:
    2963:   For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
    2964:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
    2965:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
    2966:   the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.

    2968:   For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    2970: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
    2971:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
    2972:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    2973: @*/
    2974: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
    2975: {
    2976:   PetscInt     i, N, j;
    2977:   PetscScalar *aa, **b;

    2979:   PetscFunctionBegin;
    2981:   PetscAssertPointer(a, 8);
    2983:   PetscCall(VecGetLocalSize(x, &N));
    2984:   PetscCheck(m * n * p == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 3d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n, p);
    2985:   PetscCall(VecGetArray(x, &aa));

    2987:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
    2988:   b = (PetscScalar **)((*a) + m);
    2989:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
    2990:   for (i = 0; i < m; i++)
    2991:     for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset(aa, i * n * p + j * p - pstart);
    2992:   *a -= mstart;
    2993:   PetscFunctionReturn(PETSC_SUCCESS);
    2994: }

    2996: /*@C
    2997:   VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
    2998:   processor's portion of the vector data.  You MUST call `VecRestoreArray3dWrite()`
    2999:   when you no longer need access to the array.

    3001:   Logically Collective

    3003:   Input Parameters:
    3004: + x      - the vector
    3005: . m      - first dimension of three dimensional array
    3006: . n      - second dimension of three dimensional array
    3007: . p      - third dimension of three dimensional array
    3008: . mstart - first index you will use in first coordinate direction (often 0)
    3009: . nstart - first index in the second coordinate direction (often 0)
    3010: - pstart - first index in the third coordinate direction (often 0)

    3012:   Output Parameter:
    3013: . a - location to put pointer to the array

    3015:   Level: developer

    3017:   Notes:
    3018:   For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
    3019:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
    3020:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
    3021:   the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.

    3023:   For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    3025: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
    3026:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
    3027:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    3028: @*/
    3029: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
    3030: {
    3031:   PetscInt     i, N, j;
    3032:   PetscScalar *aa, **b;

    3034:   PetscFunctionBegin;
    3036:   PetscAssertPointer(a, 8);
    3038:   PetscCall(VecGetLocalSize(x, &N));
    3039:   PetscCheck(m * n * p == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 3d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n, p);
    3040:   PetscCall(VecGetArrayWrite(x, &aa));

    3042:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
    3043:   b = (PetscScalar **)((*a) + m);
    3044:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
    3045:   for (i = 0; i < m; i++)
    3046:     for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;

    3048:   *a -= mstart;
    3049:   PetscFunctionReturn(PETSC_SUCCESS);
    3050: }

    3052: /*@C
    3053:   VecRestoreArray3d - Restores a vector after `VecGetArray3d()` has been called.

    3055:   Logically Collective

    3057:   Input Parameters:
    3058: + x      - the vector
    3059: . m      - first dimension of three dimensional array
    3060: . n      - second dimension of the three dimensional array
    3061: . p      - third dimension of the three dimensional array
    3062: . mstart - first index you will use in first coordinate direction (often 0)
    3063: . nstart - first index in the second coordinate direction (often 0)
    3064: . pstart - first index in the third coordinate direction (often 0)
    3065: - a      - location of pointer to array obtained from VecGetArray3d()

    3067:   Level: developer

    3069:   Notes:
    3070:   For regular PETSc vectors this routine does not involve any copies. For
    3071:   any special vectors that do not store local vector data in a contiguous
    3072:   array, this routine will copy the data back into the underlying
    3073:   vector data structure from the array obtained with `VecGetArray()`.

    3075:   This routine actually zeros out the `a` pointer.

    3077: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
    3078:           `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
    3079:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    3080: @*/
    3081: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
    3082: {
    3083:   void *dummy;

    3085:   PetscFunctionBegin;
    3087:   PetscAssertPointer(a, 8);
    3089:   dummy = (void *)(*a + mstart);
    3090:   PetscCall(PetscFree(dummy));
    3091:   PetscCall(VecRestoreArray(x, NULL));
    3092:   *a = NULL;
    3093:   PetscFunctionReturn(PETSC_SUCCESS);
    3094: }

    3096: /*@C
    3097:   VecRestoreArray3dWrite - Restores a vector after `VecGetArray3dWrite()` has been called.

    3099:   Logically Collective

    3101:   Input Parameters:
    3102: + x      - the vector
    3103: . m      - first dimension of three dimensional array
    3104: . n      - second dimension of the three dimensional array
    3105: . p      - third dimension of the three dimensional array
    3106: . mstart - first index you will use in first coordinate direction (often 0)
    3107: . nstart - first index in the second coordinate direction (often 0)
    3108: . pstart - first index in the third coordinate direction (often 0)
    3109: - a      - location of pointer to array obtained from VecGetArray3d()

    3111:   Level: developer

    3113:   Notes:
    3114:   For regular PETSc vectors this routine does not involve any copies. For
    3115:   any special vectors that do not store local vector data in a contiguous
    3116:   array, this routine will copy the data back into the underlying
    3117:   vector data structure from the array obtained with `VecGetArray()`.

    3119:   This routine actually zeros out the `a` pointer.

    3121: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
    3122:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
    3123:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    3124: @*/
    3125: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
    3126: {
    3127:   void *dummy;

    3129:   PetscFunctionBegin;
    3131:   PetscAssertPointer(a, 8);
    3133:   dummy = (void *)(*a + mstart);
    3134:   PetscCall(PetscFree(dummy));
    3135:   PetscCall(VecRestoreArrayWrite(x, NULL));
    3136:   *a = NULL;
    3137:   PetscFunctionReturn(PETSC_SUCCESS);
    3138: }

    3140: /*@C
    3141:   VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
    3142:   processor's portion of the vector data.  You MUST call `VecRestoreArray4d()`
    3143:   when you no longer need access to the array.

    3145:   Logically Collective

    3147:   Input Parameters:
    3148: + x      - the vector
    3149: . m      - first dimension of four dimensional array
    3150: . n      - second dimension of four dimensional array
    3151: . p      - third dimension of four dimensional array
    3152: . q      - fourth dimension of four dimensional array
    3153: . mstart - first index you will use in first coordinate direction (often 0)
    3154: . nstart - first index in the second coordinate direction (often 0)
    3155: . pstart - first index in the third coordinate direction (often 0)
    3156: - qstart - first index in the fourth coordinate direction (often 0)

    3158:   Output Parameter:
    3159: . a - location to put pointer to the array

    3161:   Level: developer

    3163:   Notes:
    3164:   For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
    3165:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
    3166:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
    3167:   the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.

    3169:   For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    3171: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
    3172:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
    3173:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
    3174: @*/
    3175: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
    3176: {
    3177:   PetscInt     i, N, j, k;
    3178:   PetscScalar *aa, ***b, **c;

    3180:   PetscFunctionBegin;
    3182:   PetscAssertPointer(a, 10);
    3184:   PetscCall(VecGetLocalSize(x, &N));
    3185:   PetscCheck(m * n * p * q == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 4d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n, p, q);
    3186:   PetscCall(VecGetArray(x, &aa));

    3188:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
    3189:   b = (PetscScalar ***)((*a) + m);
    3190:   c = (PetscScalar **)(b + m * n);
    3191:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
    3192:   for (i = 0; i < m; i++)
    3193:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
    3194:   for (i = 0; i < m; i++)
    3195:     for (j = 0; j < n; j++)
    3196:       for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
    3197:   *a -= mstart;
    3198:   PetscFunctionReturn(PETSC_SUCCESS);
    3199: }

    3201: /*@C
    3202:   VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
    3203:   processor's portion of the vector data.  You MUST call `VecRestoreArray4dWrite()`
    3204:   when you no longer need access to the array.

    3206:   Logically Collective

    3208:   Input Parameters:
    3209: + x      - the vector
    3210: . m      - first dimension of four dimensional array
    3211: . n      - second dimension of four dimensional array
    3212: . p      - third dimension of four dimensional array
    3213: . q      - fourth dimension of four dimensional array
    3214: . mstart - first index you will use in first coordinate direction (often 0)
    3215: . nstart - first index in the second coordinate direction (often 0)
    3216: . pstart - first index in the third coordinate direction (often 0)
    3217: - qstart - first index in the fourth coordinate direction (often 0)

    3219:   Output Parameter:
    3220: . a - location to put pointer to the array

    3222:   Level: developer

    3224:   Notes:
    3225:   For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
    3226:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
    3227:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
    3228:   the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.

    3230:   For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    3232: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
    3233:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
    3234:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    3235: @*/
    3236: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
    3237: {
    3238:   PetscInt     i, N, j, k;
    3239:   PetscScalar *aa, ***b, **c;

    3241:   PetscFunctionBegin;
    3243:   PetscAssertPointer(a, 10);
    3245:   PetscCall(VecGetLocalSize(x, &N));
    3246:   PetscCheck(m * n * p * q == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 4d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n, p, q);
    3247:   PetscCall(VecGetArrayWrite(x, &aa));

    3249:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
    3250:   b = (PetscScalar ***)((*a) + m);
    3251:   c = (PetscScalar **)(b + m * n);
    3252:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
    3253:   for (i = 0; i < m; i++)
    3254:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
    3255:   for (i = 0; i < m; i++)
    3256:     for (j = 0; j < n; j++)
    3257:       for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
    3258:   *a -= mstart;
    3259:   PetscFunctionReturn(PETSC_SUCCESS);
    3260: }

    3262: /*@C
    3263:   VecRestoreArray4d - Restores a vector after `VecGetArray4d()` has been called.

    3265:   Logically Collective

    3267:   Input Parameters:
    3268: + x      - the vector
    3269: . m      - first dimension of four dimensional array
    3270: . n      - second dimension of the four dimensional array
    3271: . p      - third dimension of the four dimensional array
    3272: . q      - fourth dimension of the four dimensional array
    3273: . mstart - first index you will use in first coordinate direction (often 0)
    3274: . nstart - first index in the second coordinate direction (often 0)
    3275: . pstart - first index in the third coordinate direction (often 0)
    3276: . qstart - first index in the fourth coordinate direction (often 0)
    3277: - a      - location of pointer to array obtained from VecGetArray4d()

    3279:   Level: developer

    3281:   Notes:
    3282:   For regular PETSc vectors this routine does not involve any copies. For
    3283:   any special vectors that do not store local vector data in a contiguous
    3284:   array, this routine will copy the data back into the underlying
    3285:   vector data structure from the array obtained with `VecGetArray()`.

    3287:   This routine actually zeros out the `a` pointer.

    3289: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
    3290:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
    3291:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`
    3292: @*/
    3293: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
    3294: {
    3295:   void *dummy;

    3297:   PetscFunctionBegin;
    3299:   PetscAssertPointer(a, 10);
    3301:   dummy = (void *)(*a + mstart);
    3302:   PetscCall(PetscFree(dummy));
    3303:   PetscCall(VecRestoreArray(x, NULL));
    3304:   *a = NULL;
    3305:   PetscFunctionReturn(PETSC_SUCCESS);
    3306: }

    3308: /*@C
    3309:   VecRestoreArray4dWrite - Restores a vector after `VecGetArray4dWrite()` has been called.

    3311:   Logically Collective

    3313:   Input Parameters:
    3314: + x      - the vector
    3315: . m      - first dimension of four dimensional array
    3316: . n      - second dimension of the four dimensional array
    3317: . p      - third dimension of the four dimensional array
    3318: . q      - fourth dimension of the four dimensional array
    3319: . mstart - first index you will use in first coordinate direction (often 0)
    3320: . nstart - first index in the second coordinate direction (often 0)
    3321: . pstart - first index in the third coordinate direction (often 0)
    3322: . qstart - first index in the fourth coordinate direction (often 0)
    3323: - a      - location of pointer to array obtained from `VecGetArray4d()`

    3325:   Level: developer

    3327:   Notes:
    3328:   For regular PETSc vectors this routine does not involve any copies. For
    3329:   any special vectors that do not store local vector data in a contiguous
    3330:   array, this routine will copy the data back into the underlying
    3331:   vector data structure from the array obtained with `VecGetArray()`.

    3333:   This routine actually zeros out the `a` pointer.

    3335: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
    3336:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
    3337:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    3338: @*/
    3339: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
    3340: {
    3341:   void *dummy;

    3343:   PetscFunctionBegin;
    3345:   PetscAssertPointer(a, 10);
    3347:   dummy = (void *)(*a + mstart);
    3348:   PetscCall(PetscFree(dummy));
    3349:   PetscCall(VecRestoreArrayWrite(x, NULL));
    3350:   *a = NULL;
    3351:   PetscFunctionReturn(PETSC_SUCCESS);
    3352: }

    3354: /*@C
    3355:   VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
    3356:   processor's portion of the vector data.  You MUST call `VecRestoreArray2dRead()`
    3357:   when you no longer need access to the array.

    3359:   Logically Collective

    3361:   Input Parameters:
    3362: + x      - the vector
    3363: . m      - first dimension of two dimensional array
    3364: . n      - second dimension of two dimensional array
    3365: . mstart - first index you will use in first coordinate direction (often 0)
    3366: - nstart - first index in the second coordinate direction (often 0)

    3368:   Output Parameter:
    3369: . a - location to put pointer to the array

    3371:   Level: developer

    3373:   Notes:
    3374:   For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
    3375:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
    3376:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
    3377:   the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.

    3379:   For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    3381: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
    3382:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
    3383:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    3384: @*/
    3385: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
    3386: {
    3387:   PetscInt           i, N;
    3388:   const PetscScalar *aa;

    3390:   PetscFunctionBegin;
    3392:   PetscAssertPointer(a, 6);
    3394:   PetscCall(VecGetLocalSize(x, &N));
    3395:   PetscCheck(m * n == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 2d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n);
    3396:   PetscCall(VecGetArrayRead(x, &aa));

    3398:   PetscCall(PetscMalloc1(m, a));
    3399:   for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
    3400:   *a -= mstart;
    3401:   PetscFunctionReturn(PETSC_SUCCESS);
    3402: }

    3404: /*@C
    3405:   VecRestoreArray2dRead - Restores a vector after `VecGetArray2dRead()` has been called.

    3407:   Logically Collective

    3409:   Input Parameters:
    3410: + x      - the vector
    3411: . m      - first dimension of two dimensional array
    3412: . n      - second dimension of the two dimensional array
    3413: . mstart - first index you will use in first coordinate direction (often 0)
    3414: . nstart - first index in the second coordinate direction (often 0)
    3415: - a      - location of pointer to array obtained from VecGetArray2d()

    3417:   Level: developer

    3419:   Notes:
    3420:   For regular PETSc vectors this routine does not involve any copies. For
    3421:   any special vectors that do not store local vector data in a contiguous
    3422:   array, this routine will copy the data back into the underlying
    3423:   vector data structure from the array obtained with `VecGetArray()`.

    3425:   This routine actually zeros out the `a` pointer.

    3427: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
    3428:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
    3429:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    3430: @*/
    3431: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
    3432: {
    3433:   void *dummy;

    3435:   PetscFunctionBegin;
    3437:   PetscAssertPointer(a, 6);
    3439:   dummy = (void *)(*a + mstart);
    3440:   PetscCall(PetscFree(dummy));
    3441:   PetscCall(VecRestoreArrayRead(x, NULL));
    3442:   *a = NULL;
    3443:   PetscFunctionReturn(PETSC_SUCCESS);
    3444: }

    3446: /*@C
    3447:   VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
    3448:   processor's portion of the vector data.  You MUST call `VecRestoreArray1dRead()`
    3449:   when you no longer need access to the array.

    3451:   Logically Collective

    3453:   Input Parameters:
    3454: + x      - the vector
    3455: . m      - first dimension of two dimensional array
    3456: - mstart - first index you will use in first coordinate direction (often 0)

    3458:   Output Parameter:
    3459: . a - location to put pointer to the array

    3461:   Level: developer

    3463:   Notes:
    3464:   For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
    3465:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
    3466:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

    3468:   For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    3470: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
    3471:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
    3472:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    3473: @*/
    3474: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
    3475: {
    3476:   PetscInt N;

    3478:   PetscFunctionBegin;
    3480:   PetscAssertPointer(a, 4);
    3482:   PetscCall(VecGetLocalSize(x, &N));
    3483:   PetscCheck(m == N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local array size %" PetscInt_FMT " does not match 1d array dimensions %" PetscInt_FMT, N, m);
    3484:   PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
    3485:   *a -= mstart;
    3486:   PetscFunctionReturn(PETSC_SUCCESS);
    3487: }

    3489: /*@C
    3490:   VecRestoreArray1dRead - Restores a vector after `VecGetArray1dRead()` has been called.

    3492:   Logically Collective

    3494:   Input Parameters:
    3495: + x      - the vector
    3496: . m      - first dimension of two dimensional array
    3497: . mstart - first index you will use in first coordinate direction (often 0)
    3498: - a      - location of pointer to array obtained from `VecGetArray1dRead()`

    3500:   Level: developer

    3502:   Notes:
    3503:   For regular PETSc vectors this routine does not involve any copies. For
    3504:   any special vectors that do not store local vector data in a contiguous
    3505:   array, this routine will copy the data back into the underlying
    3506:   vector data structure from the array obtained with `VecGetArray1dRead()`.

    3508:   This routine actually zeros out the `a` pointer.

    3510: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
    3511:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
    3512:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    3513: @*/
    3514: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
    3515: {
    3516:   PetscFunctionBegin;
    3519:   PetscCall(VecRestoreArrayRead(x, NULL));
    3520:   *a = NULL;
    3521:   PetscFunctionReturn(PETSC_SUCCESS);
    3522: }

    3524: /*@C
    3525:   VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
    3526:   processor's portion of the vector data.  You MUST call `VecRestoreArray3dRead()`
    3527:   when you no longer need access to the array.

    3529:   Logically Collective

    3531:   Input Parameters:
    3532: + x      - the vector
    3533: . m      - first dimension of three dimensional array
    3534: . n      - second dimension of three dimensional array
    3535: . p      - third dimension of three dimensional array
    3536: . mstart - first index you will use in first coordinate direction (often 0)
    3537: . nstart - first index in the second coordinate direction (often 0)
    3538: - pstart - first index in the third coordinate direction (often 0)

    3540:   Output Parameter:
    3541: . a - location to put pointer to the array

    3543:   Level: developer

    3545:   Notes:
    3546:   For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
    3547:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
    3548:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
    3549:   the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3dRead()`.

    3551:   For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    3553: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
    3554:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
    3555:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    3556: @*/
    3557: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
    3558: {
    3559:   PetscInt           i, N, j;
    3560:   const PetscScalar *aa;
    3561:   PetscScalar      **b;

    3563:   PetscFunctionBegin;
    3565:   PetscAssertPointer(a, 8);
    3567:   PetscCall(VecGetLocalSize(x, &N));
    3568:   PetscCheck(m * n * p == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 3d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n, p);
    3569:   PetscCall(VecGetArrayRead(x, &aa));

    3571:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
    3572:   b = (PetscScalar **)((*a) + m);
    3573:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
    3574:   for (i = 0; i < m; i++)
    3575:     for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset((PetscScalar *)aa, i * n * p + j * p - pstart);
    3576:   *a -= mstart;
    3577:   PetscFunctionReturn(PETSC_SUCCESS);
    3578: }

    3580: /*@C
    3581:   VecRestoreArray3dRead - Restores a vector after `VecGetArray3dRead()` has been called.

    3583:   Logically Collective

    3585:   Input Parameters:
    3586: + x      - the vector
    3587: . m      - first dimension of three dimensional array
    3588: . n      - second dimension of the three dimensional array
    3589: . p      - third dimension of the three dimensional array
    3590: . mstart - first index you will use in first coordinate direction (often 0)
    3591: . nstart - first index in the second coordinate direction (often 0)
    3592: . pstart - first index in the third coordinate direction (often 0)
    3593: - a      - location of pointer to array obtained from `VecGetArray3dRead()`

    3595:   Level: developer

    3597:   Notes:
    3598:   For regular PETSc vectors this routine does not involve any copies. For
    3599:   any special vectors that do not store local vector data in a contiguous
    3600:   array, this routine will copy the data back into the underlying
    3601:   vector data structure from the array obtained with `VecGetArray()`.

    3603:   This routine actually zeros out the `a` pointer.

    3605: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
    3606:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
    3607:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    3608: @*/
    3609: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
    3610: {
    3611:   void *dummy;

    3613:   PetscFunctionBegin;
    3615:   PetscAssertPointer(a, 8);
    3617:   dummy = (void *)(*a + mstart);
    3618:   PetscCall(PetscFree(dummy));
    3619:   PetscCall(VecRestoreArrayRead(x, NULL));
    3620:   *a = NULL;
    3621:   PetscFunctionReturn(PETSC_SUCCESS);
    3622: }

    3624: /*@C
    3625:   VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
    3626:   processor's portion of the vector data.  You MUST call `VecRestoreArray4dRead()`
    3627:   when you no longer need access to the array.

    3629:   Logically Collective

    3631:   Input Parameters:
    3632: + x      - the vector
    3633: . m      - first dimension of four dimensional array
    3634: . n      - second dimension of four dimensional array
    3635: . p      - third dimension of four dimensional array
    3636: . q      - fourth dimension of four dimensional array
    3637: . mstart - first index you will use in first coordinate direction (often 0)
    3638: . nstart - first index in the second coordinate direction (often 0)
    3639: . pstart - first index in the third coordinate direction (often 0)
    3640: - qstart - first index in the fourth coordinate direction (often 0)

    3642:   Output Parameter:
    3643: . a - location to put pointer to the array

    3645:   Level: beginner

    3647:   Notes:
    3648:   For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
    3649:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
    3650:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
    3651:   the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.

    3653:   For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    3655: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
    3656:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
    3657:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    3658: @*/
    3659: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
    3660: {
    3661:   PetscInt           i, N, j, k;
    3662:   const PetscScalar *aa;
    3663:   PetscScalar     ***b, **c;

    3665:   PetscFunctionBegin;
    3667:   PetscAssertPointer(a, 10);
    3669:   PetscCall(VecGetLocalSize(x, &N));
    3670:   PetscCheck(m * n * p * q == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 4d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n, p, q);
    3671:   PetscCall(VecGetArrayRead(x, &aa));

    3673:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
    3674:   b = (PetscScalar ***)((*a) + m);
    3675:   c = (PetscScalar **)(b + m * n);
    3676:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
    3677:   for (i = 0; i < m; i++)
    3678:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
    3679:   for (i = 0; i < m; i++)
    3680:     for (j = 0; j < n; j++)
    3681:       for (k = 0; k < p; k++) c[i * n * p + j * p + k] = (PetscScalar *)aa + i * n * p * q + j * p * q + k * q - qstart;
    3682:   *a -= mstart;
    3683:   PetscFunctionReturn(PETSC_SUCCESS);
    3684: }

    3686: /*@C
    3687:   VecRestoreArray4dRead - Restores a vector after `VecGetArray4d()` has been called.

    3689:   Logically Collective

    3691:   Input Parameters:
    3692: + x      - the vector
    3693: . m      - first dimension of four dimensional array
    3694: . n      - second dimension of the four dimensional array
    3695: . p      - third dimension of the four dimensional array
    3696: . q      - fourth dimension of the four dimensional array
    3697: . mstart - first index you will use in first coordinate direction (often 0)
    3698: . nstart - first index in the second coordinate direction (often 0)
    3699: . pstart - first index in the third coordinate direction (often 0)
    3700: . qstart - first index in the fourth coordinate direction (often 0)
    3701: - a      - location of pointer to array obtained from `VecGetArray4dRead()`

    3703:   Level: beginner

    3705:   Notes:
    3706:   For regular PETSc vectors this routine does not involve any copies. For
    3707:   any special vectors that do not store local vector data in a contiguous
    3708:   array, this routine will copy the data back into the underlying
    3709:   vector data structure from the array obtained with `VecGetArray()`.

    3711:   This routine actually zeros out the `a` pointer.

    3713: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
    3714:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
    3715:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
    3716: @*/
    3717: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
    3718: {
    3719:   void *dummy;

    3721:   PetscFunctionBegin;
    3723:   PetscAssertPointer(a, 10);
    3725:   dummy = (void *)(*a + mstart);
    3726:   PetscCall(PetscFree(dummy));
    3727:   PetscCall(VecRestoreArrayRead(x, NULL));
    3728:   *a = NULL;
    3729:   PetscFunctionReturn(PETSC_SUCCESS);
    3730: }

    3732: /*@
    3733:   VecLockGet - Get the current lock status of a vector

    3735:   Logically Collective

    3737:   Input Parameter:
    3738: . x - the vector

    3740:   Output Parameter:
    3741: . state - greater than zero indicates the vector is locked for read; less than zero indicates the vector is
    3742:            locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.

    3744:   Level: advanced

    3746: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
    3747: @*/
    3748: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
    3749: {
    3750:   PetscFunctionBegin;
    3752:   PetscAssertPointer(state, 2);
    3753:   *state = x->lock;
    3754:   PetscFunctionReturn(PETSC_SUCCESS);
    3755: }

    3757: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
    3758: {
    3759:   PetscFunctionBegin;
    3761:   PetscAssertPointer(file, 2);
    3762:   PetscAssertPointer(func, 3);
    3763:   PetscAssertPointer(line, 4);
    3764: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
    3765:   {
    3766:     const int index = x->lockstack.currentsize - 1;

    3768:     *file = index < 0 ? NULL : x->lockstack.file[index];
    3769:     *func = index < 0 ? NULL : x->lockstack.function[index];
    3770:     *line = index < 0 ? 0 : x->lockstack.line[index];
    3771:   }
    3772: #else
    3773:   *file = NULL;
    3774:   *func = NULL;
    3775:   *line = 0;
    3776: #endif
    3777:   PetscFunctionReturn(PETSC_SUCCESS);
    3778: }

    3780: /*@
    3781:   VecLockReadPush - Push a read-only lock on a vector to prevent it from being written to

    3783:   Logically Collective

    3785:   Input Parameter:
    3786: . x - the vector

    3788:   Level: intermediate

    3790:   Notes:
    3791:   If this is set then calls to `VecGetArray()` or `VecSetValues()` or any other routines that change the vectors values will generate an error.

    3793:   The call can be nested, i.e., called multiple times on the same vector, but each `VecLockReadPush()` has to have one matching
    3794:   `VecLockReadPop()`, which removes the latest read-only lock.

    3796: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
    3797: @*/
    3798: PetscErrorCode VecLockReadPush(Vec x)
    3799: {
    3800:   PetscFunctionBegin;
    3802:   PetscCheck(x->lock++ >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector is already locked for exclusive write access but you want to read it");
    3803: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
    3804:   {
    3805:     const char *file, *func;
    3806:     int         index, line;

    3808:     if ((index = petscstack.currentsize - 2) < 0) {
    3809:       // vec was locked "outside" of petsc, either in user-land or main. the error message will
    3810:       // now show this function as the culprit, but it will include the stacktrace
    3811:       file = "unknown user-file";
    3812:       func = "unknown_user_function";
    3813:       line = 0;
    3814:     } else {
    3815:       file = petscstack.file[index];
    3816:       func = petscstack.function[index];
    3817:       line = petscstack.line[index];
    3818:     }
    3819:     PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
    3820:   }
    3821: #endif
    3822:   PetscFunctionReturn(PETSC_SUCCESS);
    3823: }

    3825: /*@
    3826:   VecLockReadPop - Pop a read-only lock from a vector

    3828:   Logically Collective

    3830:   Input Parameter:
    3831: . x - the vector

    3833:   Level: intermediate

    3835: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
    3836: @*/
    3837: PetscErrorCode VecLockReadPop(Vec x)
    3838: {
    3839:   PetscFunctionBegin;
    3841:   PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
    3842: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
    3843:   {
    3844:     const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];

    3846:     PetscStackPop_Private(x->lockstack, previous);
    3847:   }
    3848: #endif
    3849:   PetscFunctionReturn(PETSC_SUCCESS);
    3850: }

    3852: /*@
    3853:   VecLockWriteSet - Lock or unlock a vector for exclusive read/write access

    3855:   Logically Collective

    3857:   Input Parameters:
    3858: + x   - the vector
    3859: - flg - `PETSC_TRUE` to lock the vector for exclusive read/write access; `PETSC_FALSE` to unlock it.

    3861:   Level: intermediate

    3863:   Notes:
    3864:   The function is useful in split-phase computations, which usually have a begin phase and an end phase.
    3865:   One can call `VecLockWriteSet`(x,`PETSC_TRUE`) in the begin phase to lock a vector for exclusive
    3866:   access, and call `VecLockWriteSet`(x,`PETSC_FALSE`) in the end phase to unlock the vector from exclusive
    3867:   access. In this way, one is ensured no other operations can access the vector in between. The code may like

    3869: .vb
    3870:        VecGetArray(x,&xdata); // begin phase
    3871:        VecLockWriteSet(v,PETSC_TRUE);

    3873:        Other operations, which can not access x anymore (they can access xdata, of course)

    3875:        VecRestoreArray(x,&vdata); // end phase
    3876:        VecLockWriteSet(v,PETSC_FALSE);
    3877: .ve

    3879:   The call can not be nested on the same vector, in other words, one can not call `VecLockWriteSet`(x,`PETSC_TRUE`)
    3880:   again before calling `VecLockWriteSet`(v,`PETSC_FALSE`).

    3882: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
    3883: @*/
    3884: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
    3885: {
    3886:   PetscFunctionBegin;
    3888:   if (flg) {
    3889:     PetscCheck(x->lock <= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector is already locked for read-only access but you want to write it");
    3890:     PetscCheck(x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector is already locked for exclusive write access but you want to write it");
    3891:     x->lock = -1;
    3892:   } else {
    3893:     PetscCheck(x->lock == -1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector is not locked for exclusive write access but you want to unlock it from that");
    3894:     x->lock = 0;
    3895:   }
    3896:   PetscFunctionReturn(PETSC_SUCCESS);
    3897: }