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/sfimpl.h"
  6: #include "petscsystypes.h"
  7: #include <petsc/private/vecimpl.h>

  9: PetscInt VecGetSubVectorSavedStateId = -1;

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

 21:     PetscCall(VecGetOffloadMask(vec, &mask));
 22:     if (!PetscOffloadHost(mask)) PetscFunctionReturn(PETSC_SUCCESS);
 23:     PetscCall(VecGetLocalSize(vec, &n));
 24:     PetscCall(VecGetArrayRead(vec, &x));
 25:     for (PetscInt i = 0; i < n; i++) {
 26:       if (begin) {
 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 beginning of function: Parameter number %" PetscInt_FMT, i, argnum);
 28:       } else {
 29:         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);
 30:       }
 31:     }
 32:     PetscCall(VecRestoreArrayRead(vec, &x));
 33:   }
 34:   PetscFunctionReturn(PETSC_SUCCESS);
 35: }
 36: #endif

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

 41:   Logically Collective

 43:   Input Parameters:
 44: + x - the numerators
 45: - y - the denominators

 47:   Output Parameter:
 48: . max - the result

 50:   Level: advanced

 52:   Notes:
 53:   `x` and `y` may be the same vector

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

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

 79: /*@
 80:   VecDot - Computes the vector dot product.

 82:   Collective

 84:   Input Parameters:
 85: + x - first vector
 86: - y - second vector

 88:   Output Parameter:
 89: . val - the dot product

 91:   Level: intermediate

 93:   Notes for Users of Complex Numbers:
 94:   For complex vectors, `VecDot()` computes
 95: $     val = (x,y) = y^H x,
 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: $     val = (x,y) = y^T x,
102:   where y^T denotes the transpose of y.

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

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

129: /*@
130:   VecDotRealPart - Computes the real part of the vector dot product.

132:   Collective

134:   Input Parameters:
135: + x - first vector
136: - y - second vector

138:   Output Parameter:
139: . val - the real part of the dot product;

141:   Level: intermediate

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

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

148:   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
149:   the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)

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

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

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

166: /*@
167:   VecNorm  - Computes the vector norm.

169:   Collective

171:   Input Parameters:
172: + x    - the vector
173: - type - the type of the norm requested

175:   Output Parameter:
176: . val - the norm

178:   Level: intermediate

180:   Notes:
181:   See `NormType` for descriptions of each norm.

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

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

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

201:   PetscFunctionBegin;
204:   VecCheckAssembled(x);
206:   PetscAssertPointer(val, 3);

208:   PetscCall(VecNormAvailable(x, type, &flg, val));
209:   // check that all MPI processes call this routine together and have same availability
210:   if (PetscDefined(USE_DEBUG)) {
211:     PetscBool minflg;

213:     PetscCall(MPIU_Allreduce(&flg, &minflg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)x)));
214:     PetscCheck(flg == minflg, PetscObjectComm((PetscObject)(x)), PETSC_ERR_ARG_WRONGSTATE, "Some MPI processes have cached norm, others do not. This may happen when some MPI processes call VecGetArray() and some others do not.");
215:   }
216:   if (flg) PetscFunctionReturn(PETSC_SUCCESS);

218:   PetscCall(VecLockReadPush(x));
219:   PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
220:   PetscUseTypeMethod(x, norm, type, val);
221:   PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
222:   PetscCall(VecLockReadPop(x));

224:   if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
225:   PetscFunctionReturn(PETSC_SUCCESS);
226: }

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

231:   Not Collective

233:   Input Parameters:
234: + x    - the vector
235: - 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
236:           `NORM_1_AND_2`, which computes both norms and stores them
237:           in a two element array.

239:   Output Parameters:
240: + available - `PETSC_TRUE` if the val returned is valid
241: - val       - the norm

243:   Level: intermediate

245:   Developer Notes:
246:   `PETSC_HAVE_SLOW_BLAS_NORM2` will cause a C (loop unrolled) version of the norm to be used, rather
247:   than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
248:   (as opposed to vendor provided) because the FORTRAN BLAS `NRM2()` routine is very slow.

250: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
251:           `VecNormBegin()`, `VecNormEnd()`
252: @*/
253: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
254: {
255:   PetscFunctionBegin;
258:   PetscAssertPointer(available, 3);
259:   PetscAssertPointer(val, 4);

261:   if (type == NORM_1_AND_2) {
262:     *available = PETSC_FALSE;
263:   } else {
264:     PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
265:   }
266:   PetscFunctionReturn(PETSC_SUCCESS);
267: }

269: /*@
270:   VecNormalize - Normalizes a vector by its 2-norm.

272:   Collective

274:   Input Parameter:
275: . x - the vector

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

280:   Level: intermediate

282: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
283: @*/
284: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
285: {
286:   PetscReal norm;

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

306: /*@C
307:   VecMax - Determines the vector component with maximum real part and its location.

309:   Collective

311:   Input Parameter:
312: . x - the vector

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

318:   Level: intermediate

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

323:   Returns the smallest index with the maximum value

325: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
326: @*/
327: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
328: {
329:   PetscFunctionBegin;
332:   VecCheckAssembled(x);
333:   if (p) PetscAssertPointer(p, 2);
334:   PetscAssertPointer(val, 3);
335:   PetscCall(VecLockReadPush(x));
336:   PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
337:   PetscUseTypeMethod(x, max, p, val);
338:   PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
339:   PetscCall(VecLockReadPop(x));
340:   PetscFunctionReturn(PETSC_SUCCESS);
341: }

343: /*@C
344:   VecMin - Determines the vector component with minimum real part and its location.

346:   Collective

348:   Input Parameter:
349: . x - the vector

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

355:   Level: intermediate

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

360:   This returns the smallest index with the minimum value

362: .seealso: [](ch_vectors), `Vec`, `VecMax()`
363: @*/
364: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
365: {
366:   PetscFunctionBegin;
369:   VecCheckAssembled(x);
370:   if (p) PetscAssertPointer(p, 2);
371:   PetscAssertPointer(val, 3);
372:   PetscCall(VecLockReadPush(x));
373:   PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
374:   PetscUseTypeMethod(x, min, p, val);
375:   PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
376:   PetscCall(VecLockReadPop(x));
377:   PetscFunctionReturn(PETSC_SUCCESS);
378: }

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

384:   Collective

386:   Input Parameters:
387: + x - first vector
388: - y - second vector

390:   Output Parameter:
391: . val - the dot product

393:   Level: intermediate

395:   Notes for Users of Complex Numbers:
396:   For complex vectors, `VecTDot()` computes the indefinite form
397: $     val = (x,y) = y^T x,
398:   where y^T denotes the transpose of y.

400:   Use `VecDot()` for the inner product
401: $     val = (x,y) = y^H x,
402:   where y^H denotes the conjugate transpose of y.

404: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
405: @*/
406: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
407: {
408:   PetscFunctionBegin;
411:   PetscAssertPointer(val, 3);
414:   PetscCheckSameTypeAndComm(x, 1, y, 2);
415:   VecCheckSameSize(x, 1, y, 2);
416:   VecCheckAssembled(x);
417:   VecCheckAssembled(y);

419:   PetscCall(VecLockReadPush(x));
420:   PetscCall(VecLockReadPush(y));
421:   PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
422:   PetscUseTypeMethod(x, tdot, y, val);
423:   PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
424:   PetscCall(VecLockReadPop(x));
425:   PetscCall(VecLockReadPop(y));
426:   PetscFunctionReturn(PETSC_SUCCESS);
427: }

429: PetscErrorCode VecScaleAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
430: {
431:   PetscReal   norms[4];
432:   PetscBool   flgs[4];
433:   PetscScalar one = 1.0;

435:   PetscFunctionBegin;
438:   VecCheckAssembled(x);
439:   PetscCall(VecSetErrorIfLocked(x, 1));
440:   if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);

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

445:   PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
446:   VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
447:   PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));

449:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
450:   /* put the scaled stashed norms back into the Vec */
451:   for (PetscInt i = 0; i < 4; i++) {
452:     PetscReal ar = PetscAbsScalar(alpha);
453:     if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
454:   }
455:   PetscFunctionReturn(PETSC_SUCCESS);
456: }

458: /*@
459:   VecScale - Scales a vector.

461:   Not Collective

463:   Input Parameters:
464: + x     - the vector
465: - alpha - the scalar

467:   Level: intermediate

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

472: .seealso: [](ch_vectors), `Vec`, `VecSet()`
473: @*/
474: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
475: {
476:   PetscFunctionBegin;
477:   PetscCall(VecScaleAsync_Private(x, alpha, NULL));
478:   PetscFunctionReturn(PETSC_SUCCESS);
479: }

481: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
482: {
483:   PetscFunctionBegin;
486:   VecCheckAssembled(x);
488:   PetscCall(VecSetErrorIfLocked(x, 1));

490:   if (alpha == 0) {
491:     PetscReal norm;
492:     PetscBool set;

494:     PetscCall(VecNormAvailable(x, NORM_2, &set, &norm));
495:     if (set == PETSC_TRUE && norm == 0) PetscFunctionReturn(PETSC_SUCCESS);
496:   }
497:   PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
498:   VecMethodDispatch(x, dctx, VecAsyncFnName(Set), set, (Vec, PetscScalar, PetscDeviceContext), alpha);
499:   PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
500:   PetscCall(PetscObjectStateIncrease((PetscObject)x));

502:   /*  norms can be simply set (if |alpha|*N not too large) */

504:   {
505:     PetscReal      val = PetscAbsScalar(alpha);
506:     const PetscInt N   = x->map->N;

508:     if (N == 0) {
509:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0l));
510:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
511:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
512:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
513:     } else if (val > PETSC_MAX_REAL / N) {
514:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
515:     } else {
516:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
517:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
518:       val *= PetscSqrtReal((PetscReal)N);
519:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
520:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
521:     }
522:   }
523:   PetscFunctionReturn(PETSC_SUCCESS);
524: }

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

529:   Logically Collective

531:   Input Parameters:
532: + x     - the vector
533: - alpha - the scalar

535:   Level: beginner

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

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

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

548: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
549: @*/
550: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
551: {
552:   PetscFunctionBegin;
553:   PetscCall(VecSetAsync_Private(x, alpha, NULL));
554:   PetscFunctionReturn(PETSC_SUCCESS);
555: }

557: PetscErrorCode VecAXPYAsync_Private(Vec y, PetscScalar alpha, Vec x, PetscDeviceContext dctx)
558: {
559:   PetscFunctionBegin;
564:   PetscCheckSameTypeAndComm(x, 3, y, 1);
565:   VecCheckSameSize(x, 3, y, 1);
566:   VecCheckAssembled(x);
567:   VecCheckAssembled(y);
569:   if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
570:   PetscCall(VecSetErrorIfLocked(y, 1));
571:   if (x == y) {
572:     PetscCall(VecScale(y, alpha + 1.0));
573:     PetscFunctionReturn(PETSC_SUCCESS);
574:   }
575:   PetscCall(VecLockReadPush(x));
576:   PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
577:   VecMethodDispatch(y, dctx, VecAsyncFnName(AXPY), axpy, (Vec, PetscScalar, Vec, PetscDeviceContext), alpha, x);
578:   PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
579:   PetscCall(VecLockReadPop(x));
580:   PetscCall(PetscObjectStateIncrease((PetscObject)y));
581:   PetscFunctionReturn(PETSC_SUCCESS);
582: }
583: /*@
584:   VecAXPY - Computes `y = alpha x + y`.

586:   Logically Collective

588:   Input Parameters:
589: + alpha - the scalar
590: . x     - vector scale by `alpha`
591: - y     - vector accumulated into

593:   Output Parameter:
594: . y - output vector

596:   Level: intermediate

598:   Notes:
599:   This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
600: .vb
601:     VecAXPY(y,alpha,x)                   y = alpha x           +      y
602:     VecAYPX(y,beta,x)                    y =       x           + beta y
603:     VecAXPBY(y,alpha,beta,x)             y = alpha x           + beta y
604:     VecWAXPY(w,alpha,x,y)                w = alpha x           +      y
605:     VecAXPBYPCZ(z,alpha,beta,gamma,x,y)  z = alpha x           + beta y + gamma z
606:     VecMAXPY(y,nv,alpha[],x[])           y = sum alpha[i] x[i] +      y
607: .ve

609: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
610: @*/
611: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
612: {
613:   PetscFunctionBegin;
614:   PetscCall(VecAXPYAsync_Private(y, alpha, x, NULL));
615:   PetscFunctionReturn(PETSC_SUCCESS);
616: }

