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:   /* Cached data? */
209:   PetscCall(VecNormAvailable(x, type, &flg, val));
210:   if (flg) PetscFunctionReturn(PETSC_SUCCESS);

212:   PetscCall(VecLockReadPush(x));
213:   PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
214:   PetscUseTypeMethod(x, norm, type, val);
215:   PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
216:   PetscCall(VecLockReadPop(x));

218:   if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
219:   PetscFunctionReturn(PETSC_SUCCESS);
220: }

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

225:   Not Collective

227:   Input Parameters:
228: + x    - the vector
229: - 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
230:           `NORM_1_AND_2`, which computes both norms and stores them
231:           in a two element array.

233:   Output Parameters:
234: + available - `PETSC_TRUE` if the val returned is valid
235: - val       - the norm

237:   Level: intermediate

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

244: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
245:           `VecNormBegin()`, `VecNormEnd()`
246: @*/
247: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
248: {
249:   PetscFunctionBegin;
252:   PetscAssertPointer(available, 3);
253:   PetscAssertPointer(val, 4);

255:   if (type == NORM_1_AND_2) {
256:     *available = PETSC_FALSE;
257:   } else {
258:     PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
259:   }
260:   PetscFunctionReturn(PETSC_SUCCESS);
261: }

263: /*@
264:   VecNormalize - Normalizes a vector by its 2-norm.

266:   Collective

268:   Input Parameter:
269: . x - the vector

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

274:   Level: intermediate

276: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
277: @*/
278: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
279: {
280:   PetscReal norm;

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

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

303:   Collective

305:   Input Parameter:
306: . x - the vector

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

312:   Level: intermediate

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

317:   Returns the smallest index with the maximum value

319: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
320: @*/
321: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
322: {
323:   PetscFunctionBegin;
326:   VecCheckAssembled(x);
327:   if (p) PetscAssertPointer(p, 2);
328:   PetscAssertPointer(val, 3);
329:   PetscCall(VecLockReadPush(x));
330:   PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
331:   PetscUseTypeMethod(x, max, p, val);
332:   PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
333:   PetscCall(VecLockReadPop(x));
334:   PetscFunctionReturn(PETSC_SUCCESS);
335: }

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

340:   Collective

342:   Input Parameter:
343: . x - the vector

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

349:   Level: intermediate

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

354:   This returns the smallest index with the minimum value

356: .seealso: [](ch_vectors), `Vec`, `VecMax()`
357: @*/
358: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
359: {
360:   PetscFunctionBegin;
363:   VecCheckAssembled(x);
364:   if (p) PetscAssertPointer(p, 2);
365:   PetscAssertPointer(val, 3);
366:   PetscCall(VecLockReadPush(x));
367:   PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
368:   PetscUseTypeMethod(x, min, p, val);
369:   PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
370:   PetscCall(VecLockReadPop(x));
371:   PetscFunctionReturn(PETSC_SUCCESS);
372: }

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

378:   Collective

380:   Input Parameters:
381: + x - first vector
382: - y - second vector

384:   Output Parameter:
385: . val - the dot product

387:   Level: intermediate

389:   Notes for Users of Complex Numbers:
390:   For complex vectors, `VecTDot()` computes the indefinite form
391: $     val = (x,y) = y^T x,
392:   where y^T denotes the transpose of y.

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

398: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
399: @*/
400: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
401: {
402:   PetscFunctionBegin;
405:   PetscAssertPointer(val, 3);
408:   PetscCheckSameTypeAndComm(x, 1, y, 2);
409:   VecCheckSameSize(x, 1, y, 2);
410:   VecCheckAssembled(x);
411:   VecCheckAssembled(y);

413:   PetscCall(VecLockReadPush(x));
414:   PetscCall(VecLockReadPush(y));
415:   PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
416:   PetscUseTypeMethod(x, tdot, y, val);
417:   PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
418:   PetscCall(VecLockReadPop(x));
419:   PetscCall(VecLockReadPop(y));
420:   PetscFunctionReturn(PETSC_SUCCESS);
421: }

423: PetscErrorCode VecScaleAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
424: {
425:   PetscReal   norms[4];
426:   PetscBool   flgs[4];
427:   PetscScalar one = 1.0;

429:   PetscFunctionBegin;
432:   VecCheckAssembled(x);
433:   PetscCall(VecSetErrorIfLocked(x, 1));
434:   if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);

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

439:   PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
440:   VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
441:   PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));

443:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
444:   /* put the scaled stashed norms back into the Vec */
445:   for (PetscInt i = 0; i < 4; i++) {
446:     PetscReal ar = PetscAbsScalar(alpha);
447:     if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
448:   }
449:   PetscFunctionReturn(PETSC_SUCCESS);
450: }

452: /*@
453:   VecScale - Scales a vector.

455:   Not Collective

457:   Input Parameters:
458: + x     - the vector
459: - alpha - the scalar

461:   Level: intermediate

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

466: .seealso: [](ch_vectors), `Vec`, `VecSet()`
467: @*/
468: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
469: {
470:   PetscFunctionBegin;
471:   PetscCall(VecScaleAsync_Private(x, alpha, NULL));
472:   PetscFunctionReturn(PETSC_SUCCESS);
473: }

475: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
476: {
477:   PetscFunctionBegin;
480:   VecCheckAssembled(x);
482:   PetscCall(VecSetErrorIfLocked(x, 1));

484:   if (alpha == 0) {
485:     PetscReal norm;
486:     PetscBool set;

488:     PetscCall(VecNormAvailable(x, NORM_2, &set, &norm));
489:     if (set == PETSC_TRUE && norm == 0) PetscFunctionReturn(PETSC_SUCCESS);
490:   }
491:   PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
492:   VecMethodDispatch(x, dctx, VecAsyncFnName(Set), set, (Vec, PetscScalar, PetscDeviceContext), alpha);
493:   PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
494:   PetscCall(PetscObjectStateIncrease((PetscObject)x));

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

498:   {
499:     PetscReal      val = PetscAbsScalar(alpha);
500:     const PetscInt N   = x->map->N;

502:     if (N == 0) {
503:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0l));
504:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
505:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
506:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
507:     } else if (val > PETSC_MAX_REAL / N) {
508:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
509:     } else {
510:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
511:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
512:       val *= PetscSqrtReal((PetscReal)N);
513:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
514:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
515:     }
516:   }
517:   PetscFunctionReturn(PETSC_SUCCESS);
518: }

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

523:   Logically Collective

525:   Input Parameters:
526: + x     - the vector
527: - alpha - the scalar

529:   Level: beginner

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

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

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

542: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
543: @*/
544: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
545: {
546:   PetscFunctionBegin;
547:   PetscCall(VecSetAsync_Private(x, alpha, NULL));
548:   PetscFunctionReturn(PETSC_SUCCESS);
549: }

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

580:   Logically Collective

582:   Input Parameters:
583: + alpha - the scalar
584: . x     - vector scale by `alpha`
585: - y     - vector accumulated into

587:   Output Parameter:
588: . y - output vector

590:   Level: intermediate

592:   Notes:
593:   This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
594: .vb
595:     VecAXPY(y,alpha,x)                   y = alpha x           +      y
596:     VecAYPX(y,beta,x)                    y =       x           + beta y
597:     VecAXPBY(y,alpha,beta,x)             y = alpha x           + beta y
598:     VecWAXPY(w,alpha,x,y)                w = alpha x           +      y
599:     VecAXPBYPCZ(z,alpha,beta,gamma,x,y)  z = alpha x           + beta y + gamma z
600:     VecMAXPY(y,nv,alpha[],x[])           y = sum alpha[i] x[i] +      y
601: .ve

