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:     PetscMPIInt b0 = (PetscMPIInt)flg, b1[2], b2[2];
212:     b1[0]          = -b0;
213:     b1[1]          = b0;
214:     PetscCall(MPIU_Allreduce(b1, b2, 2, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
215:     PetscCheck(-b2[0] == b2[1], 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.");
216:   }
217:   if (flg) PetscFunctionReturn(PETSC_SUCCESS);

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

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

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

232:   Not Collective

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

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

244:   Level: intermediate

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

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

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

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

273:   Collective

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

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

281:   Level: intermediate

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

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

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

310:   Collective

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

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

319:   Level: intermediate

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

324:   Returns the smallest index with the maximum value

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

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

347:   Collective

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

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

356:   Level: intermediate

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

361:   This returns the smallest index with the minimum value

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

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

385:   Collective

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

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

394:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

462:   Not Collective

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

468:   Level: intermediate

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

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

482: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
483: {
484:   PetscFunctionBegin;
487:   VecCheckAssembled(x);
489:   PetscCall(VecSetErrorIfLocked(x, 1));

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

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

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

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

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

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

530:   Logically Collective

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

536:   Level: beginner

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

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

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

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

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

587:   Logically Collective

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

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

597:   Level: intermediate

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

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

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

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

652:   Logically Collective

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

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

662:   Level: intermediate

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

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

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

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

707:   Logically Collective

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

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

718:   Level: intermediate

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

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

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

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

769:   Logically Collective

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

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

782:   Level: intermediate

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

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

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

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

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

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

837:   Logically Collective

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

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

847:   Level: intermediate

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

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

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

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

867:   Not Collective

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

876:   Level: beginner

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

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

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

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

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

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

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

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

921:   Not Collective

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

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

931:   Level: beginner

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

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

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

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

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

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

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

964:   Not Collective

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

973:   Level: intermediate

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

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

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

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

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

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

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

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

1016:   Not Collective

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

1025:   Level: intermediate

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

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

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

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

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

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

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

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

1071:   Not Collective

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

1080:   Level: intermediate

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

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

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

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

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

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

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

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

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

1155:   Collective

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

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

1165:   Level: intermediate

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

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

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

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

1189:   Collective

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

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

1199:   Level: intermediate

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

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

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

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

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

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

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

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

1260:   Logically Collective

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

1268:   Level: intermediate

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

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

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

1285:   Logically Collective

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

1294:   Level: intermediate

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

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

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

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

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

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

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

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

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

1357:   Collective

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

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

1367:   Level: advanced

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

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

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

1389:   PetscFunctionBegin;
1393:   PetscAssertPointer(Y, 3);

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

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

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

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

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

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

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

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

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

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

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

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

1509:   Collective

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

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

1518:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

1636:   Collective

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

1643:   Level: advanced

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

1652:   PetscFunctionBegin;
1655:   PetscCheckSameComm(X, 1, is, 2);
1656:   PetscAssertPointer(Y, 3);

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

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

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

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

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

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

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

1746:   Not Collective.

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

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

1754:   Level: beginner

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

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

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

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

1787:   Not Collective.

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

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

1795:   Level: beginner

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

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

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

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

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

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

1840:   Not Collective.

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

1846:   Level: beginner

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

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

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

1874:   Collective

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

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

1882:   Level: beginner

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

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

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

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

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

1922:   Logically Collective.

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

1928:   Level: beginner

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

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

1954:   Logically Collective

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

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

1962:   Level: beginner

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

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

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

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

1992:   Logically Collective

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

1998:   Level: beginner

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

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

2021:   Not Collective

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

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

2029:   Level: beginner

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

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

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

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

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

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

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

2071:   Not Collective

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

2077:   Level: beginner

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

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

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

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

2110:   Logically Collective

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

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

2118:   Level: intermediate

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

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

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

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

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

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

2152:   Logically Collective

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

2158:   Level: beginner

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

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

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

2185:   Logically Collective; No Fortran Support

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

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

2194:   Level: intermediate

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

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

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

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

2221:   Logically Collective; No Fortran Support

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

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

2234:   Level: intermediate

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

2243:   PetscFunctionBegin;
2244:   PetscAssertPointer(x, 1);
2246:   PetscAssertPointer(a, 3);

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

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

2258:   Logically Collective; No Fortran Support

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

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

2267:   Level: beginner

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

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

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

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

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

2305:   Logically Collective; No Fortran Support

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

2311:   Level: beginner

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

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

2338:   Not Collective; No Fortran Support

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

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

2347:   Level: beginner

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

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

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

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

2382:   Not Collective; No Fortran Support

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

2388:   Level: beginner

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

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

2413:   Logically Collective; No Fortran Support

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

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

2422:   Level: beginner

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

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

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

2453:   Logically Collective; No Fortran Support

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

2459:   Level: beginner

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

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

2488:   Logically Collective; No Fortran Support

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

2494:   Level: developer

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

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

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

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

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

2524:   Logically Collective; No Fortran Support

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

2530:   Level: developer

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

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

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

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

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

2558:     Collective

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

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

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

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

2583:     Level: beginner

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

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

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

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

2598:     Logically Collective

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

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

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

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

2619:     Level: beginner

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

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

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

2630:     Collective

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

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

2639:     Level: beginner

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

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

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

2653:     Logically Collective

2655:     Input Parameter:
2656: .   x - vector

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

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

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

2674:      Level: beginner

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

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

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

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

2691:     Logically Collective

2693:     Input Parameter:
2694: .   x - vector

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

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

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

2712:     Level: beginner

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

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

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

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

2727:     Logically Collective

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

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

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

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

2748:     Level: beginner

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

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

2758:   Logically Collective

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

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

2770:   Level: developer

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

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

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

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

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

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

2808:   Logically Collective

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

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

2820:   Level: developer

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

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

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

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

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

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

2856:   Logically Collective

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

2866:   Level: developer

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

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

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

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

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

2898:   Logically Collective

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

2908:   Level: developer

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

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

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

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

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

2941:   Logically Collective

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

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

2951:   Level: developer

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

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

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

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

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

2984:   Logically Collective

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

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

2994:   Level: developer

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

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

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

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

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

3025:   Logically Collective

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

3033:   Level: developer

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

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

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

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

3060:   Logically Collective

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

3068:   Level: developer

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

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

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

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

3097:   Logically Collective

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

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

3111:   Level: developer

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

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

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

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

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

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

3152:   Logically Collective

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

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

3166:   Level: developer

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

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

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

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

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

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

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

3206:   Logically Collective

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

3218:   Level: developer

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

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

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

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

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

3250:   Logically Collective

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

3262:   Level: developer

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

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

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

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

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

3296:   Logically Collective

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

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

3312:   Level: developer

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

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

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

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

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

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

3357:   Logically Collective

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

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

3373:   Level: developer

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

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

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

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

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

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

3416:   Logically Collective

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

3430:   Level: developer

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

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

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

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

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

3462:   Logically Collective

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

3476:   Level: developer

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

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

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

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

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

3510:   Logically Collective

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

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

3522:   Level: developer

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

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

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

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

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

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

3558:   Logically Collective

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

3568:   Level: developer

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

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

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

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

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

3602:   Logically Collective

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

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

3612:   Level: developer

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

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

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

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

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

3643:   Logically Collective

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

3651:   Level: developer

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

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

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

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

3680:   Logically Collective

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

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

3694:   Level: developer

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

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

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

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

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

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

3734:   Logically Collective

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

3746:   Level: developer

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

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

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

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

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

3780:   Logically Collective

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

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

3796:   Level: beginner

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

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

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

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

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

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

3840:   Logically Collective

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

3854:   Level: beginner

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

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

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

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

3883: #if defined(PETSC_USE_DEBUG)

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

3888:   Logically Collective

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

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

3897:   Level: advanced

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

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

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

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

3936:   Logically Collective

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

3941:   Level: intermediate

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

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

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

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

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

3982:   Logically Collective

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

3987:   Level: intermediate

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

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

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

4009:   Logically Collective

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

4015:   Level: intermediate

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

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

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

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

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

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

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

4057:   Level: deprecated

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

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

4072:   Level: deprecated

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

4083: #endif