618: PetscErrorCode VecAYPXAsync_Private(Vec y, PetscScalar beta, Vec x, PetscDeviceContext dctx)
619: {
620:   PetscFunctionBegin;
625:   PetscCheckSameTypeAndComm(x, 3, y, 1);
626:   VecCheckSameSize(x, 1, y, 3);
627:   VecCheckAssembled(x);
628:   VecCheckAssembled(y);
630:   PetscCall(VecSetErrorIfLocked(y, 1));
631:   if (x == y) {
632:     PetscCall(VecScale(y, beta + 1.0));
633:     PetscFunctionReturn(PETSC_SUCCESS);
634:   }
635:   PetscCall(VecLockReadPush(x));
636:   if (beta == (PetscScalar)0.0) {
637:     PetscCall(VecCopy(x, y));
638:   } else {
639:     PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
640:     VecMethodDispatch(y, dctx, VecAsyncFnName(AYPX), aypx, (Vec, PetscScalar, Vec, PetscDeviceContext), beta, x);
641:     PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
642:     PetscCall(PetscObjectStateIncrease((PetscObject)y));
643:   }
644:   PetscCall(VecLockReadPop(x));
645:   PetscFunctionReturn(PETSC_SUCCESS);
646: }

648: /*@
649:   VecAYPX - Computes `y = x + beta y`.

651:   Logically Collective

653:   Input Parameters:
654: + beta - the scalar
655: . x    - the unscaled vector
656: - y    - the vector to be scaled

658:   Output Parameter:
659: . y - output vector

661:   Level: intermediate

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

666: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
667: @*/
668: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
669: {
670:   PetscFunctionBegin;
671:   PetscCall(VecAYPXAsync_Private(y, beta, x, NULL));
672:   PetscFunctionReturn(PETSC_SUCCESS);
673: }

675: PetscErrorCode VecAXPBYAsync_Private(Vec y, PetscScalar alpha, PetscScalar beta, Vec x, PetscDeviceContext dctx)
676: {
677:   PetscFunctionBegin;
682:   PetscCheckSameTypeAndComm(x, 4, y, 1);
683:   VecCheckSameSize(y, 1, x, 4);
684:   VecCheckAssembled(x);
685:   VecCheckAssembled(y);
688:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
689:   if (x == y) {
690:     PetscCall(VecScale(y, alpha + beta));
691:     PetscFunctionReturn(PETSC_SUCCESS);
692:   }

694:   PetscCall(VecSetErrorIfLocked(y, 1));
695:   PetscCall(VecLockReadPush(x));
696:   PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
697:   VecMethodDispatch(y, dctx, VecAsyncFnName(AXPBY), axpby, (Vec, PetscScalar, PetscScalar, Vec, PetscDeviceContext), alpha, beta, x);
698:   PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
699:   PetscCall(PetscObjectStateIncrease((PetscObject)y));
700:   PetscCall(VecLockReadPop(x));
701:   PetscFunctionReturn(PETSC_SUCCESS);
702: }
703: /*@
704:   VecAXPBY - Computes `y = alpha x + beta y`.

706:   Logically Collective

708:   Input Parameters:
709: + alpha - first scalar
710: . beta  - second scalar
711: . x     - the first scaled vector
712: - y     - the second scaled vector

714:   Output Parameter:
715: . y - output vector

717:   Level: intermediate

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

722: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
723: @*/
724: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
725: {
726:   PetscFunctionBegin;
727:   PetscCall(VecAXPBYAsync_Private(y, alpha, beta, x, NULL));
728:   PetscFunctionReturn(PETSC_SUCCESS);
729: }

731: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
732: {
733:   PetscFunctionBegin;
740:   PetscCheckSameTypeAndComm(x, 5, y, 6);
741:   PetscCheckSameTypeAndComm(x, 5, z, 1);
742:   VecCheckSameSize(x, 1, y, 5);
743:   VecCheckSameSize(x, 1, z, 6);
744:   PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
745:   PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
746:   VecCheckAssembled(x);
747:   VecCheckAssembled(y);
748:   VecCheckAssembled(z);
752:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);

754:   PetscCall(VecSetErrorIfLocked(z, 1));
755:   PetscCall(VecLockReadPush(x));
756:   PetscCall(VecLockReadPush(y));
757:   PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
758:   VecMethodDispatch(z, dctx, VecAsyncFnName(AXPBYPCZ), axpbypcz, (Vec, PetscScalar, PetscScalar, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, beta, gamma, x, y);
759:   PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
760:   PetscCall(PetscObjectStateIncrease((PetscObject)z));
761:   PetscCall(VecLockReadPop(x));
762:   PetscCall(VecLockReadPop(y));
763:   PetscFunctionReturn(PETSC_SUCCESS);
764: }
765: /*@
766:   VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`

768:   Logically Collective

770:   Input Parameters:
771: + alpha - first scalar
772: . beta  - second scalar
773: . gamma - third scalar
774: . x     - first vector
775: . y     - second vector
776: - z     - third vector

778:   Output Parameter:
779: . z - output vector

781:   Level: intermediate

783:   Note:
784:   `x`, `y` and `z` must be different vectors

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

789: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
790: @*/
791: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
792: {
793:   PetscFunctionBegin;
794:   PetscCall(VecAXPBYPCZAsync_Private(z, alpha, beta, gamma, x, y, NULL));
795:   PetscFunctionReturn(PETSC_SUCCESS);
796: }