603: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
604: @*/
605: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
606: {
607:   PetscFunctionBegin;
608:   PetscCall(VecAXPYAsync_Private(y, alpha, x, NULL));
609:   PetscFunctionReturn(PETSC_SUCCESS);
610: }

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

642: /*@
643:   VecAYPX - Computes `y = x + beta y`.

645:   Logically Collective

647:   Input Parameters:
648: + beta - the scalar
649: . x    - the unscaled vector
650: - y    - the vector to be scaled

652:   Output Parameter:
653: . y - output vector

655:   Level: intermediate

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

660: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
661: @*/
662: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
663: {
664:   PetscFunctionBegin;
665:   PetscCall(VecAYPXAsync_Private(y, beta, x, NULL));
666:   PetscFunctionReturn(PETSC_SUCCESS);
667: }

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

688:   PetscCall(VecSetErrorIfLocked(y, 1));
689:   PetscCall(VecLockReadPush(x));
690:   PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
691:   VecMethodDispatch(y, dctx, VecAsyncFnName(AXPBY), axpby, (Vec, PetscScalar, PetscScalar, Vec, PetscDeviceContext), alpha, beta, x);
692:   PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
693:   PetscCall(PetscObjectStateIncrease((PetscObject)y));
694:   PetscCall(VecLockReadPop(x));
695:   PetscFunctionReturn(PETSC_SUCCESS);
696: }
697: /*@
698:   VecAXPBY - Computes `y = alpha x + beta y`.

700:   Logically Collective

702:   Input Parameters:
703: + alpha - first scalar
704: . beta  - second scalar
705: . x     - the first scaled vector
706: - y     - the second scaled vector

708:   Output Parameter:
709: . y - output vector

711:   Level: intermediate

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

716: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
717: @*/
718: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
719: {
720:   PetscFunctionBegin;
721:   PetscCall(VecAXPBYAsync_Private(y, alpha, beta, x, NULL));
722:   PetscFunctionReturn(PETSC_SUCCESS);
723: }

725: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
726: {
727:   PetscFunctionBegin;
734:   PetscCheckSameTypeAndComm(x, 5, y, 6);
735:   PetscCheckSameTypeAndComm(x, 5, z, 1);
736:   VecCheckSameSize(x, 1, y, 5);
737:   VecCheckSameSize(x, 1, z, 6);
738:   PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
739:   PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
740:   VecCheckAssembled(x);
741:   VecCheckAssembled(y);
742:   VecCheckAssembled(z);
746:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);

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

762:   Logically Collective

764:   Input Parameters:
765: + alpha - first scalar
766: . beta  - second scalar
767: . gamma - third scalar
768: . x     - first vector
769: . y     - second vector
770: - z     - third vector

772:   Output Parameter:
773: . z - output vector

775:   Level: intermediate

777:   Note:
778:   `x`, `y` and `z` must be different vectors

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

783: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
784: @*/
785: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
786: {
787:   PetscFunctionBegin;
788:   PetscCall(VecAXPBYPCZAsync_Private(z, alpha, beta, gamma, x, y, NULL));
789:   PetscFunctionReturn(PETSC_SUCCESS);
790: }