798: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
799: {
800:   PetscFunctionBegin;
807:   PetscCheckSameTypeAndComm(x, 3, y, 4);
808:   PetscCheckSameTypeAndComm(y, 4, w, 1);
809:   VecCheckSameSize(x, 3, y, 4);
810:   VecCheckSameSize(x, 3, w, 1);
811:   PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
812:   PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
813:   VecCheckAssembled(x);
814:   VecCheckAssembled(y);
816:   PetscCall(VecSetErrorIfLocked(w, 1));

818:   PetscCall(VecLockReadPush(x));
819:   PetscCall(VecLockReadPush(y));
820:   if (alpha == (PetscScalar)0.0) {
821:     PetscCall(VecCopyAsync_Private(y, w, dctx));
822:   } else {
823:     PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
824:     VecMethodDispatch(w, dctx, VecAsyncFnName(WAXPY), waxpy, (Vec, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, x, y);
825:     PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
826:     PetscCall(PetscObjectStateIncrease((PetscObject)w));
827:   }
828:   PetscCall(VecLockReadPop(x));
829:   PetscCall(VecLockReadPop(y));
830:   PetscFunctionReturn(PETSC_SUCCESS);
831: }

833: /*@
834:   VecWAXPY - Computes `w = alpha x + y`.

836:   Logically Collective

838:   Input Parameters:
839: + alpha - the scalar
840: . x     - first vector, multiplied by `alpha`
841: - y     - second vector

843:   Output Parameter:
844: . w - the result

846:   Level: intermediate

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

851:   Developer Notes:
852:   The implementation is optimized for alpha of -1.0, 0.0, and 1.0

854: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
855: @*/
856: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
857: {
858:   PetscFunctionBegin;
859:   PetscCall(VecWAXPYAsync_Private(w, alpha, x, y, NULL));
860:   PetscFunctionReturn(PETSC_SUCCESS);
861: }

863: /*@C
864:   VecSetValues - Inserts or adds values into certain locations of a vector.

866:   Not Collective

868:   Input Parameters:
869: + x    - vector to insert in
870: . ni   - number of elements to add
871: . ix   - indices where to add
872: . y    - array of values
873: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries

875:   Level: beginner

877:   Notes:
878: .vb
879:    `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
880: .ve

882:   Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
883:   options cannot be mixed without intervening calls to the assembly
884:   routines.

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

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

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

897: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
898:           `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
899: @*/
900: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
901: {
902:   PetscFunctionBeginHot;
904:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
905:   PetscAssertPointer(ix, 3);
906:   PetscAssertPointer(y, 4);

909:   PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
910:   PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
911:   PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
912:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
913:   PetscFunctionReturn(PETSC_SUCCESS);
914: }

916: /*@C
917:   VecGetValues - Gets values from certain locations of a vector. Currently
918:   can only get values on the same processor on which they are owned

920:   Not Collective

922:   Input Parameters:
923: + x  - vector to get values from
924: . ni - number of elements to get
925: - ix - indices where to get them from (in global 1d numbering)

927:   Output Parameter:
928: . y - array of values

930:   Level: beginner

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

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

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

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

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

945: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
946: @*/
947: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
948: {
949:   PetscFunctionBegin;
951:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
952:   PetscAssertPointer(ix, 3);
953:   PetscAssertPointer(y, 4);
955:   VecCheckAssembled(x);
956:   PetscUseTypeMethod(x, getvalues, ni, ix, y);
957:   PetscFunctionReturn(PETSC_SUCCESS);
958: }

960: /*@C
961:   VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.

963:   Not Collective

965:   Input Parameters:
966: + x    - vector to insert in
967: . ni   - number of blocks to add
968: . ix   - indices where to add in block count, rather than element count
969: . y    - array of values
970: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries

972:   Level: intermediate

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

978:   Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
979:   options cannot be mixed without intervening calls to the assembly
980:   routines.

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

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

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

992: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
993:           `VecSetValues()`
994: @*/
995: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
996: {
997:   PetscFunctionBeginHot;
999:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1000:   PetscAssertPointer(ix, 3);
1001:   PetscAssertPointer(y, 4);

1004:   PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1005:   PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1006:   PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1007:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
1008:   PetscFunctionReturn(PETSC_SUCCESS);
1009: }

1011: /*@C
1012:   VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1013:   using a local ordering of the nodes.

1015:   Not Collective

1017:   Input Parameters:
1018: + x    - vector to insert in
1019: . ni   - number of elements to add
1020: . ix   - indices where to add
1021: . y    - array of values
1022: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries

1024:   Level: intermediate

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

1029:   Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
1030:   options cannot be mixed without intervening calls to the assembly
1031:   routines.

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

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

1038: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1039:           `VecSetValuesBlockedLocal()`
1040: @*/
1041: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1042: {
1043:   PetscInt lixp[128], *lix = lixp;

1045:   PetscFunctionBeginHot;
1047:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1048:   PetscAssertPointer(ix, 3);
1049:   PetscAssertPointer(y, 4);

1052:   PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1053:   if (!x->ops->setvalueslocal) {
1054:     if (x->map->mapping) {
1055:       if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1056:       PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1057:       PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1058:       if (ni > 128) PetscCall(PetscFree(lix));
1059:     } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1060:   } else PetscUseTypeMethod(x, setvalueslocal, ni, ix, y, iora);
1061:   PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1062:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
1063:   PetscFunctionReturn(PETSC_SUCCESS);
1064: }

1066: /*@
1067:   VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1068:   using a local ordering of the nodes.

1070:   Not Collective

1072:   Input Parameters:
1073: + x    - vector to insert in
1074: . ni   - number of blocks to add
1075: . ix   - indices where to add in block count, not element count
1076: . y    - array of values
1077: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries

1079:   Level: intermediate

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

1085:   Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1086:   options cannot be mixed without intervening calls to the assembly
1087:   routines.

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

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

1094: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1095:           `VecSetLocalToGlobalMapping()`
1096: @*/
1097: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1098: {
1099:   PetscInt lixp[128], *lix = lixp;

1101:   PetscFunctionBeginHot;
1103:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1104:   PetscAssertPointer(ix, 3);
1105:   PetscAssertPointer(y, 4);
1107:   PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1108:   if (x->map->mapping) {
1109:     if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1110:     PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1111:     PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1112:     if (ni > 128) PetscCall(PetscFree(lix));
1113:   } else {
1114:     PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1115:   }
1116:   PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1117:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
1118:   PetscFunctionReturn(PETSC_SUCCESS);
1119: }

1121: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1122: {
1123:   PetscFunctionBegin;
1126:   VecCheckAssembled(x);
1128:   if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1129:   PetscAssertPointer(y, 3);
1130:   for (PetscInt i = 0; i < nv; ++i) {
1133:     PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1134:     VecCheckSameSize(x, 1, y[i], 3);
1135:     VecCheckAssembled(y[i]);
1136:     PetscCall(VecLockReadPush(y[i]));
1137:   }
1138:   PetscAssertPointer(result, 4);

1141:   PetscCall(VecLockReadPush(x));
1142:   PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1143:   PetscCall((*mxdot)(x, nv, y, result));
1144:   PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1145:   PetscCall(VecLockReadPop(x));
1146:   for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1147:   PetscFunctionReturn(PETSC_SUCCESS);
1148: }

1150: /*@
1151:   VecMTDot - Computes indefinite vector multiple dot products.
1152:   That is, it does NOT use the complex conjugate.

1154:   Collective

1156:   Input Parameters:
1157: + x  - one vector
1158: . nv - number of vectors
1159: - y  - array of vectors.  Note that vectors are pointers

1161:   Output Parameter:
1162: . val - array of the dot products

1164:   Level: intermediate

1166:   Notes for Users of Complex Numbers:
1167:   For complex vectors, `VecMTDot()` computes the indefinite form
1168: $      val = (x,y) = y^T x,
1169:   where y^T denotes the transpose of y.

1171:   Use `VecMDot()` for the inner product
1172: $      val = (x,y) = y^H x,
1173:   where y^H denotes the conjugate transpose of y.

1175: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1176: @*/
1177: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1178: {
1179:   PetscFunctionBegin;
1181:   PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1182:   PetscFunctionReturn(PETSC_SUCCESS);
1183: }

1185: /*@
1186:   VecMDot - Computes multiple vector dot products.

1188:   Collective

1190:   Input Parameters:
1191: + x  - one vector
1192: . nv - number of vectors
1193: - y  - array of vectors.

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

1198:   Level: intermediate

1200:   Notes for Users of Complex Numbers:
1201:   For complex vectors, `VecMDot()` computes
1202: $     val = (x,y) = y^H x,
1203:   where y^H denotes the conjugate transpose of y.

1205:   Use `VecMTDot()` for the indefinite form
1206: $     val = (x,y) = y^T x,
1207:   where y^T denotes the transpose of y.

1209: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`
1210: @*/
1211: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1212: {
1213:   PetscFunctionBegin;
1215:   PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1216:   PetscFunctionReturn(PETSC_SUCCESS);
1217: }

1219: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
1220: {
1221:   PetscFunctionBegin;
1223:   VecCheckAssembled(y);
1225:   PetscCall(VecSetErrorIfLocked(y, 1));
1226:   PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1227:   if (nv) {
1228:     PetscInt zeros = 0;

1230:     PetscAssertPointer(alpha, 3);
1231:     PetscAssertPointer(x, 4);
1232:     for (PetscInt i = 0; i < nv; ++i) {
1236:       PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1237:       VecCheckSameSize(y, 1, x[i], 4);
1238:       PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1239:       VecCheckAssembled(x[i]);
1240:       PetscCall(VecLockReadPush(x[i]));
1241:       zeros += alpha[i] == (PetscScalar)0.0;
1242:     }

1244:     if (zeros < nv) {
1245:       PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1246:       VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1247:       PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1248:       PetscCall(PetscObjectStateIncrease((PetscObject)y));
1249:     }

1251:     for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1252:   }
1253:   PetscFunctionReturn(PETSC_SUCCESS);
1254: }

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

1259:   Logically Collective

1261:   Input Parameters:
1262: + nv    - number of scalars and x-vectors
1263: . alpha - array of scalars
1264: . y     - one vector
1265: - x     - array of vectors

1267:   Level: intermediate

1269:   Note:
1270:   `y` cannot be any of the `x` vectors

1272: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`,`VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1273: @*/
1274: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1275: {
1276:   PetscFunctionBegin;
1277:   PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1278:   PetscFunctionReturn(PETSC_SUCCESS);
1279: }

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

1284:   Logically Collective

1286:   Input Parameters:
1287: + nv    - number of scalars and x-vectors
1288: . alpha - array of scalars
1289: . beta  - scalar
1290: . y     - one vector
1291: - x     - array of vectors

1293:   Level: intermediate

1295:   Note:
1296:   `y` cannot be any of the `x` vectors.

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

1301: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1302: @*/
1303: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
1304: {
1305:   PetscFunctionBegin;
1307:   VecCheckAssembled(y);
1309:   PetscCall(VecSetErrorIfLocked(y, 1));
1310:   PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);

1313:   if (y->ops->maxpby) {
1314:     PetscInt zeros = 0;

1316:     if (nv) {
1317:       PetscAssertPointer(alpha, 3);
1318:       PetscAssertPointer(x, 5);
1319:     }

1321:     for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
1325:       PetscCheckSameTypeAndComm(y, 1, x[i], 5);
1326:       VecCheckSameSize(y, 1, x[i], 5);
1327:       PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1328:       VecCheckAssembled(x[i]);
1329:       PetscCall(VecLockReadPush(x[i]));
1330:       zeros += alpha[i] == (PetscScalar)0.0;
1331:     }

1333:     if (zeros < nv) { // has nonzero alpha
1334:       PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1335:       PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1336:       PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1337:       PetscCall(PetscObjectStateIncrease((PetscObject)y));
1338:     } else {
1339:       PetscCall(VecScale(y, beta));
1340:     }

1342:     for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1343:   } else { // no maxpby
1344:     if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1345:     else PetscCall(VecScale(y, beta));
1346:     PetscCall(VecMAXPY(y, nv, alpha, x));
1347:   }
1348:   PetscFunctionReturn(PETSC_SUCCESS);
1349: }

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

1356:   Collective

1358:   Input Parameters:
1359: + nx - number of vectors to be concatenated
1360: - X  - array containing the vectors to be concatenated in the order of concatenation

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

1366:   Level: advanced

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

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

1378: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1379: @*/
1380: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1381: {
1382:   MPI_Comm comm;
1383:   VecType  vec_type;
1384:   Vec      Ytmp, Xtmp;
1385:   IS      *is_tmp;
1386:   PetscInt i, shift = 0, Xnl, Xng, Xbegin;

1388:   PetscFunctionBegin;
1392:   PetscAssertPointer(Y, 3);

1394:   if ((*X)->ops->concatenate) {
1395:     /* use the dedicated concatenation function if available */
1396:     PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1397:   } else {
1398:     /* loop over vectors and start creating IS */
1399:     comm = PetscObjectComm((PetscObject)(*X));
1400:     PetscCall(VecGetType(*X, &vec_type));
1401:     PetscCall(PetscMalloc1(nx, &is_tmp));
1402:     for (i = 0; i < nx; i++) {
1403:       PetscCall(VecGetSize(X[i], &Xng));
1404:       PetscCall(VecGetLocalSize(X[i], &Xnl));
1405:       PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1406:       PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1407:       shift += Xng;
1408:     }
1409:     /* create the concatenated vector */
1410:     PetscCall(VecCreate(comm, &Ytmp));
1411:     PetscCall(VecSetType(Ytmp, vec_type));
1412:     PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1413:     PetscCall(VecSetUp(Ytmp));
1414:     /* copy data from X array to Y and return */
1415:     for (i = 0; i < nx; i++) {
1416:       PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1417:       PetscCall(VecCopy(X[i], Xtmp));
1418:       PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1419:     }
1420:     *Y = Ytmp;
1421:     if (x_is) {
1422:       *x_is = is_tmp;
1423:     } else {
1424:       for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1425:       PetscCall(PetscFree(is_tmp));
1426:     }
1427:   }
1428:   PetscFunctionReturn(PETSC_SUCCESS);
1429: }

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

1434:     Input Parameters:
1435: +   X - the original vector
1436: -   is - the index set of the subvector

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

1443: */
1444: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1445: {
1446:   PetscInt  gstart, gend, lstart;
1447:   PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1448:   PetscInt  n, N, ibs, vbs, bs = -1;

1450:   PetscFunctionBegin;
1451:   PetscCall(ISGetLocalSize(is, &n));
1452:   PetscCall(ISGetSize(is, &N));
1453:   PetscCall(ISGetBlockSize(is, &ibs));
1454:   PetscCall(VecGetBlockSize(X, &vbs));
1455:   PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1456:   PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1457:   /* block size is given by IS if ibs > 1; otherwise, check the vector */
1458:   if (ibs > 1) {
1459:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1460:     bs = ibs;
1461:   } else {
1462:     if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1463:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1464:     if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1465:   }

1467:   *contig    = red[0];
1468:   *start     = lstart;
1469:   *blocksize = bs;
1470:   PetscFunctionReturn(PETSC_SUCCESS);
1471: }

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

1475:     Input Parameters:
1476: +   X - the original vector
1477: .   is - the index set of the subvector
1478: -   bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()

1480:     Output Parameter:
1481: .   Z  - the subvector, which will compose the VecScatter context on output
1482: */
1483: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1484: {
1485:   PetscInt   n, N;
1486:   VecScatter vscat;
1487:   Vec        Y;

1489:   PetscFunctionBegin;
1490:   PetscCall(ISGetLocalSize(is, &n));
1491:   PetscCall(ISGetSize(is, &N));
1492:   PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1493:   PetscCall(VecSetSizes(Y, n, N));
1494:   PetscCall(VecSetBlockSize(Y, bs));
1495:   PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1496:   PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1497:   PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1498:   PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1499:   PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1500:   PetscCall(VecScatterDestroy(&vscat));
1501:   *Z = Y;
1502:   PetscFunctionReturn(PETSC_SUCCESS);
1503: }

1505: /*@
1506:   VecGetSubVector - Gets a vector representing part of another vector

1508:   Collective

1510:   Input Parameters:
1511: + X  - vector from which to extract a subvector
1512: - is - index set representing portion of `X` to extract

1514:   Output Parameter:
1515: . Y - subvector corresponding to `is`

1517:   Level: advanced

1519:   Notes:
1520:   The subvector `Y` should be returned with `VecRestoreSubVector()`.
1521:   `X` and must be defined on the same communicator

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

1526:   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`.

1528: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1529: @*/
1530: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1531: {
1532:   Vec Z;

1534:   PetscFunctionBegin;
1537:   PetscCheckSameComm(X, 1, is, 2);
1538:   PetscAssertPointer(Y, 3);
1539:   if (X->ops->getsubvector) {
1540:     PetscUseTypeMethod(X, getsubvector, is, &Z);
1541:   } else { /* Default implementation currently does no caching */
1542:     PetscBool contig;
1543:     PetscInt  n, N, start, bs;

1545:     PetscCall(ISGetLocalSize(is, &n));
1546:     PetscCall(ISGetSize(is, &N));
1547:     PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1548:     if (contig) { /* We can do a no-copy implementation */
1549:       const PetscScalar *x;
1550:       PetscInt           state = 0;
1551:       PetscBool          isstd, iscuda, iship;

1553:       PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1554:       PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1555:       PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1556:       if (iscuda) {
1557: #if defined(PETSC_HAVE_CUDA)
1558:         const PetscScalar *x_d;
1559:         PetscMPIInt        size;
1560:         PetscOffloadMask   flg;

1562:         PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1563:         PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1564:         PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1565:         if (x) x += start;
1566:         if (x_d) x_d += start;
1567:         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1568:         if (size == 1) {
1569:           PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1570:         } else {
1571:           PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1572:         }
1573:         Z->offloadmask = flg;
1574: #endif
1575:       } else if (iship) {
1576: #if defined(PETSC_HAVE_HIP)
1577:         const PetscScalar *x_d;
1578:         PetscMPIInt        size;
1579:         PetscOffloadMask   flg;

1581:         PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1582:         PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1583:         PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1584:         if (x) x += start;
1585:         if (x_d) x_d += start;
1586:         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1587:         if (size == 1) {
1588:           PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1589:         } else {
1590:           PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1591:         }
1592:         Z->offloadmask = flg;
1593: #endif
1594:       } else if (isstd) {
1595:         PetscMPIInt size;

1597:         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1598:         PetscCall(VecGetArrayRead(X, &x));
1599:         if (x) x += start;
1600:         if (size == 1) {
1601:           PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1602:         } else {
1603:           PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1604:         }
1605:         PetscCall(VecRestoreArrayRead(X, &x));
1606:       } else { /* default implementation: use place array */
1607:         PetscCall(VecGetArrayRead(X, &x));
1608:         PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1609:         PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1610:         PetscCall(VecSetSizes(Z, n, N));
1611:         PetscCall(VecSetBlockSize(Z, bs));
1612:         PetscCall(VecPlaceArray(Z, x ? x + start : NULL));
1613:         PetscCall(VecRestoreArrayRead(X, &x));
1614:       }

1616:       /* this is relevant only in debug mode */
1617:       PetscCall(VecLockGet(X, &state));
1618:       if (state) PetscCall(VecLockReadPush(Z));
1619:       Z->ops->placearray   = NULL;
1620:       Z->ops->replacearray = NULL;
1621:     } else { /* Have to create a scatter and do a copy */
1622:       PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1623:     }
1624:   }
1625:   /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1626:   if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1627:   PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1628:   *Y = Z;
1629:   PetscFunctionReturn(PETSC_SUCCESS);
1630: }

1632: /*@
1633:   VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`

1635:   Collective

1637:   Input Parameters:
1638: + X  - vector from which subvector was obtained
1639: . is - index set representing the subset of `X`
1640: - Y  - subvector being restored

1642:   Level: advanced

1644: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1645: @*/
1646: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1647: {
1648:   PETSC_UNUSED PetscObjectState dummystate = 0;
1649:   PetscBool                     unchanged;

1651:   PetscFunctionBegin;
1654:   PetscCheckSameComm(X, 1, is, 2);
1655:   PetscAssertPointer(Y, 3);

1658:   if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1659:   else {
1660:     PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1661:     if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1662:       VecScatter scatter;
1663:       PetscInt   state;

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

1668:       PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1669:       if (scatter) {
1670:         PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1671:         PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1672:       } else {
1673:         PetscBool iscuda, iship;
1674:         PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1675:         PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));

1677:         if (iscuda) {
1678: #if defined(PETSC_HAVE_CUDA)
1679:           PetscOffloadMask ymask = (*Y)->offloadmask;

1681:           /* The offloadmask of X dictates where to move memory
1682:               If X GPU data is valid, then move Y data on GPU if needed
1683:               Otherwise, move back to the CPU */
1684:           switch (X->offloadmask) {
1685:           case PETSC_OFFLOAD_BOTH:
1686:             if (ymask == PETSC_OFFLOAD_CPU) {
1687:               PetscCall(VecCUDAResetArray(*Y));
1688:             } else if (ymask == PETSC_OFFLOAD_GPU) {
1689:               X->offloadmask = PETSC_OFFLOAD_GPU;
1690:             }
1691:             break;
1692:           case PETSC_OFFLOAD_GPU:
1693:             if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1694:             break;
1695:           case PETSC_OFFLOAD_CPU:
1696:             if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1697:             break;
1698:           case PETSC_OFFLOAD_UNALLOCATED:
1699:           case PETSC_OFFLOAD_KOKKOS:
1700:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1701:           }
1702: #endif
1703:         } else if (iship) {
1704: #if defined(PETSC_HAVE_HIP)
1705:           PetscOffloadMask ymask = (*Y)->offloadmask;

1707:           /* The offloadmask of X dictates where to move memory
1708:               If X GPU data is valid, then move Y data on GPU if needed
1709:               Otherwise, move back to the CPU */
1710:           switch (X->offloadmask) {
1711:           case PETSC_OFFLOAD_BOTH:
1712:             if (ymask == PETSC_OFFLOAD_CPU) {
1713:               PetscCall(VecHIPResetArray(*Y));
1714:             } else if (ymask == PETSC_OFFLOAD_GPU) {
1715:               X->offloadmask = PETSC_OFFLOAD_GPU;
1716:             }
1717:             break;
1718:           case PETSC_OFFLOAD_GPU:
1719:             if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1720:             break;
1721:           case PETSC_OFFLOAD_CPU:
1722:             if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1723:             break;
1724:           case PETSC_OFFLOAD_UNALLOCATED:
1725:           case PETSC_OFFLOAD_KOKKOS:
1726:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1727:           }
1728: #endif
1729:         } else {
1730:           /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1731:           PetscCall(VecResetArray(*Y));
1732:         }
1733:         PetscCall(PetscObjectStateIncrease((PetscObject)X));
1734:       }
1735:     }
1736:   }
1737:   PetscCall(VecDestroy(Y));
1738:   PetscFunctionReturn(PETSC_SUCCESS);
1739: }

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

1745:   Not Collective.

1747:   Input Parameter:
1748: . v - The vector for which the local vector is desired.

1750:   Output Parameter:
1751: . w - Upon exit this contains the local vector.

1753:   Level: beginner

1755: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1756: @*/
1757: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1758: {
1759:   PetscMPIInt size;

1761:   PetscFunctionBegin;
1763:   PetscAssertPointer(w, 2);
1764:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1765:   if (size == 1) PetscCall(VecDuplicate(v, w));
1766:   else if (v->ops->createlocalvector) PetscUseTypeMethod(v, createlocalvector, w);
1767:   else {
1768:     VecType  type;
1769:     PetscInt n;

1771:     PetscCall(VecCreate(PETSC_COMM_SELF, w));
1772:     PetscCall(VecGetLocalSize(v, &n));
1773:     PetscCall(VecSetSizes(*w, n, n));
1774:     PetscCall(VecGetBlockSize(v, &n));
1775:     PetscCall(VecSetBlockSize(*w, n));
1776:     PetscCall(VecGetType(v, &type));
1777:     PetscCall(VecSetType(*w, type));
1778:   }
1779:   PetscFunctionReturn(PETSC_SUCCESS);
1780: }

1782: /*@
1783:   VecGetLocalVectorRead - Maps the local portion of a vector into a
1784:   vector.

1786:   Not Collective.

1788:   Input Parameter:
1789: . v - The vector for which the local vector is desired.

1791:   Output Parameter:
1792: . w - Upon exit this contains the local vector.

1794:   Level: beginner

1796:   Notes:
1797:   You must call `VecRestoreLocalVectorRead()` when the local
1798:   vector is no longer needed.

1800:   This function is similar to `VecGetArrayRead()` which maps the local
1801:   portion into a raw pointer.  `VecGetLocalVectorRead()` is usually
1802:   almost as efficient as `VecGetArrayRead()` but in certain circumstances
1803:   `VecGetLocalVectorRead()` can be much more efficient than
1804:   `VecGetArrayRead()`.  This is because the construction of a contiguous
1805:   array representing the vector data required by `VecGetArrayRead()` can
1806:   be an expensive operation for certain vector types.  For example, for
1807:   GPU vectors `VecGetArrayRead()` requires that the data between device
1808:   and host is synchronized.

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

1813: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1814: @*/
1815: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1816: {
1817:   PetscFunctionBegin;
1820:   VecCheckSameLocalSize(v, 1, w, 2);
1821:   if (v->ops->getlocalvectorread) {
1822:     PetscUseTypeMethod(v, getlocalvectorread, w);
1823:   } else {
1824:     PetscScalar *a;

1826:     PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1827:     PetscCall(VecPlaceArray(w, a));
1828:   }
1829:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
1830:   PetscCall(VecLockReadPush(v));
1831:   PetscCall(VecLockReadPush(w));
1832:   PetscFunctionReturn(PETSC_SUCCESS);
1833: }

1835: /*@
1836:   VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1837:   previously mapped into a vector using `VecGetLocalVectorRead()`.

1839:   Not Collective.

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

1845:   Level: beginner

1847: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1848: @*/
1849: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1850: {
1851:   PetscFunctionBegin;
1854:   if (v->ops->restorelocalvectorread) {
1855:     PetscUseTypeMethod(v, restorelocalvectorread, w);
1856:   } else {
1857:     const PetscScalar *a;

1859:     PetscCall(VecGetArrayRead(w, &a));
1860:     PetscCall(VecRestoreArrayRead(v, &a));
1861:     PetscCall(VecResetArray(w));
1862:   }
1863:   PetscCall(VecLockReadPop(v));
1864:   PetscCall(VecLockReadPop(w));
1865:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
1866:   PetscFunctionReturn(PETSC_SUCCESS);
1867: }

1869: /*@
1870:   VecGetLocalVector - Maps the local portion of a vector into a
1871:   vector.

1873:   Collective

1875:   Input Parameter:
1876: . v - The vector for which the local vector is desired.

1878:   Output Parameter:
1879: . w - Upon exit this contains the local vector.

1881:   Level: beginner

1883:   Notes:
1884:   You must call `VecRestoreLocalVector()` when the local
1885:   vector is no longer needed.

1887:   This function is similar to `VecGetArray()` which maps the local
1888:   portion into a raw pointer.  `VecGetLocalVector()` is usually about as
1889:   efficient as `VecGetArray()` but in certain circumstances
1890:   `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1891:   This is because the construction of a contiguous array representing
1892:   the vector data required by `VecGetArray()` can be an expensive
1893:   operation for certain vector types.  For example, for GPU vectors
1894:   `VecGetArray()` requires that the data between device and host is
1895:   synchronized.

1897: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1898: @*/
1899: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1900: {
1901:   PetscFunctionBegin;
1904:   VecCheckSameLocalSize(v, 1, w, 2);
1905:   if (v->ops->getlocalvector) {
1906:     PetscUseTypeMethod(v, getlocalvector, w);
1907:   } else {
1908:     PetscScalar *a;

1910:     PetscCall(VecGetArray(v, &a));
1911:     PetscCall(VecPlaceArray(w, a));
1912:   }
1913:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
1914:   PetscFunctionReturn(PETSC_SUCCESS);
1915: }

1917: /*@
1918:   VecRestoreLocalVector - Unmaps the local portion of a vector
1919:   previously mapped into a vector using `VecGetLocalVector()`.

1921:   Logically Collective.

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

1927:   Level: beginner

1929: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1930: @*/
1931: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1932: {
1933:   PetscFunctionBegin;
1936:   if (v->ops->restorelocalvector) {
1937:     PetscUseTypeMethod(v, restorelocalvector, w);
1938:   } else {
1939:     PetscScalar *a;
1940:     PetscCall(VecGetArray(w, &a));
1941:     PetscCall(VecRestoreArray(v, &a));
1942:     PetscCall(VecResetArray(w));
1943:   }
1944:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
1945:   PetscCall(PetscObjectStateIncrease((PetscObject)v));
1946:   PetscFunctionReturn(PETSC_SUCCESS);
1947: }

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

1953:   Logically Collective

1955:   Input Parameter:
1956: . x - the vector

1958:   Output Parameter:
1959: . a - location to put pointer to the array

1961:   Level: beginner

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

1969:   Fortran Notes:
1970:   `VecGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayF90()`