792: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
793: {
794:   PetscFunctionBegin;
801:   PetscCheckSameTypeAndComm(x, 3, y, 4);
802:   PetscCheckSameTypeAndComm(y, 4, w, 1);
803:   VecCheckSameSize(x, 3, y, 4);
804:   VecCheckSameSize(x, 3, w, 1);
805:   PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
806:   PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
807:   VecCheckAssembled(x);
808:   VecCheckAssembled(y);
810:   PetscCall(VecSetErrorIfLocked(w, 1));

812:   PetscCall(VecLockReadPush(x));
813:   PetscCall(VecLockReadPush(y));
814:   if (alpha == (PetscScalar)0.0) {
815:     PetscCall(VecCopyAsync_Private(y, w, dctx));
816:   } else {
817:     PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
818:     VecMethodDispatch(w, dctx, VecAsyncFnName(WAXPY), waxpy, (Vec, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, x, y);
819:     PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
820:     PetscCall(PetscObjectStateIncrease((PetscObject)w));
821:   }
822:   PetscCall(VecLockReadPop(x));
823:   PetscCall(VecLockReadPop(y));
824:   PetscFunctionReturn(PETSC_SUCCESS);
825: }

827: /*@
828:   VecWAXPY - Computes `w = alpha x + y`.

830:   Logically Collective

832:   Input Parameters:
833: + alpha - the scalar
834: . x     - first vector, multiplied by `alpha`
835: - y     - second vector

837:   Output Parameter:
838: . w - the result

840:   Level: intermediate

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

845:   Developer Notes:
846:   The implementation is optimized for alpha of -1.0, 0.0, and 1.0

848: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
849: @*/
850: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
851: {
852:   PetscFunctionBegin;
853:   PetscCall(VecWAXPYAsync_Private(w, alpha, x, y, NULL));
854:   PetscFunctionReturn(PETSC_SUCCESS);
855: }

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

860:   Not Collective

862:   Input Parameters:
863: + x    - vector to insert in
864: . ni   - number of elements to add
865: . ix   - indices where to add
866: . y    - array of values
867: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries

869:   Level: beginner

871:   Notes:
872: .vb
873:    `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
874: .ve

876:   Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
877:   options cannot be mixed without intervening calls to the assembly
878:   routines.

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

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

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

891: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
892:           `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
893: @*/
894: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
895: {
896:   PetscFunctionBeginHot;
898:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
899:   PetscAssertPointer(ix, 3);
900:   PetscAssertPointer(y, 4);

903:   PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
904:   PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
905:   PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
906:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
907:   PetscFunctionReturn(PETSC_SUCCESS);
908: }

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

914:   Not Collective

916:   Input Parameters:
917: + x  - vector to get values from
918: . ni - number of elements to get
919: - ix - indices where to get them from (in global 1d numbering)

921:   Output Parameter:
922: . y - array of values

924:   Level: beginner

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

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

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

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

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

939: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
940: @*/
941: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
942: {
943:   PetscFunctionBegin;
945:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
946:   PetscAssertPointer(ix, 3);
947:   PetscAssertPointer(y, 4);
949:   VecCheckAssembled(x);
950:   PetscUseTypeMethod(x, getvalues, ni, ix, y);
951:   PetscFunctionReturn(PETSC_SUCCESS);
952: }

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

957:   Not Collective

959:   Input Parameters:
960: + x    - vector to insert in
961: . ni   - number of blocks to add
962: . ix   - indices where to add in block count, rather than element count
963: . y    - array of values
964: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries

966:   Level: intermediate

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

972:   Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
973:   options cannot be mixed without intervening calls to the assembly
974:   routines.

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

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

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

986: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
987:           `VecSetValues()`
988: @*/
989: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
990: {
991:   PetscFunctionBeginHot;
993:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
994:   PetscAssertPointer(ix, 3);
995:   PetscAssertPointer(y, 4);

998:   PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
999:   PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1000:   PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1001:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
1002:   PetscFunctionReturn(PETSC_SUCCESS);
1003: }

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

1009:   Not Collective

1011:   Input Parameters:
1012: + x    - vector to insert in
1013: . ni   - number of elements to add
1014: . ix   - indices where to add
1015: . y    - array of values
1016: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries

1018:   Level: intermediate

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

1023:   Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
1024:   options cannot be mixed without intervening calls to the assembly
1025:   routines.

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

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

1032: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1033:           `VecSetValuesBlockedLocal()`
1034: @*/
1035: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1036: {
1037:   PetscInt lixp[128], *lix = lixp;

1039:   PetscFunctionBeginHot;
1041:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1042:   PetscAssertPointer(ix, 3);
1043:   PetscAssertPointer(y, 4);

1046:   PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1047:   if (!x->ops->setvalueslocal) {
1048:     if (x->map->mapping) {
1049:       if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1050:       PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1051:       PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1052:       if (ni > 128) PetscCall(PetscFree(lix));
1053:     } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1054:   } else PetscUseTypeMethod(x, setvalueslocal, ni, ix, y, iora);
1055:   PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1056:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
1057:   PetscFunctionReturn(PETSC_SUCCESS);
1058: }

1060: /*@
1061:   VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1062:   using a local ordering of the nodes.

1064:   Not Collective

1066:   Input Parameters:
1067: + x    - vector to insert in
1068: . ni   - number of blocks to add
1069: . ix   - indices where to add in block count, not element count
1070: . y    - array of values
1071: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries

1073:   Level: intermediate

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

1079:   Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1080:   options cannot be mixed without intervening calls to the assembly
1081:   routines.

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

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

1088: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1089:           `VecSetLocalToGlobalMapping()`
1090: @*/
1091: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1092: {
1093:   PetscInt lixp[128], *lix = lixp;

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

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

1135:   PetscCall(VecLockReadPush(x));
1136:   PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1137:   PetscCall((*mxdot)(x, nv, y, result));
1138:   PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1139:   PetscCall(VecLockReadPop(x));
1140:   for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1141:   PetscFunctionReturn(PETSC_SUCCESS);
1142: }

1144: /*@
1145:   VecMTDot - Computes indefinite vector multiple dot products.
1146:   That is, it does NOT use the complex conjugate.

1148:   Collective

1150:   Input Parameters:
1151: + x  - one vector
1152: . nv - number of vectors
1153: - y  - array of vectors.  Note that vectors are pointers

1155:   Output Parameter:
1156: . val - array of the dot products

1158:   Level: intermediate

1160:   Notes for Users of Complex Numbers:
1161:   For complex vectors, `VecMTDot()` computes the indefinite form
1162: $      val = (x,y) = y^T x,
1163:   where y^T denotes the transpose of y.

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

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

1179: /*@
1180:   VecMDot - Computes multiple vector dot products.

1182:   Collective

1184:   Input Parameters:
1185: + x  - one vector
1186: . nv - number of vectors
1187: - y  - array of vectors.

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

1192:   Level: intermediate

1194:   Notes for Users of Complex Numbers:
1195:   For complex vectors, `VecMDot()` computes
1196: $     val = (x,y) = y^H x,
1197:   where y^H denotes the conjugate transpose of y.

1199:   Use `VecMTDot()` for the indefinite form
1200: $     val = (x,y) = y^T x,
1201:   where y^T denotes the transpose of y.

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

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

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

1238:     if (zeros < nv) {
1239:       PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1240:       VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1241:       PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1242:       PetscCall(PetscObjectStateIncrease((PetscObject)y));
1243:     }

1245:     for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1246:   }
1247:   PetscFunctionReturn(PETSC_SUCCESS);
1248: }

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

1253:   Logically Collective

1255:   Input Parameters:
1256: + nv    - number of scalars and x-vectors
1257: . alpha - array of scalars
1258: . y     - one vector
1259: - x     - array of vectors

1261:   Level: intermediate

1263:   Note:
1264:   `y` cannot be any of the `x` vectors

1266: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`,`VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1267: @*/
1268: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1269: {
1270:   PetscFunctionBegin;
1271:   PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1272:   PetscFunctionReturn(PETSC_SUCCESS);
1273: }

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

1278:   Logically Collective

1280:   Input Parameters:
1281: + nv    - number of scalars and x-vectors
1282: . alpha - array of scalars
1283: . beta  - scalar
1284: . y     - one vector
1285: - x     - array of vectors

1287:   Level: intermediate

1289:   Note:
1290:   `y` cannot be any of the `x` vectors.

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

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

1307:   if (y->ops->maxpby) {
1308:     PetscInt zeros = 0;

1310:     if (nv) {
1311:       PetscAssertPointer(alpha, 3);
1312:       PetscAssertPointer(x, 5);
1313:     }

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

1327:     if (zeros < nv) { // has nonzero alpha
1328:       PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1329:       PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1330:       PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1331:       PetscCall(PetscObjectStateIncrease((PetscObject)y));
1332:     } else {
1333:       PetscCall(VecScale(y, beta));
1334:     }

1336:     for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1337:   } else { // no maxpby
1338:     if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1339:     else PetscCall(VecScale(y, beta));
1340:     PetscCall(VecMAXPY(y, nv, alpha, x));
1341:   }
1342:   PetscFunctionReturn(PETSC_SUCCESS);
1343: }

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

1350:   Collective

1352:   Input Parameters:
1353: + nx - number of vectors to be concatenated
1354: - X  - array containing the vectors to be concatenated in the order of concatenation

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

1360:   Level: advanced

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

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