1972: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
1973:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
1974: @*/
1975: PetscErrorCode VecGetArray(Vec x, PetscScalar **a)
1976: {
1977:   PetscFunctionBegin;
1979:   PetscCall(VecSetErrorIfLocked(x, 1));
1980:   if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
1981:     PetscUseTypeMethod(x, getarray, a);
1982:   } else if (x->petscnative) { /* VECSTANDARD */
1983:     *a = *((PetscScalar **)x->data);
1984:   } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
1985:   PetscFunctionReturn(PETSC_SUCCESS);
1986: }

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

1991:   Logically Collective

1993:   Input Parameters:
1994: + x - the vector
1995: - a - location of pointer to array obtained from `VecGetArray()`

1997:   Level: beginner

1999:   Fortran Notes:
2000:   `VecRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayF90()`

2002: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2003:           `VecGetArrayPair()`, `VecRestoreArrayPair()`
2004: @*/
2005: PetscErrorCode VecRestoreArray(Vec x, PetscScalar **a)
2006: {
2007:   PetscFunctionBegin;
2009:   if (a) PetscAssertPointer(a, 2);
2010:   if (x->ops->restorearray) {
2011:     PetscUseTypeMethod(x, restorearray, a);
2012:   } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
2013:   if (a) *a = NULL;
2014:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
2015:   PetscFunctionReturn(PETSC_SUCCESS);
2016: }
2017: /*@C
2018:   VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.

2020:   Not Collective

2022:   Input Parameter:
2023: . x - the vector

2025:   Output Parameter:
2026: . a - the array

2028:   Level: beginner

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

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

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

2039:   Fortran Notes:
2040:   `VecGetArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayReadF90()`

2042: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2043: @*/
2044: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar **a)
2045: {
2046:   PetscFunctionBegin;
2048:   PetscAssertPointer(a, 2);
2049:   if (x->ops->getarrayread) {
2050:     PetscUseTypeMethod(x, getarrayread, a);
2051:   } else if (x->ops->getarray) {
2052:     PetscObjectState state;

2054:     /* VECNEST, VECCUDA, VECKOKKOS etc */
2055:     // x->ops->getarray may bump the object state, but since we know this is a read-only get
2056:     // we can just undo that
2057:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2058:     PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
2059:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
2060:   } else if (x->petscnative) {
2061:     /* VECSTANDARD */
2062:     *a = *((PetscScalar **)x->data);
2063:   } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
2064:   PetscFunctionReturn(PETSC_SUCCESS);
2065: }

2067: /*@C
2068:   VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`

2070:   Not Collective

2072:   Input Parameters:
2073: + x - the vector
2074: - a - the array

2076:   Level: beginner

2078:   Fortran Notes:
2079:   `VecRestoreArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayReadF90()`

2081: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2082: @*/
2083: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar **a)
2084: {
2085:   PetscFunctionBegin;
2087:   if (a) PetscAssertPointer(a, 2);
2088:   if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
2089:     /* nothing */
2090:   } else if (x->ops->restorearrayread) { /* VECNEST */
2091:     PetscUseTypeMethod(x, restorearrayread, a);
2092:   } else { /* No one? */
2093:     PetscObjectState state;

2095:     // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2096:     // we can just undo that
2097:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2098:     PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2099:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
2100:   }
2101:   if (a) *a = NULL;
2102:   PetscFunctionReturn(PETSC_SUCCESS);
2103: }

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

2109:   Logically Collective

2111:   Input Parameter:
2112: . x - the vector

2114:   Output Parameter:
2115: . a - location to put pointer to the array

2117:   Level: intermediate

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

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

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

2128:   Fortran Notes:
2129:   `VecGetArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayWriteF90()`

2131: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteF90()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
2132:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
2133: @*/
2134: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar **a)
2135: {
2136:   PetscFunctionBegin;
2138:   PetscAssertPointer(a, 2);
2139:   PetscCall(VecSetErrorIfLocked(x, 1));
2140:   if (x->ops->getarraywrite) {
2141:     PetscUseTypeMethod(x, getarraywrite, a);
2142:   } else {
2143:     PetscCall(VecGetArray(x, a));
2144:   }
2145:   PetscFunctionReturn(PETSC_SUCCESS);
2146: }

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

2151:   Logically Collective

2153:   Input Parameters:
2154: + x - the vector
2155: - a - location of pointer to array obtained from `VecGetArray()`

2157:   Level: beginner

2159:   Fortran Notes:
2160:   `VecRestoreArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayWriteF90()`

2162: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteF90()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2163:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2164: @*/
2165: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar **a)
2166: {
2167:   PetscFunctionBegin;
2169:   if (a) PetscAssertPointer(a, 2);
2170:   if (x->ops->restorearraywrite) {
2171:     PetscUseTypeMethod(x, restorearraywrite, a);
2172:   } else if (x->ops->restorearray) {
2173:     PetscUseTypeMethod(x, restorearray, a);
2174:   }
2175:   if (a) *a = NULL;
2176:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
2177:   PetscFunctionReturn(PETSC_SUCCESS);
2178: }

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

2184:   Logically Collective; No Fortran Support

2186:   Input Parameters:
2187: + x - the vectors
2188: - n - the number of vectors

2190:   Output Parameter:
2191: . a - location to put pointer to the array

2193:   Level: intermediate

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

2198: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2199: @*/
2200: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2201: {
2202:   PetscInt      i;
2203:   PetscScalar **q;

2205:   PetscFunctionBegin;
2206:   PetscAssertPointer(x, 1);
2208:   PetscAssertPointer(a, 3);
2209:   PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2210:   PetscCall(PetscMalloc1(n, &q));
2211:   for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2212:   *a = q;
2213:   PetscFunctionReturn(PETSC_SUCCESS);
2214: }

2216: /*@C
2217:   VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2218:   has been called.

2220:   Logically Collective; No Fortran Support

2222:   Input Parameters:
2223: + x - the vector
2224: . n - the number of vectors
2225: - a - location of pointer to arrays obtained from `VecGetArrays()`

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

2233:   Level: intermediate

2235: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2236: @*/
2237: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2238: {
2239:   PetscInt      i;
2240:   PetscScalar **q = *a;

2242:   PetscFunctionBegin;
2243:   PetscAssertPointer(x, 1);
2245:   PetscAssertPointer(a, 3);

2247:   for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2248:   PetscCall(PetscFree(q));
2249:   PetscFunctionReturn(PETSC_SUCCESS);
2250: }

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

2257:   Logically Collective; No Fortran Support

2259:   Input Parameter:
2260: . x - the vector

2262:   Output Parameters:
2263: + a     - location to put pointer to the array
2264: - mtype - memory type of the array

2266:   Level: beginner

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

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

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

2279: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`,
2280:           `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2281: @*/
2282: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2283: {
2284:   PetscFunctionBegin;
2287:   PetscAssertPointer(a, 2);
2288:   if (mtype) PetscAssertPointer(mtype, 3);
2289:   PetscCall(VecSetErrorIfLocked(x, 1));
2290:   if (x->ops->getarrayandmemtype) {
2291:     /* VECCUDA, VECKOKKOS etc */
2292:     PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2293:   } else {
2294:     /* VECSTANDARD, VECNEST, VECVIENNACL */
2295:     PetscCall(VecGetArray(x, a));
2296:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2297:   }
2298:   PetscFunctionReturn(PETSC_SUCCESS);
2299: }

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

2304:   Logically Collective; No Fortran Support

2306:   Input Parameters:
2307: + x - the vector
2308: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`

2310:   Level: beginner

2312: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`,
2313:           `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2314: @*/
2315: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar **a)
2316: {
2317:   PetscFunctionBegin;
2320:   if (a) PetscAssertPointer(a, 2);
2321:   if (x->ops->restorearrayandmemtype) {
2322:     /* VECCUDA, VECKOKKOS etc */
2323:     PetscUseTypeMethod(x, restorearrayandmemtype, a);
2324:   } else {
2325:     /* VECNEST, VECVIENNACL */
2326:     PetscCall(VecRestoreArray(x, a));
2327:   } /* VECSTANDARD does nothing */
2328:   if (a) *a = NULL;
2329:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
2330:   PetscFunctionReturn(PETSC_SUCCESS);
2331: }

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

2337:   Not Collective; No Fortran Support

2339:   Input Parameter:
2340: . x - the vector

2342:   Output Parameters:
2343: + a     - the array
2344: - mtype - memory type of the array

2346:   Level: beginner

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

2351: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2352: @*/
2353: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar **a, PetscMemType *mtype)
2354: {
2355:   PetscFunctionBegin;
2358:   PetscAssertPointer(a, 2);
2359:   if (mtype) PetscAssertPointer(mtype, 3);
2360:   if (x->ops->getarrayreadandmemtype) {
2361:     /* VECCUDA/VECHIP though they are also petscnative */
2362:     PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2363:   } else if (x->ops->getarrayandmemtype) {
2364:     /* VECKOKKOS */
2365:     PetscObjectState state;

2367:     // see VecGetArrayRead() for why
2368:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2369:     PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2370:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
2371:   } else {
2372:     PetscCall(VecGetArrayRead(x, a));
2373:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2374:   }
2375:   PetscFunctionReturn(PETSC_SUCCESS);
2376: }

2378: /*@C
2379:   VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`

2381:   Not Collective; No Fortran Support

2383:   Input Parameters:
2384: + x - the vector
2385: - a - the array

2387:   Level: beginner

2389: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2390: @*/
2391: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar **a)
2392: {
2393:   PetscFunctionBegin;
2396:   if (a) PetscAssertPointer(a, 2);
2397:   if (x->ops->restorearrayreadandmemtype) {
2398:     /* VECCUDA/VECHIP */
2399:     PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2400:   } else if (!x->petscnative) {
2401:     /* VECNEST */
2402:     PetscCall(VecRestoreArrayRead(x, a));
2403:   }
2404:   if (a) *a = NULL;
2405:   PetscFunctionReturn(PETSC_SUCCESS);
2406: }

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

2412:   Logically Collective; No Fortran Support

2414:   Input Parameter:
2415: . x - the vector

2417:   Output Parameters:
2418: + a     - the array
2419: - mtype - memory type of the array

2421:   Level: beginner

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

2426: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2427: @*/
2428: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2429: {
2430:   PetscFunctionBegin;
2433:   PetscCall(VecSetErrorIfLocked(x, 1));
2434:   PetscAssertPointer(a, 2);
2435:   if (mtype) PetscAssertPointer(mtype, 3);
2436:   if (x->ops->getarraywriteandmemtype) {
2437:     /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2438:     PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2439:   } else if (x->ops->getarrayandmemtype) {
2440:     PetscCall(VecGetArrayAndMemType(x, a, mtype));
2441:   } else {
2442:     /* VECNEST, VECVIENNACL */
2443:     PetscCall(VecGetArrayWrite(x, a));
2444:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2445:   }
2446:   PetscFunctionReturn(PETSC_SUCCESS);
2447: }

2449: /*@C
2450:   VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`

2452:   Logically Collective; No Fortran Support

2454:   Input Parameters:
2455: + x - the vector
2456: - a - the array

2458:   Level: beginner

2460: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2461: @*/
2462: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar **a)
2463: {
2464:   PetscFunctionBegin;
2467:   PetscCall(VecSetErrorIfLocked(x, 1));
2468:   if (a) PetscAssertPointer(a, 2);
2469:   if (x->ops->restorearraywriteandmemtype) {
2470:     /* VECCUDA/VECHIP */
2471:     PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2472:     PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2473:   } else if (x->ops->restorearrayandmemtype) {
2474:     PetscCall(VecRestoreArrayAndMemType(x, a));
2475:   } else {
2476:     PetscCall(VecRestoreArray(x, a));
2477:   }
2478:   if (a) *a = NULL;
2479:   PetscFunctionReturn(PETSC_SUCCESS);
2480: }

2482: /*@
2483:   VecPlaceArray - Allows one to replace the array in a vector with an
2484:   array provided by the user. This is useful to avoid copying an array
2485:   into a vector.

2487:   Logically Collective; No Fortran Support

2489:   Input Parameters:
2490: + vec   - the vector
2491: - array - the array

2493:   Level: developer

2495:   Notes:
2496:   Use `VecReplaceArray()` instead to permanently replace the array

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

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

2505: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2506: @*/
2507: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2508: {
2509:   PetscFunctionBegin;
2512:   if (array) PetscAssertPointer(array, 2);
2513:   PetscUseTypeMethod(vec, placearray, array);
2514:   PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2515:   PetscFunctionReturn(PETSC_SUCCESS);
2516: }

2518: /*@C
2519:   VecReplaceArray - Allows one to replace the array in a vector with an
2520:   array provided by the user. This is useful to avoid copying an array
2521:   into a vector.

2523:   Logically Collective; No Fortran Support

2525:   Input Parameters:
2526: + vec   - the vector
2527: - array - the array

2529:   Level: developer

2531:   Notes:
2532:   This permanently replaces the array and frees the memory associated
2533:   with the old array. Use `VecPlaceArray()` to temporarily replace the array.

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

2538: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2539: @*/
2540: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2541: {
2542:   PetscFunctionBegin;
2545:   PetscUseTypeMethod(vec, replacearray, array);
2546:   PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2547:   PetscFunctionReturn(PETSC_SUCCESS);
2548: }

2550: /*MC
2551:     VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
2552:     and makes them accessible via a Fortran pointer.

2554:     Synopsis:
2555:     VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)

2557:     Collective

2559:     Input Parameters:
2560: +   x - a vector to mimic
2561: -   n - the number of vectors to obtain

2563:     Output Parameters:
2564: +   y - Fortran pointer to the array of vectors
2565: -   ierr - error code

2567:     Example of Usage:
2568: .vb
2569: #include <petsc/finclude/petscvec.h>
2570:     use petscvec

2572:     Vec x
2573:     Vec, pointer :: y(:)
2574:     ....
2575:     call VecDuplicateVecsF90(x,2,y,ierr)
2576:     call VecSet(y(2),alpha,ierr)
2577:     call VecSet(y(2),alpha,ierr)
2578:     ....
2579:     call VecDestroyVecsF90(2,y,ierr)
2580: .ve

2582:     Level: beginner

2584:     Note:
2585:     Use `VecDestroyVecsF90()` to free the space.

2587: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecsF90()`, `VecDuplicateVecs()`
2588: M*/

2590: /*MC
2591:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2592:     `VecGetArrayF90()`.

2594:     Synopsis:
2595:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2597:     Logically Collective

2599:     Input Parameters:
2600: +   x - vector
2601: -   xx_v - the Fortran pointer to the array

2603:     Output Parameter:
2604: .   ierr - error code

2606:     Example of Usage:
2607: .vb
2608: #include <petsc/finclude/petscvec.h>
2609:     use petscvec

2611:     PetscScalar, pointer :: xx_v(:)
2612:     ....
2613:     call VecGetArrayF90(x,xx_v,ierr)
2614:     xx_v(3) = a
2615:     call VecRestoreArrayF90(x,xx_v,ierr)
2616: .ve

2618:     Level: beginner

2620: .seealso: [](ch_vectors), `Vec`, `VecGetArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrayReadF90()`
2621: M*/

2623: /*MC
2624:     VecDestroyVecsF90 - Frees a block of vectors obtained with `VecDuplicateVecsF90()`.

2626:     Synopsis:
2627:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

2629:     Collective

2631:     Input Parameters:
2632: +   n - the number of vectors previously obtained
2633: -   x - pointer to array of vector pointers

2635:     Output Parameter:
2636: .   ierr - error code

2638:     Level: beginner

2640: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecs()`, `VecDuplicateVecsF90()`
2641: M*/

2643: /*MC
2644:     VecGetArrayF90 - Accesses a vector array from Fortran. For default PETSc
2645:     vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2646:     this routine is implementation dependent. You MUST call `VecRestoreArrayF90()`
2647:     when you no longer need access to the array.

2649:     Synopsis:
2650:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2652:     Logically Collective

2654:     Input Parameter:
2655: .   x - vector

2657:     Output Parameters:
2658: +   xx_v - the Fortran pointer to the array
2659: -   ierr - error code

2661:     Example of Usage:
2662: .vb
2663: #include <petsc/finclude/petscvec.h>
2664:     use petscvec

2666:     PetscScalar, pointer :: xx_v(:)
2667:     ....
2668:     call VecGetArrayF90(x,xx_v,ierr)
2669:     xx_v(3) = a
2670:     call VecRestoreArrayF90(x,xx_v,ierr)
2671: .ve

2673:      Level: beginner

2675:     Note:
2676:     If you ONLY intend to read entries from the array and not change any entries you should use `VecGetArrayReadF90()`.

2678: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayReadF90()`
2679: M*/

2681: /*MC
2682:     VecGetArrayReadF90 - Accesses a read only array from Fortran. For default PETSc
2683:     vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2684:     this routine is implementation dependent. You MUST call `VecRestoreArrayReadF90()`
2685:     when you no longer need access to the array.

2687:     Synopsis:
2688:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2690:     Logically Collective

2692:     Input Parameter:
2693: .   x - vector

2695:     Output Parameters:
2696: +   xx_v - the Fortran pointer to the array
2697: -   ierr - error code

2699:     Example of Usage:
2700: .vb
2701: #include <petsc/finclude/petscvec.h>
2702:     use petscvec

2704:     PetscScalar, pointer :: xx_v(:)
2705:     ....
2706:     call VecGetArrayReadF90(x,xx_v,ierr)
2707:     a = xx_v(3)
2708:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2709: .ve

2711:     Level: beginner

2713:     Note:
2714:     If you intend to write entries into the array you must use `VecGetArrayF90()`.

2716: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecGetArrayF90()`
2717: M*/

2719: /*MC
2720:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2721:     `VecGetArrayReadF90()`.

2723:     Synopsis:
2724:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2726:     Logically Collective

2728:     Input Parameters:
2729: +   x - vector
2730: -   xx_v - the Fortran pointer to the array

2732:     Output Parameter:
2733: .   ierr - error code

2735:     Example of Usage:
2736: .vb
2737: #include <petsc/finclude/petscvec.h>
2738:     use petscvec

2740:     PetscScalar, pointer :: xx_v(:)
2741:     ....
2742:     call VecGetArrayReadF90(x,xx_v,ierr)
2743:     a = xx_v(3)
2744:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2745: .ve

2747:     Level: beginner

2749: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecRestoreArrayF90()`
2750: M*/

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

2757:   Logically Collective

2759:   Input Parameters:
2760: + x      - the vector
2761: . m      - first dimension of two dimensional array
2762: . n      - second dimension of two dimensional array
2763: . mstart - first index you will use in first coordinate direction (often 0)
2764: - nstart - first index in the second coordinate direction (often 0)

2766:   Output Parameter:
2767: . a - location to put pointer to the array

2769:   Level: developer

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

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

2779: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2780:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2781:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2782: @*/
2783: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2784: {
2785:   PetscInt     i, N;
2786:   PetscScalar *aa;

2788:   PetscFunctionBegin;
2790:   PetscAssertPointer(a, 6);
2792:   PetscCall(VecGetLocalSize(x, &N));
2793:   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);
2794:   PetscCall(VecGetArray(x, &aa));

2796:   PetscCall(PetscMalloc1(m, a));
2797:   for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2798:   *a -= mstart;
2799:   PetscFunctionReturn(PETSC_SUCCESS);
2800: }

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

2807:   Logically Collective

2809:   Input Parameters:
2810: + x      - the vector
2811: . m      - first dimension of two dimensional array
2812: . n      - second dimension of two dimensional array
2813: . mstart - first index you will use in first coordinate direction (often 0)
2814: - nstart - first index in the second coordinate direction (often 0)

2816:   Output Parameter:
2817: . a - location to put pointer to the array

2819:   Level: developer

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

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

2829: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2830:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2831:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2832: @*/
2833: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2834: {
2835:   PetscInt     i, N;
2836:   PetscScalar *aa;

2838:   PetscFunctionBegin;
2840:   PetscAssertPointer(a, 6);
2842:   PetscCall(VecGetLocalSize(x, &N));
2843:   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);
2844:   PetscCall(VecGetArrayWrite(x, &aa));

2846:   PetscCall(PetscMalloc1(m, a));
2847:   for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2848:   *a -= mstart;
2849:   PetscFunctionReturn(PETSC_SUCCESS);
2850: }

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

2855:   Logically Collective

2857:   Input Parameters:
2858: + x      - the vector
2859: . m      - first dimension of two dimensional array
2860: . n      - second dimension of the two dimensional array
2861: . mstart - first index you will use in first coordinate direction (often 0)
2862: . nstart - first index in the second coordinate direction (often 0)
2863: - a      - location of pointer to array obtained from `VecGetArray2d()`

2865:   Level: developer

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

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

2875: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2876:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2877:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2878: @*/
2879: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2880: {
2881:   void *dummy;

2883:   PetscFunctionBegin;
2885:   PetscAssertPointer(a, 6);
2887:   dummy = (void *)(*a + mstart);
2888:   PetscCall(PetscFree(dummy));
2889:   PetscCall(VecRestoreArray(x, NULL));
2890:   *a = NULL;
2891:   PetscFunctionReturn(PETSC_SUCCESS);
2892: }

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

2897:   Logically Collective

2899:   Input Parameters:
2900: + x      - the vector
2901: . m      - first dimension of two dimensional array
2902: . n      - second dimension of the two dimensional array
2903: . mstart - first index you will use in first coordinate direction (often 0)
2904: . nstart - first index in the second coordinate direction (often 0)
2905: - a      - location of pointer to array obtained from `VecGetArray2d()`

2907:   Level: developer

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

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

2917: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2918:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2919:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2920: @*/
2921: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2922: {
2923:   void *dummy;

2925:   PetscFunctionBegin;
2927:   PetscAssertPointer(a, 6);
2929:   dummy = (void *)(*a + mstart);
2930:   PetscCall(PetscFree(dummy));
2931:   PetscCall(VecRestoreArrayWrite(x, NULL));
2932:   PetscFunctionReturn(PETSC_SUCCESS);
2933: }

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

2940:   Logically Collective

2942:   Input Parameters:
2943: + x      - the vector
2944: . m      - first dimension of two dimensional array
2945: - mstart - first index you will use in first coordinate direction (often 0)

2947:   Output Parameter:
2948: . a - location to put pointer to the array

2950:   Level: developer

2952:   Notes:
2953:   For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2954:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2955:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

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

2959: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2960:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2961:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2962: @*/
2963: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2964: {
2965:   PetscInt N;

2967:   PetscFunctionBegin;
2969:   PetscAssertPointer(a, 4);
2971:   PetscCall(VecGetLocalSize(x, &N));
2972:   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);
2973:   PetscCall(VecGetArray(x, a));
2974:   *a -= mstart;
2975:   PetscFunctionReturn(PETSC_SUCCESS);
2976: }

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

2983:   Logically Collective

2985:   Input Parameters:
2986: + x      - the vector
2987: . m      - first dimension of two dimensional array
2988: - mstart - first index you will use in first coordinate direction (often 0)

2990:   Output Parameter:
2991: . a - location to put pointer to the array

2993:   Level: developer

2995:   Notes:
2996:   For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2997:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2998:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

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

3002: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3003:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3004:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3005: @*/
3006: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3007: {
3008:   PetscInt N;

3010:   PetscFunctionBegin;
3012:   PetscAssertPointer(a, 4);
3014:   PetscCall(VecGetLocalSize(x, &N));
3015:   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);
3016:   PetscCall(VecGetArrayWrite(x, a));
3017:   *a -= mstart;
3018:   PetscFunctionReturn(PETSC_SUCCESS);
3019: }

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