1372: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1373: @*/
1374: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1375: {
1376:   MPI_Comm comm;
1377:   VecType  vec_type;
1378:   Vec      Ytmp, Xtmp;
1379:   IS      *is_tmp;
1380:   PetscInt i, shift = 0, Xnl, Xng, Xbegin;

1382:   PetscFunctionBegin;
1386:   PetscAssertPointer(Y, 3);

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

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

1428:     Input Parameters:
1429: +   X - the original vector
1430: -   is - the index set of the subvector

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

1437: */
1438: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1439: {
1440:   PetscInt  gstart, gend, lstart;
1441:   PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1442:   PetscInt  n, N, ibs, vbs, bs = -1;

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

1461:   *contig    = red[0];
1462:   *start     = lstart;
1463:   *blocksize = bs;
1464:   PetscFunctionReturn(PETSC_SUCCESS);
1465: }

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

1469:     Input Parameters:
1470: +   X - the original vector
1471: .   is - the index set of the subvector
1472: -   bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()

1474:     Output Parameter:
1475: .   Z  - the subvector, which will compose the VecScatter context on output
1476: */
1477: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1478: {
1479:   PetscInt   n, N;
1480:   VecScatter vscat;
1481:   Vec        Y;

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

1499: /*@
1500:   VecGetSubVector - Gets a vector representing part of another vector

1502:   Collective

1504:   Input Parameters:
1505: + X  - vector from which to extract a subvector
1506: - is - index set representing portion of `X` to extract

1508:   Output Parameter:
1509: . Y - subvector corresponding to `is`

1511:   Level: advanced

1513:   Notes:
1514:   The subvector `Y` should be returned with `VecRestoreSubVector()`.
1515:   `X` and must be defined on the same communicator

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

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

1522: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1523: @*/
1524: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1525: {
1526:   Vec Z;

1528:   PetscFunctionBegin;
1531:   PetscCheckSameComm(X, 1, is, 2);
1532:   PetscAssertPointer(Y, 3);
1533:   if (X->ops->getsubvector) {
1534:     PetscUseTypeMethod(X, getsubvector, is, &Z);
1535:   } else { /* Default implementation currently does no caching */
1536:     PetscBool contig;
1537:     PetscInt  n, N, start, bs;

1539:     PetscCall(ISGetLocalSize(is, &n));
1540:     PetscCall(ISGetSize(is, &N));
1541:     PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1542:     if (contig) { /* We can do a no-copy implementation */
1543:       const PetscScalar *x;
1544:       PetscInt           state = 0;
1545:       PetscBool          isstd, iscuda, iship;

1547:       PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1548:       PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1549:       PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1550:       if (iscuda) {
1551: #if defined(PETSC_HAVE_CUDA)
1552:         const PetscScalar *x_d;
1553:         PetscMPIInt        size;
1554:         PetscOffloadMask   flg;

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

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

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

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

1626: /*@
1627:   VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`

1629:   Collective

1631:   Input Parameters:
1632: + X  - vector from which subvector was obtained
1633: . is - index set representing the subset of `X`
1634: - Y  - subvector being restored

1636:   Level: advanced

1638: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1639: @*/
1640: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1641: {
1642:   PETSC_UNUSED PetscObjectState dummystate = 0;
1643:   PetscBool                     unchanged;

1645:   PetscFunctionBegin;
1648:   PetscCheckSameComm(X, 1, is, 2);
1649:   PetscAssertPointer(Y, 3);

1652:   if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1653:   else {
1654:     PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1655:     if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1656:       VecScatter scatter;
1657:       PetscInt   state;

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

1662:       PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1663:       if (scatter) {
1664:         PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1665:         PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1666:       } else {
1667:         PetscBool iscuda, iship;
1668:         PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1669:         PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));

1671:         if (iscuda) {
1672: #if defined(PETSC_HAVE_CUDA)
1673:           PetscOffloadMask ymask = (*Y)->offloadmask;

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

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

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

1739:   Not Collective.

1741:   Input Parameter:
1742: . v - The vector for which the local vector is desired.

1744:   Output Parameter:
1745: . w - Upon exit this contains the local vector.

1747:   Level: beginner

1749: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1750: @*/
1751: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1752: {
1753:   PetscMPIInt size;

1755:   PetscFunctionBegin;
1757:   PetscAssertPointer(w, 2);
1758:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1759:   if (size == 1) PetscCall(VecDuplicate(v, w));
1760:   else if (v->ops->createlocalvector) PetscUseTypeMethod(v, createlocalvector, w);
1761:   else {
1762:     VecType  type;
1763:     PetscInt n;

1765:     PetscCall(VecCreate(PETSC_COMM_SELF, w));
1766:     PetscCall(VecGetLocalSize(v, &n));
1767:     PetscCall(VecSetSizes(*w, n, n));
1768:     PetscCall(VecGetBlockSize(v, &n));
1769:     PetscCall(VecSetBlockSize(*w, n));
1770:     PetscCall(VecGetType(v, &type));
1771:     PetscCall(VecSetType(*w, type));
1772:   }
1773:   PetscFunctionReturn(PETSC_SUCCESS);
1774: }

1776: /*@
1777:   VecGetLocalVectorRead - Maps the local portion of a vector into a
1778:   vector.

1780:   Not Collective.

1782:   Input Parameter:
1783: . v - The vector for which the local vector is desired.

1785:   Output Parameter:
1786: . w - Upon exit this contains the local vector.

1788:   Level: beginner

1790:   Notes:
1791:   You must call `VecRestoreLocalVectorRead()` when the local
1792:   vector is no longer needed.

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

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

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

1820:     PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1821:     PetscCall(VecPlaceArray(w, a));
1822:   }
1823:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
1824:   PetscCall(VecLockReadPush(v));
1825:   PetscCall(VecLockReadPush(w));
1826:   PetscFunctionReturn(PETSC_SUCCESS);
1827: }

1829: /*@
1830:   VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1831:   previously mapped into a vector using `VecGetLocalVectorRead()`.

1833:   Not Collective.

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

1839:   Level: beginner

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

1853:     PetscCall(VecGetArrayRead(w, &a));
1854:     PetscCall(VecRestoreArrayRead(v, &a));
1855:     PetscCall(VecResetArray(w));
1856:   }
1857:   PetscCall(VecLockReadPop(v));
1858:   PetscCall(VecLockReadPop(w));
1859:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
1860:   PetscFunctionReturn(PETSC_SUCCESS);
1861: }

1863: /*@
1864:   VecGetLocalVector - Maps the local portion of a vector into a
1865:   vector.

1867:   Collective

1869:   Input Parameter:
1870: . v - The vector for which the local vector is desired.

1872:   Output Parameter:
1873: . w - Upon exit this contains the local vector.

1875:   Level: beginner

1877:   Notes:
1878:   You must call `VecRestoreLocalVector()` when the local
1879:   vector is no longer needed.

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

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

1904:     PetscCall(VecGetArray(v, &a));
1905:     PetscCall(VecPlaceArray(w, a));
1906:   }
1907:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
1908:   PetscFunctionReturn(PETSC_SUCCESS);
1909: }

1911: /*@
1912:   VecRestoreLocalVector - Unmaps the local portion of a vector
1913:   previously mapped into a vector using `VecGetLocalVector()`.

1915:   Logically Collective.

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

1921:   Level: beginner

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

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

1947:   Logically Collective

1949:   Input Parameter:
1950: . x - the vector

1952:   Output Parameter:
1953: . a - location to put pointer to the array

1955:   Level: beginner

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

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

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

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

1985:   Logically Collective

1987:   Input Parameters:
1988: + x - the vector
1989: - a - location of pointer to array obtained from `VecGetArray()`

1991:   Level: beginner

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

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

2014:   Not Collective

2016:   Input Parameter:
2017: . x - the vector

2019:   Output Parameter:
2020: . a - the array

2022:   Level: beginner

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

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

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

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

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

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

2061: /*@C
2062:   VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`

2064:   Not Collective

2066:   Input Parameters:
2067: + x - the vector
2068: - a - the array

2070:   Level: beginner

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

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

2089:     // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2090:     // we can just undo that
2091:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2092:     PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2093:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
2094:   }
2095:   if (a) *a = NULL;
2096:   PetscFunctionReturn(PETSC_SUCCESS);
2097: }

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

2103:   Logically Collective

2105:   Input Parameter:
2106: . x - the vector

2108:   Output Parameter:
2109: . a - location to put pointer to the array

2111:   Level: intermediate

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

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

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

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

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

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

2145:   Logically Collective

2147:   Input Parameters:
2148: + x - the vector
2149: - a - location of pointer to array obtained from `VecGetArray()`

2151:   Level: beginner

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

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

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

2178:   Logically Collective; No Fortran Support

2180:   Input Parameters:
2181: + x - the vectors
2182: - n - the number of vectors

2184:   Output Parameter:
2185: . a - location to put pointer to the array

2187:   Level: intermediate

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

2192: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2193: @*/
2194: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2195: {
2196:   PetscInt      i;
2197:   PetscScalar **q;

2199:   PetscFunctionBegin;
2200:   PetscAssertPointer(x, 1);
2202:   PetscAssertPointer(a, 3);
2203:   PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2204:   PetscCall(PetscMalloc1(n, &q));
2205:   for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2206:   *a = q;
2207:   PetscFunctionReturn(PETSC_SUCCESS);
2208: }

2210: /*@C
2211:   VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2212:   has been called.

2214:   Logically Collective; No Fortran Support

2216:   Input Parameters:
2217: + x - the vector
2218: . n - the number of vectors
2219: - a - location of pointer to arrays obtained from `VecGetArrays()`

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

2227:   Level: intermediate

2229: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2230: @*/
2231: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2232: {
2233:   PetscInt      i;
2234:   PetscScalar **q = *a;

2236:   PetscFunctionBegin;
2237:   PetscAssertPointer(x, 1);
2239:   PetscAssertPointer(a, 3);

2241:   for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2242:   PetscCall(PetscFree(q));
2243:   PetscFunctionReturn(PETSC_SUCCESS);
2244: }

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

2251:   Logically Collective; No Fortran Support

2253:   Input Parameter:
2254: . x - the vector

2256:   Output Parameters:
2257: + a     - location to put pointer to the array
2258: - mtype - memory type of the array

2260:   Level: beginner

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

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

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

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

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

2298:   Logically Collective; No Fortran Support

2300:   Input Parameters:
2301: + x - the vector
2302: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`

2304:   Level: beginner

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

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

2331:   Not Collective; No Fortran Support

2333:   Input Parameter:
2334: . x - the vector

2336:   Output Parameters:
2337: + a     - the array
2338: - mtype - memory type of the array

2340:   Level: beginner

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

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

2361:     // see VecGetArrayRead() for why
2362:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2363:     PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2364:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
2365:   } else {
2366:     PetscCall(VecGetArrayRead(x, a));
2367:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2368:   }
2369:   PetscFunctionReturn(PETSC_SUCCESS);
2370: }

2372: /*@C
2373:   VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`

2375:   Not Collective; No Fortran Support

2377:   Input Parameters:
2378: + x - the vector
2379: - a - the array

2381:   Level: beginner

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

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

2406:   Not Collective; No Fortran Support

2408:   Input Parameter:
2409: . x - the vector

2411:   Output Parameters:
2412: + a     - the array
2413: - mtype - memory type of the array

2415:   Level: beginner

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

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

2443: /*@C
2444:   VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`

2446:   Not Collective; No Fortran Support

2448:   Input Parameters:
2449: + x - the vector
2450: - a - the array

2452:   Level: beginner

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

2476: /*@
2477:   VecPlaceArray - Allows one to replace the array in a vector with an
2478:   array provided by the user. This is useful to avoid copying an array
2479:   into a vector.

2481:   Not Collective; No Fortran Support

2483:   Input Parameters:
2484: + vec   - the vector
2485: - array - the array

2487:   Level: developer

2489:   Notes:
2490:   Use `VecReplaceArray()` instead to permanently replace the array

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

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

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

2512: /*@C
2513:   VecReplaceArray - Allows one to replace the array in a vector with an
2514:   array provided by the user. This is useful to avoid copying an array
2515:   into a vector.

2517:   Not Collective; No Fortran Support

2519:   Input Parameters:
2520: + vec   - the vector
2521: - array - the array

2523:   Level: developer

2525:   Notes:
2526:   This permanently replaces the array and frees the memory associated
2527:   with the old array. Use `VecPlaceArray()` to temporarily replace the array.

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

2532: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2533: @*/
2534: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2535: {
2536:   PetscFunctionBegin;
2539:   PetscUseTypeMethod(vec, replacearray, array);
2540:   PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2541:   PetscFunctionReturn(PETSC_SUCCESS);
2542: }

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

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

2551:     Collective

2553:     Input Parameters:
2554: +   x - a vector to mimic
2555: -   n - the number of vectors to obtain

2557:     Output Parameters:
2558: +   y - Fortran pointer to the array of vectors
2559: -   ierr - error code

2561:     Example of Usage:
2562: .vb
2563: #include <petsc/finclude/petscvec.h>
2564:     use petscvec

2566:     Vec x
2567:     Vec, pointer :: y(:)
2568:     ....
2569:     call VecDuplicateVecsF90(x,2,y,ierr)
2570:     call VecSet(y(2),alpha,ierr)
2571:     call VecSet(y(2),alpha,ierr)
2572:     ....
2573:     call VecDestroyVecsF90(2,y,ierr)
2574: .ve

2576:     Level: beginner

2578:     Note:
2579:     Use `VecDestroyVecsF90()` to free the space.

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

2584: /*MC
2585:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2586:     `VecGetArrayF90()`.

2588:     Synopsis:
2589:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2591:     Logically Collective

2593:     Input Parameters:
2594: +   x - vector
2595: -   xx_v - the Fortran pointer to the array

2597:     Output Parameter:
2598: .   ierr - error code

2600:     Example of Usage:
2601: .vb
2602: #include <petsc/finclude/petscvec.h>
2603:     use petscvec

2605:     PetscScalar, pointer :: xx_v(:)
2606:     ....
2607:     call VecGetArrayF90(x,xx_v,ierr)
2608:     xx_v(3) = a
2609:     call VecRestoreArrayF90(x,xx_v,ierr)
2610: .ve

2612:     Level: beginner

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

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

2620:     Synopsis:
2621:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

2623:     Collective

2625:     Input Parameters:
2626: +   n - the number of vectors previously obtained
2627: -   x - pointer to array of vector pointers

2629:     Output Parameter:
2630: .   ierr - error code

2632:     Level: beginner

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

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

2643:     Synopsis:
2644:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2646:     Logically Collective

2648:     Input Parameter:
2649: .   x - vector

2651:     Output Parameters:
2652: +   xx_v - the Fortran pointer to the array
2653: -   ierr - error code

2655:     Example of Usage:
2656: .vb
2657: #include <petsc/finclude/petscvec.h>
2658:     use petscvec

2660:     PetscScalar, pointer :: xx_v(:)
2661:     ....
2662:     call VecGetArrayF90(x,xx_v,ierr)
2663:     xx_v(3) = a
2664:     call VecRestoreArrayF90(x,xx_v,ierr)
2665: .ve

2667:      Level: beginner

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

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

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

2681:     Synopsis:
2682:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2684:     Logically Collective

2686:     Input Parameter:
2687: .   x - vector

2689:     Output Parameters:
2690: +   xx_v - the Fortran pointer to the array
2691: -   ierr - error code

2693:     Example of Usage:
2694: .vb
2695: #include <petsc/finclude/petscvec.h>
2696:     use petscvec

2698:     PetscScalar, pointer :: xx_v(:)
2699:     ....
2700:     call VecGetArrayReadF90(x,xx_v,ierr)
2701:     a = xx_v(3)
2702:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2703: .ve

2705:     Level: beginner

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

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

2713: /*MC
2714:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2715:     `VecGetArrayReadF90()`.

2717:     Synopsis:
2718:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2720:     Logically Collective

2722:     Input Parameters:
2723: +   x - vector
2724: -   xx_v - the Fortran pointer to the array

2726:     Output Parameter:
2727: .   ierr - error code

2729:     Example of Usage:
2730: .vb
2731: #include <petsc/finclude/petscvec.h>
2732:     use petscvec

2734:     PetscScalar, pointer :: xx_v(:)
2735:     ....
2736:     call VecGetArrayReadF90(x,xx_v,ierr)
2737:     a = xx_v(3)
2738:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2739: .ve

2741:     Level: beginner

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

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

2751:   Logically Collective

2753:   Input Parameters:
2754: + x      - the vector
2755: . m      - first dimension of two dimensional array
2756: . n      - second dimension of two dimensional array
2757: . mstart - first index you will use in first coordinate direction (often 0)
2758: - nstart - first index in the second coordinate direction (often 0)

2760:   Output Parameter:
2761: . a - location to put pointer to the array

2763:   Level: developer

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

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

2773: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2774:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2775:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2776: @*/
2777: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2778: {
2779:   PetscInt     i, N;
2780:   PetscScalar *aa;

2782:   PetscFunctionBegin;
2784:   PetscAssertPointer(a, 6);
2786:   PetscCall(VecGetLocalSize(x, &N));
2787:   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);
2788:   PetscCall(VecGetArray(x, &aa));

2790:   PetscCall(PetscMalloc1(m, a));
2791:   for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2792:   *a -= mstart;
2793:   PetscFunctionReturn(PETSC_SUCCESS);
2794: }

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

2801:   Logically Collective

2803:   Input Parameters:
2804: + x      - the vector
2805: . m      - first dimension of two dimensional array
2806: . n      - second dimension of two dimensional array
2807: . mstart - first index you will use in first coordinate direction (often 0)
2808: - nstart - first index in the second coordinate direction (often 0)

2810:   Output Parameter:
2811: . a - location to put pointer to the array

2813:   Level: developer

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

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

2823: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2824:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2825:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2826: @*/
2827: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2828: {
2829:   PetscInt     i, N;
2830:   PetscScalar *aa;

2832:   PetscFunctionBegin;
2834:   PetscAssertPointer(a, 6);
2836:   PetscCall(VecGetLocalSize(x, &N));
2837:   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);
2838:   PetscCall(VecGetArrayWrite(x, &aa));

2840:   PetscCall(PetscMalloc1(m, a));
2841:   for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2842:   *a -= mstart;
2843:   PetscFunctionReturn(PETSC_SUCCESS);
2844: }

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

2849:   Logically Collective

2851:   Input Parameters:
2852: + x      - the vector
2853: . m      - first dimension of two dimensional array
2854: . n      - second dimension of the two dimensional array
2855: . mstart - first index you will use in first coordinate direction (often 0)
2856: . nstart - first index in the second coordinate direction (often 0)
2857: - a      - location of pointer to array obtained from `VecGetArray2d()`

2859:   Level: developer

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

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

2869: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2870:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2871:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2872: @*/
2873: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2874: {
2875:   void *dummy;

2877:   PetscFunctionBegin;
2879:   PetscAssertPointer(a, 6);
2881:   dummy = (void *)(*a + mstart);
2882:   PetscCall(PetscFree(dummy));
2883:   PetscCall(VecRestoreArray(x, NULL));
2884:   *a = NULL;
2885:   PetscFunctionReturn(PETSC_SUCCESS);
2886: }

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

2891:   Logically Collective

2893:   Input Parameters:
2894: + x      - the vector
2895: . m      - first dimension of two dimensional array
2896: . n      - second dimension of the two dimensional array
2897: . mstart - first index you will use in first coordinate direction (often 0)
2898: . nstart - first index in the second coordinate direction (often 0)
2899: - a      - location of pointer to array obtained from `VecGetArray2d()`

2901:   Level: developer

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

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

2911: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2912:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2913:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2914: @*/
2915: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2916: {
2917:   void *dummy;

2919:   PetscFunctionBegin;
2921:   PetscAssertPointer(a, 6);
2923:   dummy = (void *)(*a + mstart);
2924:   PetscCall(PetscFree(dummy));
2925:   PetscCall(VecRestoreArrayWrite(x, NULL));
2926:   PetscFunctionReturn(PETSC_SUCCESS);
2927: }

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

2934:   Logically Collective

2936:   Input Parameters:
2937: + x      - the vector
2938: . m      - first dimension of two dimensional array
2939: - mstart - first index you will use in first coordinate direction (often 0)

2941:   Output Parameter:
2942: . a - location to put pointer to the array

2944:   Level: developer

2946:   Notes:
2947:   For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2948:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2949:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

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

2953: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2954:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2955:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2956: @*/
2957: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2958: {
2959:   PetscInt N;

2961:   PetscFunctionBegin;
2963:   PetscAssertPointer(a, 4);
2965:   PetscCall(VecGetLocalSize(x, &N));
2966:   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);
2967:   PetscCall(VecGetArray(x, a));
2968:   *a -= mstart;
2969:   PetscFunctionReturn(PETSC_SUCCESS);
2970: }

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

2977:   Logically Collective

2979:   Input Parameters:
2980: + x      - the vector
2981: . m      - first dimension of two dimensional array
2982: - mstart - first index you will use in first coordinate direction (often 0)

2984:   Output Parameter:
2985: . a - location to put pointer to the array

2987:   Level: developer

2989:   Notes:
2990:   For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2991:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2992:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

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

2996: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2997:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2998:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2999: @*/
3000: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3001: {
3002:   PetscInt N;

3004:   PetscFunctionBegin;
3006:   PetscAssertPointer(a, 4);
3008:   PetscCall(VecGetLocalSize(x, &N));
3009:   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);
3010:   PetscCall(VecGetArrayWrite(x, a));
3011:   *a -= mstart;
3012:   PetscFunctionReturn(PETSC_SUCCESS);
3013: }

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

3018:   Logically Collective

3020:   Input Parameters:
3021: + x      - the vector
3022: . m      - first dimension of two dimensional array
3023: . mstart - first index you will use in first coordinate direction (often 0)
3024: - a      - location of pointer to array obtained from `VecGetArray1d()`

3026:   Level: developer

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

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

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

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

3053:   Logically Collective

3055:   Input Parameters:
3056: + x      - the vector
3057: . m      - first dimension of two dimensional array
3058: . mstart - first index you will use in first coordinate direction (often 0)
3059: - a      - location of pointer to array obtained from `VecGetArray1d()`

3061:   Level: developer

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

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

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

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

3090:   Logically Collective

3092:   Input Parameters:
3093: + x      - the vector
3094: . m      - first dimension of three dimensional array
3095: . n      - second dimension of three dimensional array
3096: . p      - third dimension of three dimensional array
3097: . mstart - first index you will use in first coordinate direction (often 0)
3098: . nstart - first index in the second coordinate direction (often 0)
3099: - pstart - first index in the third coordinate direction (often 0)

3101:   Output Parameter:
3102: . a - location to put pointer to the array

3104:   Level: developer

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

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

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

3123:   PetscFunctionBegin;
3125:   PetscAssertPointer(a, 8);
3127:   PetscCall(VecGetLocalSize(x, &N));
3128:   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);
3129:   PetscCall(VecGetArray(x, &aa));

3131:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3132:   b = (PetscScalar **)((*a) + m);
3133:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3134:   for (i = 0; i < m; i++)
3135:     for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3136:   *a -= mstart;
3137:   PetscFunctionReturn(PETSC_SUCCESS);
3138: }

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

3145:   Logically Collective

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

3156:   Output Parameter:
3157: . a - location to put pointer to the array

3159:   Level: developer

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

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

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

3178:   PetscFunctionBegin;
3180:   PetscAssertPointer(a, 8);
3182:   PetscCall(VecGetLocalSize(x, &N));
3183:   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);
3184:   PetscCall(VecGetArrayWrite(x, &aa));

3186:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3187:   b = (PetscScalar **)((*a) + m);
3188:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3189:   for (i = 0; i < m; i++)
3190:     for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;

3192:   *a -= mstart;
3193:   PetscFunctionReturn(PETSC_SUCCESS);
3194: }

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