3024:   Logically Collective

3026:   Input Parameters:
3027: + x      - the vector
3028: . m      - first dimension of two dimensional array
3029: . mstart - first index you will use in first coordinate direction (often 0)
3030: - a      - location of pointer to array obtained from `VecGetArray1d()`

3032:   Level: developer

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

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

3042: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3043:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3044:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3045: @*/
3046: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3047: {
3048:   PetscFunctionBegin;
3051:   PetscCall(VecRestoreArray(x, NULL));
3052:   *a = NULL;
3053:   PetscFunctionReturn(PETSC_SUCCESS);
3054: }

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

3059:   Logically Collective

3061:   Input Parameters:
3062: + x      - the vector
3063: . m      - first dimension of two dimensional array
3064: . mstart - first index you will use in first coordinate direction (often 0)
3065: - a      - location of pointer to array obtained from `VecGetArray1d()`

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 `VecGetArray1d()`.

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

3077: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3078:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3079:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3080: @*/
3081: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3082: {
3083:   PetscFunctionBegin;
3086:   PetscCall(VecRestoreArrayWrite(x, NULL));
3087:   *a = NULL;
3088:   PetscFunctionReturn(PETSC_SUCCESS);
3089: }

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

3096:   Logically Collective

3098:   Input Parameters:
3099: + x      - the vector
3100: . m      - first dimension of three dimensional array
3101: . n      - second dimension of three dimensional array
3102: . p      - third dimension of three dimensional array
3103: . mstart - first index you will use in first coordinate direction (often 0)
3104: . nstart - first index in the second coordinate direction (often 0)
3105: - pstart - first index in the third coordinate direction (often 0)

3107:   Output Parameter:
3108: . a - location to put pointer to the array

3110:   Level: developer

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

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

3120: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3121:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
3122:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3123: @*/
3124: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3125: {
3126:   PetscInt     i, N, j;
3127:   PetscScalar *aa, **b;

3129:   PetscFunctionBegin;
3131:   PetscAssertPointer(a, 8);
3133:   PetscCall(VecGetLocalSize(x, &N));
3134:   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);
3135:   PetscCall(VecGetArray(x, &aa));

3137:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3138:   b = (PetscScalar **)((*a) + m);
3139:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3140:   for (i = 0; i < m; i++)
3141:     for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3142:   *a -= mstart;
3143:   PetscFunctionReturn(PETSC_SUCCESS);
3144: }

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

3151:   Logically Collective

3153:   Input Parameters:
3154: + x      - the vector
3155: . m      - first dimension of three dimensional array
3156: . n      - second dimension of three dimensional array
3157: . p      - third dimension of three dimensional array
3158: . mstart - first index you will use in first coordinate direction (often 0)
3159: . nstart - first index in the second coordinate direction (often 0)
3160: - pstart - first index in the third coordinate direction (often 0)

3162:   Output Parameter:
3163: . a - location to put pointer to the array

3165:   Level: developer

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

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

3175: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3176:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3177:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3178: @*/
3179: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3180: {
3181:   PetscInt     i, N, j;
3182:   PetscScalar *aa, **b;

3184:   PetscFunctionBegin;
3186:   PetscAssertPointer(a, 8);
3188:   PetscCall(VecGetLocalSize(x, &N));
3189:   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);
3190:   PetscCall(VecGetArrayWrite(x, &aa));

3192:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3193:   b = (PetscScalar **)((*a) + m);
3194:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3195:   for (i = 0; i < m; i++)
3196:     for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;

3198:   *a -= mstart;
3199:   PetscFunctionReturn(PETSC_SUCCESS);
3200: }

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

3205:   Logically Collective

3207:   Input Parameters:
3208: + x      - the vector
3209: . m      - first dimension of three dimensional array
3210: . n      - second dimension of the three dimensional array
3211: . p      - third dimension of the three dimensional array
3212: . mstart - first index you will use in first coordinate direction (often 0)
3213: . nstart - first index in the second coordinate direction (often 0)
3214: . pstart - first index in the third coordinate direction (often 0)
3215: - a      - location of pointer to array obtained from VecGetArray3d()

3217:   Level: developer

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

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

3227: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3228:           `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3229:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3230: @*/
3231: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3232: {
3233:   void *dummy;

3235:   PetscFunctionBegin;
3237:   PetscAssertPointer(a, 8);
3239:   dummy = (void *)(*a + mstart);
3240:   PetscCall(PetscFree(dummy));
3241:   PetscCall(VecRestoreArray(x, NULL));
3242:   *a = NULL;
3243:   PetscFunctionReturn(PETSC_SUCCESS);
3244: }

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

3249:   Logically Collective

3251:   Input Parameters:
3252: + x      - the vector
3253: . m      - first dimension of three dimensional array
3254: . n      - second dimension of the three dimensional array
3255: . p      - third dimension of the three dimensional array
3256: . mstart - first index you will use in first coordinate direction (often 0)
3257: . nstart - first index in the second coordinate direction (often 0)
3258: . pstart - first index in the third coordinate direction (often 0)
3259: - a      - location of pointer to array obtained from VecGetArray3d()

3261:   Level: developer

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

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

3271: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3272:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3273:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3274: @*/
3275: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3276: {
3277:   void *dummy;

3279:   PetscFunctionBegin;
3281:   PetscAssertPointer(a, 8);
3283:   dummy = (void *)(*a + mstart);
3284:   PetscCall(PetscFree(dummy));
3285:   PetscCall(VecRestoreArrayWrite(x, NULL));
3286:   *a = NULL;
3287:   PetscFunctionReturn(PETSC_SUCCESS);
3288: }

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

3295:   Logically Collective

3297:   Input Parameters:
3298: + x      - the vector
3299: . m      - first dimension of four dimensional array
3300: . n      - second dimension of four dimensional array
3301: . p      - third dimension of four dimensional array
3302: . q      - fourth dimension of four dimensional array
3303: . mstart - first index you will use in first coordinate direction (often 0)
3304: . nstart - first index in the second coordinate direction (often 0)
3305: . pstart - first index in the third coordinate direction (often 0)
3306: - qstart - first index in the fourth coordinate direction (often 0)

3308:   Output Parameter:
3309: . a - location to put pointer to the array

3311:   Level: developer

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

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

3321: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3322:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3323:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3324: @*/
3325: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3326: {
3327:   PetscInt     i, N, j, k;
3328:   PetscScalar *aa, ***b, **c;

3330:   PetscFunctionBegin;
3332:   PetscAssertPointer(a, 10);
3334:   PetscCall(VecGetLocalSize(x, &N));
3335:   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);
3336:   PetscCall(VecGetArray(x, &aa));

3338:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3339:   b = (PetscScalar ***)((*a) + m);
3340:   c = (PetscScalar **)(b + m * n);
3341:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3342:   for (i = 0; i < m; i++)
3343:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3344:   for (i = 0; i < m; i++)
3345:     for (j = 0; j < n; j++)
3346:       for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3347:   *a -= mstart;
3348:   PetscFunctionReturn(PETSC_SUCCESS);
3349: }

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

3356:   Logically Collective

3358:   Input Parameters:
3359: + x      - the vector
3360: . m      - first dimension of four dimensional array
3361: . n      - second dimension of four dimensional array
3362: . p      - third dimension of four dimensional array
3363: . q      - fourth dimension of four dimensional array
3364: . mstart - first index you will use in first coordinate direction (often 0)
3365: . nstart - first index in the second coordinate direction (often 0)
3366: . pstart - first index in the third coordinate direction (often 0)
3367: - qstart - first index in the fourth coordinate direction (often 0)

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

3372:   Level: developer

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

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

3382: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3383:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3384:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3385: @*/
3386: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3387: {
3388:   PetscInt     i, N, j, k;
3389:   PetscScalar *aa, ***b, **c;

3391:   PetscFunctionBegin;
3393:   PetscAssertPointer(a, 10);
3395:   PetscCall(VecGetLocalSize(x, &N));
3396:   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);
3397:   PetscCall(VecGetArrayWrite(x, &aa));

3399:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3400:   b = (PetscScalar ***)((*a) + m);
3401:   c = (PetscScalar **)(b + m * n);
3402:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3403:   for (i = 0; i < m; i++)
3404:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3405:   for (i = 0; i < m; i++)
3406:     for (j = 0; j < n; j++)
3407:       for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3408:   *a -= mstart;
3409:   PetscFunctionReturn(PETSC_SUCCESS);
3410: }

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

3415:   Logically Collective

3417:   Input Parameters:
3418: + x      - the vector
3419: . m      - first dimension of four dimensional array
3420: . n      - second dimension of the four dimensional array
3421: . p      - third dimension of the four dimensional array
3422: . q      - fourth dimension of the four dimensional array
3423: . mstart - first index you will use in first coordinate direction (often 0)
3424: . nstart - first index in the second coordinate direction (often 0)
3425: . pstart - first index in the third coordinate direction (often 0)
3426: . qstart - first index in the fourth coordinate direction (often 0)
3427: - a      - location of pointer to array obtained from VecGetArray4d()

3429:   Level: developer

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

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

3439: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3440:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3441:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecGet`
3442: @*/
3443: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3444: {
3445:   void *dummy;

3447:   PetscFunctionBegin;
3449:   PetscAssertPointer(a, 10);
3451:   dummy = (void *)(*a + mstart);
3452:   PetscCall(PetscFree(dummy));
3453:   PetscCall(VecRestoreArray(x, NULL));
3454:   *a = NULL;
3455:   PetscFunctionReturn(PETSC_SUCCESS);
3456: }

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

3461:   Logically Collective

3463:   Input Parameters:
3464: + x      - the vector
3465: . m      - first dimension of four dimensional array
3466: . n      - second dimension of the four dimensional array
3467: . p      - third dimension of the four dimensional array
3468: . q      - fourth dimension of the four dimensional array
3469: . mstart - first index you will use in first coordinate direction (often 0)
3470: . nstart - first index in the second coordinate direction (often 0)
3471: . pstart - first index in the third coordinate direction (often 0)
3472: . qstart - first index in the fourth coordinate direction (often 0)
3473: - a      - location of pointer to array obtained from `VecGetArray4d()`

3475:   Level: developer

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

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

3485: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3486:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3487:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3488: @*/
3489: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3490: {
3491:   void *dummy;

3493:   PetscFunctionBegin;
3495:   PetscAssertPointer(a, 10);
3497:   dummy = (void *)(*a + mstart);
3498:   PetscCall(PetscFree(dummy));
3499:   PetscCall(VecRestoreArrayWrite(x, NULL));
3500:   *a = NULL;
3501:   PetscFunctionReturn(PETSC_SUCCESS);
3502: }

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

3509:   Logically Collective

3511:   Input Parameters:
3512: + x      - the vector
3513: . m      - first dimension of two dimensional array
3514: . n      - second dimension of two dimensional array
3515: . mstart - first index you will use in first coordinate direction (often 0)
3516: - nstart - first index in the second coordinate direction (often 0)

3518:   Output Parameter:
3519: . a - location to put pointer to the array

3521:   Level: developer

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

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

3531: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3532:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3533:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3534: @*/
3535: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3536: {
3537:   PetscInt           i, N;
3538:   const PetscScalar *aa;

3540:   PetscFunctionBegin;
3542:   PetscAssertPointer(a, 6);
3544:   PetscCall(VecGetLocalSize(x, &N));
3545:   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);
3546:   PetscCall(VecGetArrayRead(x, &aa));

3548:   PetscCall(PetscMalloc1(m, a));
3549:   for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3550:   *a -= mstart;
3551:   PetscFunctionReturn(PETSC_SUCCESS);
3552: }

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

3557:   Logically Collective