3199:   Logically Collective

3201:   Input Parameters:
3202: + x      - the vector
3203: . m      - first dimension of three dimensional array
3204: . n      - second dimension of the three dimensional array
3205: . p      - third dimension of the three dimensional array
3206: . mstart - first index you will use in first coordinate direction (often 0)
3207: . nstart - first index in the second coordinate direction (often 0)
3208: . pstart - first index in the third coordinate direction (often 0)
3209: - a      - location of pointer to array obtained from VecGetArray3d()

3211:   Level: developer

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

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

3221: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3222:           `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3223:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3224: @*/
3225: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3226: {
3227:   void *dummy;

3229:   PetscFunctionBegin;
3231:   PetscAssertPointer(a, 8);
3233:   dummy = (void *)(*a + mstart);
3234:   PetscCall(PetscFree(dummy));
3235:   PetscCall(VecRestoreArray(x, NULL));
3236:   *a = NULL;
3237:   PetscFunctionReturn(PETSC_SUCCESS);
3238: }

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

3243:   Logically Collective

3245:   Input Parameters:
3246: + x      - the vector
3247: . m      - first dimension of three dimensional array
3248: . n      - second dimension of the three dimensional array
3249: . p      - third dimension of the three dimensional array
3250: . mstart - first index you will use in first coordinate direction (often 0)
3251: . nstart - first index in the second coordinate direction (often 0)
3252: . pstart - first index in the third coordinate direction (often 0)
3253: - a      - location of pointer to array obtained from VecGetArray3d()

3255:   Level: developer

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

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

3265: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3266:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3267:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3268: @*/
3269: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3270: {
3271:   void *dummy;

3273:   PetscFunctionBegin;
3275:   PetscAssertPointer(a, 8);
3277:   dummy = (void *)(*a + mstart);
3278:   PetscCall(PetscFree(dummy));
3279:   PetscCall(VecRestoreArrayWrite(x, NULL));
3280:   *a = NULL;
3281:   PetscFunctionReturn(PETSC_SUCCESS);
3282: }

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

3289:   Logically Collective

3291:   Input Parameters:
3292: + x      - the vector
3293: . m      - first dimension of four dimensional array
3294: . n      - second dimension of four dimensional array
3295: . p      - third dimension of four dimensional array
3296: . q      - fourth dimension of four dimensional array
3297: . mstart - first index you will use in first coordinate direction (often 0)
3298: . nstart - first index in the second coordinate direction (often 0)
3299: . pstart - first index in the third coordinate direction (often 0)
3300: - qstart - first index in the fourth coordinate direction (often 0)

3302:   Output Parameter:
3303: . a - location to put pointer to the array

3305:   Level: developer

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

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

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

3324:   PetscFunctionBegin;
3326:   PetscAssertPointer(a, 10);
3328:   PetscCall(VecGetLocalSize(x, &N));
3329:   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);
3330:   PetscCall(VecGetArray(x, &aa));

3332:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3333:   b = (PetscScalar ***)((*a) + m);
3334:   c = (PetscScalar **)(b + m * n);
3335:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3336:   for (i = 0; i < m; i++)
3337:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3338:   for (i = 0; i < m; i++)
3339:     for (j = 0; j < n; j++)
3340:       for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3341:   *a -= mstart;
3342:   PetscFunctionReturn(PETSC_SUCCESS);
3343: }

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

3350:   Logically Collective

3352:   Input Parameters:
3353: + x      - the vector
3354: . m      - first dimension of four dimensional array
3355: . n      - second dimension of four dimensional array
3356: . p      - third dimension of four dimensional array
3357: . q      - fourth dimension of four dimensional array
3358: . mstart - first index you will use in first coordinate direction (often 0)
3359: . nstart - first index in the second coordinate direction (often 0)
3360: . pstart - first index in the third coordinate direction (often 0)
3361: - qstart - first index in the fourth coordinate direction (often 0)

3363:   Output Parameter:
3364: . a - location to put pointer to the array

3366:   Level: developer

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

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

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

3385:   PetscFunctionBegin;
3387:   PetscAssertPointer(a, 10);
3389:   PetscCall(VecGetLocalSize(x, &N));
3390:   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);
3391:   PetscCall(VecGetArrayWrite(x, &aa));

3393:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3394:   b = (PetscScalar ***)((*a) + m);
3395:   c = (PetscScalar **)(b + m * n);
3396:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3397:   for (i = 0; i < m; i++)
3398:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3399:   for (i = 0; i < m; i++)
3400:     for (j = 0; j < n; j++)
3401:       for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3402:   *a -= mstart;
3403:   PetscFunctionReturn(PETSC_SUCCESS);
3404: }

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

3409:   Logically Collective

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

3423:   Level: developer

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

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

3433: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3434:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3435:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecGet`
3436: @*/
3437: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3438: {
3439:   void *dummy;

3441:   PetscFunctionBegin;
3443:   PetscAssertPointer(a, 10);
3445:   dummy = (void *)(*a + mstart);
3446:   PetscCall(PetscFree(dummy));
3447:   PetscCall(VecRestoreArray(x, NULL));
3448:   *a = NULL;
3449:   PetscFunctionReturn(PETSC_SUCCESS);
3450: }

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

3455:   Logically Collective

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

3469:   Level: developer

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

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

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

3487:   PetscFunctionBegin;
3489:   PetscAssertPointer(a, 10);
3491:   dummy = (void *)(*a + mstart);
3492:   PetscCall(PetscFree(dummy));
3493:   PetscCall(VecRestoreArrayWrite(x, NULL));
3494:   *a = NULL;
3495:   PetscFunctionReturn(PETSC_SUCCESS);
3496: }

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

3503:   Logically Collective

3505:   Input Parameters:
3506: + x      - the vector
3507: . m      - first dimension of two dimensional array
3508: . n      - second dimension of two dimensional array
3509: . mstart - first index you will use in first coordinate direction (often 0)
3510: - nstart - first index in the second coordinate direction (often 0)

3512:   Output Parameter:
3513: . a - location to put pointer to the array

3515:   Level: developer

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

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

3525: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3526:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3527:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3528: @*/
3529: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3530: {
3531:   PetscInt           i, N;
3532:   const PetscScalar *aa;

3534:   PetscFunctionBegin;
3536:   PetscAssertPointer(a, 6);
3538:   PetscCall(VecGetLocalSize(x, &N));
3539:   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);
3540:   PetscCall(VecGetArrayRead(x, &aa));

3542:   PetscCall(PetscMalloc1(m, a));
3543:   for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3544:   *a -= mstart;
3545:   PetscFunctionReturn(PETSC_SUCCESS);
3546: }

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