3559:   Input Parameters:
3560: + x      - the vector
3561: . m      - first dimension of two dimensional array
3562: . n      - second dimension of the two dimensional array
3563: . mstart - first index you will use in first coordinate direction (often 0)
3564: . nstart - first index in the second coordinate direction (often 0)
3565: - a      - location of pointer to array obtained from VecGetArray2d()

3567:   Level: developer

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

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

3577: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3578:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3579:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3580: @*/
3581: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3582: {
3583:   void *dummy;

3585:   PetscFunctionBegin;
3587:   PetscAssertPointer(a, 6);
3589:   dummy = (void *)(*a + mstart);
3590:   PetscCall(PetscFree(dummy));
3591:   PetscCall(VecRestoreArrayRead(x, NULL));
3592:   *a = NULL;
3593:   PetscFunctionReturn(PETSC_SUCCESS);
3594: }

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

3601:   Logically Collective

3603:   Input Parameters:
3604: + x      - the vector
3605: . m      - first dimension of two dimensional array
3606: - mstart - first index you will use in first coordinate direction (often 0)

3608:   Output Parameter:
3609: . a - location to put pointer to the array

3611:   Level: developer

3613:   Notes:
3614:   For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3615:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3616:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

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

3620: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3621:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3622:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3623: @*/
3624: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3625: {
3626:   PetscInt N;

3628:   PetscFunctionBegin;
3630:   PetscAssertPointer(a, 4);
3632:   PetscCall(VecGetLocalSize(x, &N));
3633:   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);
3634:   PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3635:   *a -= mstart;
3636:   PetscFunctionReturn(PETSC_SUCCESS);
3637: }

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

3642:   Logically Collective

3644:   Input Parameters:
3645: + x      - the vector
3646: . m      - first dimension of two dimensional array
3647: . mstart - first index you will use in first coordinate direction (often 0)
3648: - a      - location of pointer to array obtained from `VecGetArray1dRead()`

3650:   Level: developer

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

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

3660: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3661:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3662:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3663: @*/
3664: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3665: {
3666:   PetscFunctionBegin;
3669:   PetscCall(VecRestoreArrayRead(x, NULL));
3670:   *a = NULL;
3671:   PetscFunctionReturn(PETSC_SUCCESS);
3672: }

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

3679:   Logically Collective

3681:   Input Parameters:
3682: + x      - the vector
3683: . m      - first dimension of three dimensional array
3684: . n      - second dimension of three dimensional array
3685: . p      - third dimension of three dimensional array
3686: . mstart - first index you will use in first coordinate direction (often 0)
3687: . nstart - first index in the second coordinate direction (often 0)
3688: - pstart - first index in the third coordinate direction (often 0)

3690:   Output Parameter:
3691: . a - location to put pointer to the array

3693:   Level: developer

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

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

3703: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3704:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3705:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3706: @*/
3707: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3708: {
3709:   PetscInt           i, N, j;
3710:   const PetscScalar *aa;
3711:   PetscScalar      **b;

3713:   PetscFunctionBegin;
3715:   PetscAssertPointer(a, 8);
3717:   PetscCall(VecGetLocalSize(x, &N));
3718:   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);
3719:   PetscCall(VecGetArrayRead(x, &aa));

3721:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3722:   b = (PetscScalar **)((*a) + m);
3723:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3724:   for (i = 0; i < m; i++)
3725:     for (j = 0; j < n; j++) b[i * n + j] = (PetscScalar *)aa + i * n * p + j * p - pstart;
3726:   *a -= mstart;
3727:   PetscFunctionReturn(PETSC_SUCCESS);
3728: }

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

3733:   Logically Collective

3735:   Input Parameters:
3736: + x      - the vector
3737: . m      - first dimension of three dimensional array
3738: . n      - second dimension of the three dimensional array
3739: . p      - third dimension of the three dimensional array
3740: . mstart - first index you will use in first coordinate direction (often 0)
3741: . nstart - first index in the second coordinate direction (often 0)
3742: . pstart - first index in the third coordinate direction (often 0)
3743: - a      - location of pointer to array obtained from `VecGetArray3dRead()`

3745:   Level: developer

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

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

3755: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3756:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3757:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3758: @*/
3759: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3760: {
3761:   void *dummy;

3763:   PetscFunctionBegin;
3765:   PetscAssertPointer(a, 8);
3767:   dummy = (void *)(*a + mstart);
3768:   PetscCall(PetscFree(dummy));
3769:   PetscCall(VecRestoreArrayRead(x, NULL));
3770:   *a = NULL;
3771:   PetscFunctionReturn(PETSC_SUCCESS);
3772: }

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

3779:   Logically Collective

3781:   Input Parameters:
3782: + x      - the vector
3783: . m      - first dimension of four dimensional array
3784: . n      - second dimension of four dimensional array
3785: . p      - third dimension of four dimensional array
3786: . q      - fourth dimension of four dimensional array
3787: . mstart - first index you will use in first coordinate direction (often 0)
3788: . nstart - first index in the second coordinate direction (often 0)
3789: . pstart - first index in the third coordinate direction (often 0)
3790: - qstart - first index in the fourth coordinate direction (often 0)

3792:   Output Parameter:
3793: . a - location to put pointer to the array

3795:   Level: beginner

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

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

3805: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3806:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3807:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3808: @*/
3809: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3810: {
3811:   PetscInt           i, N, j, k;
3812:   const PetscScalar *aa;
3813:   PetscScalar     ***b, **c;

3815:   PetscFunctionBegin;
3817:   PetscAssertPointer(a, 10);
3819:   PetscCall(VecGetLocalSize(x, &N));
3820:   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);
3821:   PetscCall(VecGetArrayRead(x, &aa));

3823:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3824:   b = (PetscScalar ***)((*a) + m);
3825:   c = (PetscScalar **)(b + m * n);
3826:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3827:   for (i = 0; i < m; i++)
3828:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3829:   for (i = 0; i < m; i++)
3830:     for (j = 0; j < n; j++)
3831:       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;
3832:   *a -= mstart;
3833:   PetscFunctionReturn(PETSC_SUCCESS);
3834: }

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

3839:   Logically Collective

3841:   Input Parameters:
3842: + x      - the vector
3843: . m      - first dimension of four dimensional array
3844: . n      - second dimension of the four dimensional array
3845: . p      - third dimension of the four dimensional array
3846: . q      - fourth dimension of the four dimensional array
3847: . mstart - first index you will use in first coordinate direction (often 0)
3848: . nstart - first index in the second coordinate direction (often 0)
3849: . pstart - first index in the third coordinate direction (often 0)
3850: . qstart - first index in the fourth coordinate direction (often 0)
3851: - a      - location of pointer to array obtained from `VecGetArray4dRead()`

3853:   Level: beginner

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

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

3863: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3864:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3865:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3866: @*/
3867: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3868: {
3869:   void *dummy;

3871:   PetscFunctionBegin;
3873:   PetscAssertPointer(a, 10);
3875:   dummy = (void *)(*a + mstart);
3876:   PetscCall(PetscFree(dummy));
3877:   PetscCall(VecRestoreArrayRead(x, NULL));
3878:   *a = NULL;
3879:   PetscFunctionReturn(PETSC_SUCCESS);
3880: }

3882: #if defined(PETSC_USE_DEBUG)

3884: /*@
3885:   VecLockGet  - Gets the current lock status of a vector

3887:   Logically Collective

3889:   Input Parameter:
3890: . x - the vector

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

3896:   Level: advanced

3898: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3899: @*/
3900: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3901: {
3902:   PetscFunctionBegin;
3904:   *state = x->lock;
3905:   PetscFunctionReturn(PETSC_SUCCESS);
3906: }

3908: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3909: {
3910:   PetscFunctionBegin;
3912:   PetscAssertPointer(file, 2);
3913:   PetscAssertPointer(func, 3);
3914:   PetscAssertPointer(line, 4);
3915:   #if !PetscDefined(HAVE_THREADSAFETY)
3916:   {
3917:     const int index = x->lockstack.currentsize - 1;

3919:     PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupted vec lock stack, have negative index %d", index);
3920:     *file = x->lockstack.file[index];
3921:     *func = x->lockstack.function[index];
3922:     *line = x->lockstack.line[index];
3923:   }
3924:   #else
3925:   *file = NULL;
3926:   *func = NULL;
3927:   *line = 0;
3928:   #endif
3929:   PetscFunctionReturn(PETSC_SUCCESS);
3930: }

3932: /*@
3933:   VecLockReadPush  - Pushes a read-only lock on a vector to prevent it from being written to

3935:   Logically Collective

3937:   Input Parameter:
3938: . x - the vector

3940:   Level: intermediate

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

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

3948: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3949: @*/
3950: PetscErrorCode VecLockReadPush(Vec x)
3951: {
3952:   PetscFunctionBegin;
3954:   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");
3955:   #if !PetscDefined(HAVE_THREADSAFETY)
3956:   {
3957:     const char *file, *func;
3958:     int         index, line;

3960:     if ((index = petscstack.currentsize - 2) == -1) {
3961:       // vec was locked "outside" of petsc, either in user-land or main. the error message will
3962:       // now show this function as the culprit, but it will include the stacktrace
3963:       file = "unknown user-file";
3964:       func = "unknown_user_function";
3965:       line = 0;
3966:     } else {
3967:       PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected petscstack, have negative index %d", index);
3968:       file = petscstack.file[index];
3969:       func = petscstack.function[index];
3970:       line = petscstack.line[index];
3971:     }
3972:     PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
3973:   }
3974:   #endif
3975:   PetscFunctionReturn(PETSC_SUCCESS);
3976: }

3978: /*@
3979:   VecLockReadPop  - Pops a read-only lock from a vector

3981:   Logically Collective

3983:   Input Parameter:
3984: . x - the vector

3986:   Level: intermediate

3988: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
3989: @*/
3990: PetscErrorCode VecLockReadPop(Vec x)
3991: {
3992:   PetscFunctionBegin;
3994:   PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
3995:   #if !PetscDefined(HAVE_THREADSAFETY)
3996:   {
3997:     const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];

3999:     PetscStackPop_Private(x->lockstack, previous);
4000:   }
4001:   #endif
4002:   PetscFunctionReturn(PETSC_SUCCESS);
4003: }

4005: /*@C
4006:   VecLockWriteSet  - Lock or unlock a vector for exclusive read/write access

4008:   Logically Collective

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

4014:   Level: intermediate

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

4022: .vb
4023:        VecGetArray(x,&xdata); // begin phase
4024:        VecLockWriteSet(v,PETSC_TRUE);

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

4028:        VecRestoreArray(x,&vdata); // end phase
4029:        VecLockWriteSet(v,PETSC_FALSE);
4030: .ve

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

4035: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
4036: @*/
4037: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
4038: {
4039:   PetscFunctionBegin;
4041:   if (flg) {
4042:     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");
4043:     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");
4044:     x->lock = -1;
4045:   } else {
4046:     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");
4047:     x->lock = 0;
4048:   }
4049:   PetscFunctionReturn(PETSC_SUCCESS);
4050: }

4052: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
4053: /*@
4054:   VecLockPush  - Pushes a read-only lock on a vector to prevent it from being written to

4056:   Level: deprecated

4058: .seealso: [](ch_vectors), `Vec`, `VecLockReadPush()`
4059: @*/
4060: PetscErrorCode VecLockPush(Vec x)
4061: {
4062:   PetscFunctionBegin;
4063:   PetscCall(VecLockReadPush(x));
4064:   PetscFunctionReturn(PETSC_SUCCESS);
4065: }

4067: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
4068: /*@
4069:   VecLockPop  - Pops a read-only lock from a vector

4071:   Level: deprecated

4073: .seealso: [](ch_vectors), `Vec`, `VecLockReadPop()`
4074: @*/
4075: PetscErrorCode VecLockPop(Vec x)
4076: {
4077:   PetscFunctionBegin;
4078:   PetscCall(VecLockReadPop(x));
4079:   PetscFunctionReturn(PETSC_SUCCESS);
4080: }

4082: #endif