3551:   Logically Collective

3553:   Input Parameters:
3554: + x      - the vector
3555: . m      - first dimension of two dimensional array
3556: . n      - second dimension of the two dimensional array
3557: . mstart - first index you will use in first coordinate direction (often 0)
3558: . nstart - first index in the second coordinate direction (often 0)
3559: - a      - location of pointer to array obtained from VecGetArray2d()

3561:   Level: developer

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

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

3571: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3572:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3573:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3574: @*/
3575: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3576: {
3577:   void *dummy;

3579:   PetscFunctionBegin;
3581:   PetscAssertPointer(a, 6);
3583:   dummy = (void *)(*a + mstart);
3584:   PetscCall(PetscFree(dummy));
3585:   PetscCall(VecRestoreArrayRead(x, NULL));
3586:   *a = NULL;
3587:   PetscFunctionReturn(PETSC_SUCCESS);
3588: }

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

3595:   Logically Collective

3597:   Input Parameters:
3598: + x      - the vector
3599: . m      - first dimension of two dimensional array
3600: - mstart - first index you will use in first coordinate direction (often 0)

3602:   Output Parameter:
3603: . a - location to put pointer to the array

3605:   Level: developer

3607:   Notes:
3608:   For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3609:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3610:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

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

3614: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3615:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3616:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3617: @*/
3618: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3619: {
3620:   PetscInt N;

3622:   PetscFunctionBegin;
3624:   PetscAssertPointer(a, 4);
3626:   PetscCall(VecGetLocalSize(x, &N));
3627:   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);
3628:   PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3629:   *a -= mstart;
3630:   PetscFunctionReturn(PETSC_SUCCESS);
3631: }

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

3636:   Logically Collective

3638:   Input Parameters:
3639: + x      - the vector
3640: . m      - first dimension of two dimensional array
3641: . mstart - first index you will use in first coordinate direction (often 0)
3642: - a      - location of pointer to array obtained from `VecGetArray1dRead()`

3644:   Level: developer

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

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

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

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

3673:   Logically Collective

3675:   Input Parameters:
3676: + x      - the vector
3677: . m      - first dimension of three dimensional array
3678: . n      - second dimension of three dimensional array
3679: . p      - third dimension of three dimensional array
3680: . mstart - first index you will use in first coordinate direction (often 0)
3681: . nstart - first index in the second coordinate direction (often 0)
3682: - pstart - first index in the third coordinate direction (often 0)

3684:   Output Parameter:
3685: . a - location to put pointer to the array

3687:   Level: developer

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

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

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

3707:   PetscFunctionBegin;
3709:   PetscAssertPointer(a, 8);
3711:   PetscCall(VecGetLocalSize(x, &N));
3712:   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);
3713:   PetscCall(VecGetArrayRead(x, &aa));

3715:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3716:   b = (PetscScalar **)((*a) + m);
3717:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3718:   for (i = 0; i < m; i++)
3719:     for (j = 0; j < n; j++) b[i * n + j] = (PetscScalar *)aa + i * n * p + j * p - pstart;
3720:   *a -= mstart;
3721:   PetscFunctionReturn(PETSC_SUCCESS);
3722: }

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

3727:   Logically Collective

3729:   Input Parameters:
3730: + x      - the vector
3731: . m      - first dimension of three dimensional array
3732: . n      - second dimension of the three dimensional array
3733: . p      - third dimension of the three dimensional array
3734: . mstart - first index you will use in first coordinate direction (often 0)
3735: . nstart - first index in the second coordinate direction (often 0)
3736: . pstart - first index in the third coordinate direction (often 0)
3737: - a      - location of pointer to array obtained from `VecGetArray3dRead()`

3739:   Level: developer

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

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

3749: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3750:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3751:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3752: @*/
3753: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3754: {
3755:   void *dummy;

3757:   PetscFunctionBegin;
3759:   PetscAssertPointer(a, 8);
3761:   dummy = (void *)(*a + mstart);
3762:   PetscCall(PetscFree(dummy));
3763:   PetscCall(VecRestoreArrayRead(x, NULL));
3764:   *a = NULL;
3765:   PetscFunctionReturn(PETSC_SUCCESS);
3766: }

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

3773:   Logically Collective

3775:   Input Parameters:
3776: + x      - the vector
3777: . m      - first dimension of four dimensional array
3778: . n      - second dimension of four dimensional array
3779: . p      - third dimension of four dimensional array
3780: . q      - fourth dimension of four dimensional array
3781: . mstart - first index you will use in first coordinate direction (often 0)
3782: . nstart - first index in the second coordinate direction (often 0)
3783: . pstart - first index in the third coordinate direction (often 0)
3784: - qstart - first index in the fourth coordinate direction (often 0)

3786:   Output Parameter:
3787: . a - location to put pointer to the array

3789:   Level: beginner

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

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

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

3809:   PetscFunctionBegin;
3811:   PetscAssertPointer(a, 10);
3813:   PetscCall(VecGetLocalSize(x, &N));
3814:   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);
3815:   PetscCall(VecGetArrayRead(x, &aa));

3817:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3818:   b = (PetscScalar ***)((*a) + m);
3819:   c = (PetscScalar **)(b + m * n);
3820:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3821:   for (i = 0; i < m; i++)
3822:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3823:   for (i = 0; i < m; i++)
3824:     for (j = 0; j < n; j++)
3825:       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;
3826:   *a -= mstart;
3827:   PetscFunctionReturn(PETSC_SUCCESS);
3828: }

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

3833:   Logically Collective

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

3847:   Level: beginner

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

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

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

3865:   PetscFunctionBegin;
3867:   PetscAssertPointer(a, 10);
3869:   dummy = (void *)(*a + mstart);
3870:   PetscCall(PetscFree(dummy));
3871:   PetscCall(VecRestoreArrayRead(x, NULL));
3872:   *a = NULL;
3873:   PetscFunctionReturn(PETSC_SUCCESS);
3874: }

3876: #if defined(PETSC_USE_DEBUG)

3878: /*@
3879:   VecLockGet  - Gets the current lock status of a vector

3881:   Logically Collective

3883:   Input Parameter:
3884: . x - the vector

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

3890:   Level: advanced

3892: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3893: @*/
3894: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3895: {
3896:   PetscFunctionBegin;
3898:   *state = x->lock;
3899:   PetscFunctionReturn(PETSC_SUCCESS);
3900: }

3902: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3903: {
3904:   PetscFunctionBegin;
3906:   PetscAssertPointer(file, 2);
3907:   PetscAssertPointer(func, 3);
3908:   PetscAssertPointer(line, 4);
3909:   #if !PetscDefined(HAVE_THREADSAFETY)
3910:   {
3911:     const int index = x->lockstack.currentsize - 1;

3913:     PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupted vec lock stack, have negative index %d", index);
3914:     *file = x->lockstack.file[index];
3915:     *func = x->lockstack.function[index];
3916:     *line = x->lockstack.line[index];
3917:   }
3918:   #else
3919:   *file = NULL;
3920:   *func = NULL;
3921:   *line = 0;
3922:   #endif
3923:   PetscFunctionReturn(PETSC_SUCCESS);
3924: }

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

3929:   Logically Collective

3931:   Input Parameter:
3932: . x - the vector

3934:   Level: intermediate

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

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

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

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

3972: /*@
3973:   VecLockReadPop  - Pops a read-only lock from a vector

3975:   Logically Collective

3977:   Input Parameter:
3978: . x - the vector

3980:   Level: intermediate

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

3993:     PetscStackPop_Private(x->lockstack, previous);
3994:   }
3995:   #endif
3996:   PetscFunctionReturn(PETSC_SUCCESS);
3997: }

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

4002:   Logically Collective

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

4008:   Level: intermediate

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

4016: .vb
4017:        VecGetArray(x,&xdata); // begin phase
4018:        VecLockWriteSet(v,PETSC_TRUE);

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

4022:        VecRestoreArray(x,&vdata); // end phase
4023:        VecLockWriteSet(v,PETSC_FALSE);
4024: .ve

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

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

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

4050:   Level: deprecated

4052: .seealso: [](ch_vectors), `Vec`, `VecLockReadPush()`
4053: @*/
4054: PetscErrorCode VecLockPush(Vec x)
4055: {
4056:   PetscFunctionBegin;
4057:   PetscCall(VecLockReadPush(x));
4058:   PetscFunctionReturn(PETSC_SUCCESS);
4059: }

4061: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
4062: /*@
4063:   VecLockPop  - Pops a read-only lock from a vector

4065:   Level: deprecated

4067: .seealso: [](ch_vectors), `Vec`, `VecLockReadPop()`
4068: @*/
4069: PetscErrorCode VecLockPop(Vec x)
4070: {
4071:   PetscFunctionBegin;
4072:   PetscCall(VecLockReadPop(x));
4073:   PetscFunctionReturn(PETSC_SUCCESS);
4074: }

4076: #endif