Actual source code: rvector.c

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

  7: PetscInt VecGetSubVectorSavedStateId = -1;

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

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

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

 39:   Logically Collective

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

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

 48:   Level: advanced

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

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

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

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

 80:   Collective

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

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

 89:   Level: intermediate

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

 98:   Use `VecTDot()` for the indefinite form
 99: $     val = (x,y) = y^T x,
100:   where y^T denotes the transpose of y.

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

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

127: /*@
128:   VecDotRealPart - Computes the real part of the vector dot product.

130:   Collective

132:   Input Parameters:
133: + x - first vector
134: - y - second vector

136:   Output Parameter:
137: . val - the real part of the dot product;

139:   Level: intermediate

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

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

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

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

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

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

164: /*@
165:   VecNorm  - Computes the vector norm.

167:   Collective

169:   Input Parameters:
170: + x    - the vector
171: - type - the type of the norm requested

173:   Output Parameter:
174: . val - the norm

176:   Level: intermediate

178:   Notes:
179:   See `NormType` for descriptions of each norm.

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

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

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

199:   PetscFunctionBegin;
202:   VecCheckAssembled(x);
204:   PetscAssertPointer(val, 3);

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

224:   PetscCall(VecLockReadPush(x));
225:   PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
226:   PetscUseTypeMethod(x, norm, type, val);
227:   PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
228:   PetscCall(VecLockReadPop(x));

230:   if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
231:   PetscFunctionReturn(PETSC_SUCCESS);
232: }

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

237:   Not Collective

239:   Input Parameters:
240: + x    - the vector
241: - 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
242:           `NORM_1_AND_2`, which computes both norms and stores them
243:           in a two element array.

245:   Output Parameters:
246: + available - `PETSC_TRUE` if the val returned is valid
247: - val       - the norm

249:   Level: intermediate

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

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

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

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

278:   Collective

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

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

286:   Level: intermediate

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

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

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

315:   Collective

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

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

324:   Level: intermediate

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

329:   Returns the smallest index with the maximum value

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

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

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

355:   Collective

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

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

364:   Level: intermediate

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

369:   This returns the smallest index with the minimum value

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

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

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

396:   Collective

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

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

405:   Level: intermediate

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

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

416: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
417: @*/
418: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
419: {
420:   PetscFunctionBegin;
423:   PetscAssertPointer(val, 3);
426:   PetscCheckSameTypeAndComm(x, 1, y, 2);
427:   VecCheckSameSize(x, 1, y, 2);
428:   VecCheckAssembled(x);
429:   VecCheckAssembled(y);

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

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

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

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

458:   PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
459:   VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
460:   PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));

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

471: /*@
472:   VecScale - Scales a vector.

474:   Logically Collective

476:   Input Parameters:
477: + x     - the vector
478: - alpha - the scalar

480:   Level: intermediate

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

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

494: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
495: {
496:   PetscFunctionBegin;
499:   VecCheckAssembled(x);
501:   PetscCall(VecSetErrorIfLocked(x, 1));

503:   if (alpha == 0) {
504:     PetscReal norm;
505:     PetscBool set;

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

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

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

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

541:   Logically Collective

543:   Input Parameters:
544: + x     - the vector
545: - alpha - the scalar

547:   Level: beginner

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

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

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

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

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

598:   Logically Collective

600:   Input Parameters:
601: + alpha - the scalar
602: . x     - vector scale by `alpha`
603: - y     - vector accumulated into

605:   Output Parameter:
606: . y - output vector

608:   Level: intermediate

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

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

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

660: /*@
661:   VecAYPX - Computes `y = x + beta y`.

663:   Logically Collective

665:   Input Parameters:
666: + beta - the scalar
667: . x    - the unscaled vector
668: - y    - the vector to be scaled

670:   Output Parameter:
671: . y - output vector

673:   Level: intermediate

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

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

687: PetscErrorCode VecAXPBYAsync_Private(Vec y, PetscScalar alpha, PetscScalar beta, Vec x, PetscDeviceContext dctx)
688: {
689:   PetscFunctionBegin;
694:   PetscCheckSameTypeAndComm(x, 4, y, 1);
695:   VecCheckSameSize(y, 1, x, 4);
696:   VecCheckAssembled(x);
697:   VecCheckAssembled(y);
700:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
701:   if (x == y) {
702:     PetscCall(VecScale(y, alpha + beta));
703:     PetscFunctionReturn(PETSC_SUCCESS);
704:   }

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

716: /*@
717:   VecAXPBY - Computes `y = alpha x + beta y`.

719:   Logically Collective

721:   Input Parameters:
722: + alpha - first scalar
723: . beta  - second scalar
724: . x     - the first scaled vector
725: - y     - the second scaled vector

727:   Output Parameter:
728: . y - output vector

730:   Level: intermediate

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

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

744: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
745: {
746:   PetscFunctionBegin;
753:   PetscCheckSameTypeAndComm(x, 5, y, 6);
754:   PetscCheckSameTypeAndComm(x, 5, z, 1);
755:   VecCheckSameSize(x, 5, y, 6);
756:   VecCheckSameSize(x, 5, z, 1);
757:   PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
758:   PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
759:   VecCheckAssembled(x);
760:   VecCheckAssembled(y);
761:   VecCheckAssembled(z);
765:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);

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

781:   Logically Collective

783:   Input Parameters:
784: + alpha - first scalar
785: . beta  - second scalar
786: . gamma - third scalar
787: . x     - first vector
788: . y     - second vector
789: - z     - third vector

791:   Output Parameter:
792: . z - output vector

794:   Level: intermediate

796:   Note:
797:   `x`, `y` and `z` must be different vectors

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

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

811: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
812: {
813:   PetscFunctionBegin;
820:   PetscCheckSameTypeAndComm(x, 3, y, 4);
821:   PetscCheckSameTypeAndComm(y, 4, w, 1);
822:   VecCheckSameSize(x, 3, y, 4);
823:   VecCheckSameSize(x, 3, w, 1);
824:   PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
825:   PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
826:   VecCheckAssembled(x);
827:   VecCheckAssembled(y);
829:   PetscCall(VecSetErrorIfLocked(w, 1));

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

846: /*@
847:   VecWAXPY - Computes `w = alpha x + y`.

849:   Logically Collective

851:   Input Parameters:
852: + alpha - the scalar
853: . x     - first vector, multiplied by `alpha`
854: - y     - second vector

856:   Output Parameter:
857: . w - the result

859:   Level: intermediate

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

864:   Developer Notes:
865:   The implementation is optimized for alpha of -1.0, 0.0, and 1.0

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

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

879:   Not Collective

881:   Input Parameters:
882: + x    - vector to insert in
883: . ni   - number of elements to add
884: . ix   - indices where to add
885: . y    - array of values
886: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries

888:   Level: beginner

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

895:   Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
896:   options cannot be mixed without intervening calls to the assembly
897:   routines.

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

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

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

910:   Fortran Note:
911:   If any of `ix` and `y` are scalars pass them using, for example,
912: .vb
913:   VecSetValues(mat, one, [ix], [y], INSERT_VALUES)
914: .ve

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

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

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

939:   Not Collective

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

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

949:   Level: beginner

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

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

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

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

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

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

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

982:   Not Collective

984:   Input Parameters:
985: + x    - vector to insert in
986: . ni   - number of blocks to add
987: . ix   - indices where to add in block count, rather than element count
988: . y    - array of values
989: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries

991:   Level: intermediate

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

997:   Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
998:   options cannot be mixed without intervening calls to the assembly
999:   routines.

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

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

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

1011:   Fortran Note:
1012:   If any of `ix` and `y` are scalars pass them using, for example,
1013: .vb
1014:   VecSetValuesBlocked(mat, one, [ix], [y], INSERT_VALUES)
1015: .ve

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

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

1036: /*@
1037:   VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1038:   using a local ordering of the nodes.

1040:   Not Collective

1042:   Input Parameters:
1043: + x    - vector to insert in
1044: . ni   - number of elements to add
1045: . ix   - indices where to add
1046: . y    - array of values
1047: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries

1049:   Level: intermediate

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

1054:   Calls to `VecSetValuesLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1055:   options cannot be mixed without intervening calls to the assembly
1056:   routines.

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

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

1063:   Fortran Note:
1064:   If any of `ix` and `y` are scalars pass them using, for example,
1065: .vb
1066:   VecSetValuesLocal(mat, one, [ix], [y], INSERT_VALUES)
1067: .ve

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

1076:   PetscFunctionBeginHot;
1078:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1079:   PetscAssertPointer(ix, 3);
1080:   PetscAssertPointer(y, 4);

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

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

1102:   Not Collective

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

1111:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

1193:   Collective

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

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

1203:   Level: intermediate

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

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

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

1224: /*@
1225:   VecMDot - Computes multiple vector dot products.

1227:   Collective

1229:   Input Parameters:
1230: + x  - one vector
1231: . nv - number of vectors
1232: - y  - array of vectors.

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

1237:   Level: intermediate

1239:   Notes for Users of Complex Numbers:
1240:   For complex vectors, `VecMDot()` computes
1241: $     val = (x,y) = y^H x,
1242:   where y^H denotes the conjugate transpose of y.

1244:   Use `VecMTDot()` for the indefinite form
1245: $     val = (x,y) = y^T x,
1246:   where y^T denotes the transpose of y.

1248: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`
1249: @*/
1250: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1251: {
1252:   PetscFunctionBegin;
1254:   PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1255:   PetscFunctionReturn(PETSC_SUCCESS);
1256: }

1258: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
1259: {
1260:   PetscFunctionBegin;
1262:   VecCheckAssembled(y);
1264:   PetscCall(VecSetErrorIfLocked(y, 1));
1265:   PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1266:   if (nv) {
1267:     PetscInt zeros = 0;

1269:     PetscAssertPointer(alpha, 3);
1270:     PetscAssertPointer(x, 4);
1271:     for (PetscInt i = 0; i < nv; ++i) {
1275:       PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1276:       VecCheckSameSize(y, 1, x[i], 4);
1277:       PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1278:       VecCheckAssembled(x[i]);
1279:       PetscCall(VecLockReadPush(x[i]));
1280:       zeros += alpha[i] == (PetscScalar)0.0;
1281:     }

1283:     if (zeros < nv) {
1284:       PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1285:       VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1286:       PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1287:       PetscCall(PetscObjectStateIncrease((PetscObject)y));
1288:     }

1290:     for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1291:   }
1292:   PetscFunctionReturn(PETSC_SUCCESS);
1293: }

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

1298:   Logically Collective

1300:   Input Parameters:
1301: + nv    - number of scalars and x-vectors
1302: . alpha - array of scalars
1303: . y     - one vector
1304: - x     - array of vectors

1306:   Level: intermediate

1308:   Note:
1309:   `y` cannot be any of the `x` vectors

1311: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`,`VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1312: @*/
1313: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1314: {
1315:   PetscFunctionBegin;
1316:   PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1317:   PetscFunctionReturn(PETSC_SUCCESS);
1318: }

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

1323:   Logically Collective

1325:   Input Parameters:
1326: + nv    - number of scalars and x-vectors
1327: . alpha - array of scalars
1328: . beta  - scalar
1329: . y     - one vector
1330: - x     - array of vectors

1332:   Level: intermediate

1334:   Note:
1335:   `y` cannot be any of the `x` vectors.

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

1340: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1341: @*/
1342: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
1343: {
1344:   PetscFunctionBegin;
1346:   VecCheckAssembled(y);
1348:   PetscCall(VecSetErrorIfLocked(y, 1));
1349:   PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);

1352:   if (y->ops->maxpby) {
1353:     PetscInt zeros = 0;

1355:     if (nv) {
1356:       PetscAssertPointer(alpha, 3);
1357:       PetscAssertPointer(x, 5);
1358:     }

1360:     for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
1364:       PetscCheckSameTypeAndComm(y, 1, x[i], 5);
1365:       VecCheckSameSize(y, 1, x[i], 5);
1366:       PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1367:       VecCheckAssembled(x[i]);
1368:       PetscCall(VecLockReadPush(x[i]));
1369:       zeros += alpha[i] == (PetscScalar)0.0;
1370:     }

1372:     if (zeros < nv) { // has nonzero alpha
1373:       PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1374:       PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1375:       PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1376:       PetscCall(PetscObjectStateIncrease((PetscObject)y));
1377:     } else {
1378:       PetscCall(VecScale(y, beta));
1379:     }

1381:     for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1382:   } else { // no maxpby
1383:     if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1384:     else PetscCall(VecScale(y, beta));
1385:     PetscCall(VecMAXPY(y, nv, alpha, x));
1386:   }
1387:   PetscFunctionReturn(PETSC_SUCCESS);
1388: }

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

1395:   Collective

1397:   Input Parameters:
1398: + nx - number of vectors to be concatenated
1399: - X  - array containing the vectors to be concatenated in the order of concatenation

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

1405:   Level: advanced

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

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

1417: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1418: @*/
1419: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1420: {
1421:   MPI_Comm comm;
1422:   VecType  vec_type;
1423:   Vec      Ytmp, Xtmp;
1424:   IS      *is_tmp;
1425:   PetscInt i, shift = 0, Xnl, Xng, Xbegin;

1427:   PetscFunctionBegin;
1431:   PetscAssertPointer(Y, 3);

1433:   if ((*X)->ops->concatenate) {
1434:     /* use the dedicated concatenation function if available */
1435:     PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1436:   } else {
1437:     /* loop over vectors and start creating IS */
1438:     comm = PetscObjectComm((PetscObject)*X);
1439:     PetscCall(VecGetType(*X, &vec_type));
1440:     PetscCall(PetscMalloc1(nx, &is_tmp));
1441:     for (i = 0; i < nx; i++) {
1442:       PetscCall(VecGetSize(X[i], &Xng));
1443:       PetscCall(VecGetLocalSize(X[i], &Xnl));
1444:       PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1445:       PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1446:       shift += Xng;
1447:     }
1448:     /* create the concatenated vector */
1449:     PetscCall(VecCreate(comm, &Ytmp));
1450:     PetscCall(VecSetType(Ytmp, vec_type));
1451:     PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1452:     PetscCall(VecSetUp(Ytmp));
1453:     /* copy data from X array to Y and return */
1454:     for (i = 0; i < nx; i++) {
1455:       PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1456:       PetscCall(VecCopy(X[i], Xtmp));
1457:       PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1458:     }
1459:     *Y = Ytmp;
1460:     if (x_is) {
1461:       *x_is = is_tmp;
1462:     } else {
1463:       for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1464:       PetscCall(PetscFree(is_tmp));
1465:     }
1466:   }
1467:   PetscFunctionReturn(PETSC_SUCCESS);
1468: }

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

1473:     Input Parameters:
1474: +   X - the original vector
1475: -   is - the index set of the subvector

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

1482: */
1483: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1484: {
1485:   PetscInt  gstart, gend, lstart;
1486:   PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1487:   PetscInt  n, N, ibs, vbs, bs = -1;

1489:   PetscFunctionBegin;
1490:   PetscCall(ISGetLocalSize(is, &n));
1491:   PetscCall(ISGetSize(is, &N));
1492:   PetscCall(ISGetBlockSize(is, &ibs));
1493:   PetscCall(VecGetBlockSize(X, &vbs));
1494:   PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1495:   PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1496:   /* block size is given by IS if ibs > 1; otherwise, check the vector */
1497:   if (ibs > 1) {
1498:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1499:     bs = ibs;
1500:   } else {
1501:     if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1502:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1503:     if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1504:   }

1506:   *contig    = red[0];
1507:   *start     = lstart;
1508:   *blocksize = bs;
1509:   PetscFunctionReturn(PETSC_SUCCESS);
1510: }

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

1514:     Input Parameters:
1515: +   X - the original vector
1516: .   is - the index set of the subvector
1517: -   bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()

1519:     Output Parameter:
1520: .   Z  - the subvector, which will compose the VecScatter context on output
1521: */
1522: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1523: {
1524:   PetscInt   n, N;
1525:   VecScatter vscat;
1526:   Vec        Y;

1528:   PetscFunctionBegin;
1529:   PetscCall(ISGetLocalSize(is, &n));
1530:   PetscCall(ISGetSize(is, &N));
1531:   PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1532:   PetscCall(VecSetSizes(Y, n, N));
1533:   PetscCall(VecSetBlockSize(Y, bs));
1534:   PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1535:   PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1536:   PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1537:   PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1538:   PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1539:   PetscCall(VecScatterDestroy(&vscat));
1540:   *Z = Y;
1541:   PetscFunctionReturn(PETSC_SUCCESS);
1542: }

1544: /*@
1545:   VecGetSubVector - Gets a vector representing part of another vector

1547:   Collective

1549:   Input Parameters:
1550: + X  - vector from which to extract a subvector
1551: - is - index set representing portion of `X` to extract

1553:   Output Parameter:
1554: . Y - subvector corresponding to `is`

1556:   Level: advanced

1558:   Notes:
1559:   The subvector `Y` should be returned with `VecRestoreSubVector()`.
1560:   `X` and `is` must be defined on the same communicator

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

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

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

1569: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1570: @*/
1571: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1572: {
1573:   Vec Z;

1575:   PetscFunctionBegin;
1578:   PetscCheckSameComm(X, 1, is, 2);
1579:   PetscAssertPointer(Y, 3);
1580:   if (X->ops->getsubvector) {
1581:     PetscUseTypeMethod(X, getsubvector, is, &Z);
1582:   } else { /* Default implementation currently does no caching */
1583:     PetscBool contig;
1584:     PetscInt  n, N, start, bs;

1586:     PetscCall(ISGetLocalSize(is, &n));
1587:     PetscCall(ISGetSize(is, &N));
1588:     PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1589:     if (contig) { /* We can do a no-copy implementation */
1590:       const PetscScalar *x;
1591:       PetscInt           state = 0;
1592:       PetscBool          isstd, iscuda, iship;

1594:       PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1595:       PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1596:       PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1597:       if (iscuda) {
1598: #if defined(PETSC_HAVE_CUDA)
1599:         const PetscScalar *x_d;
1600:         PetscMPIInt        size;
1601:         PetscOffloadMask   flg;

1603:         PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1604:         PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1605:         PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1606:         if (x) x += start;
1607:         if (x_d) x_d += start;
1608:         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1609:         if (size == 1) {
1610:           PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1611:         } else {
1612:           PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1613:         }
1614:         Z->offloadmask = flg;
1615: #endif
1616:       } else if (iship) {
1617: #if defined(PETSC_HAVE_HIP)
1618:         const PetscScalar *x_d;
1619:         PetscMPIInt        size;
1620:         PetscOffloadMask   flg;

1622:         PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1623:         PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1624:         PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1625:         if (x) x += start;
1626:         if (x_d) x_d += start;
1627:         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1628:         if (size == 1) {
1629:           PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1630:         } else {
1631:           PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1632:         }
1633:         Z->offloadmask = flg;
1634: #endif
1635:       } else if (isstd) {
1636:         PetscMPIInt size;

1638:         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1639:         PetscCall(VecGetArrayRead(X, &x));
1640:         if (x) x += start;
1641:         if (size == 1) {
1642:           PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1643:         } else {
1644:           PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1645:         }
1646:         PetscCall(VecRestoreArrayRead(X, &x));
1647:       } else { /* default implementation: use place array */
1648:         PetscCall(VecGetArrayRead(X, &x));
1649:         PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1650:         PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1651:         PetscCall(VecSetSizes(Z, n, N));
1652:         PetscCall(VecSetBlockSize(Z, bs));
1653:         PetscCall(VecPlaceArray(Z, PetscSafePointerPlusOffset(x, start)));
1654:         PetscCall(VecRestoreArrayRead(X, &x));
1655:       }

1657:       /* this is relevant only in debug mode */
1658:       PetscCall(VecLockGet(X, &state));
1659:       if (state) PetscCall(VecLockReadPush(Z));
1660:       Z->ops->placearray   = NULL;
1661:       Z->ops->replacearray = NULL;
1662:     } else { /* Have to create a scatter and do a copy */
1663:       PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1664:     }
1665:   }
1666:   /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1667:   if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1668:   PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1669:   *Y = Z;
1670:   PetscFunctionReturn(PETSC_SUCCESS);
1671: }

1673: /*@
1674:   VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`

1676:   Collective

1678:   Input Parameters:
1679: + X  - vector from which subvector was obtained
1680: . is - index set representing the subset of `X`
1681: - Y  - subvector being restored

1683:   Level: advanced

1685: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1686: @*/
1687: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1688: {
1689:   PETSC_UNUSED PetscObjectState dummystate = 0;
1690:   PetscBool                     unchanged;

1692:   PetscFunctionBegin;
1695:   PetscCheckSameComm(X, 1, is, 2);
1696:   PetscAssertPointer(Y, 3);

1699:   if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1700:   else {
1701:     PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1702:     if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1703:       VecScatter scatter;
1704:       PetscInt   state;

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

1709:       PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1710:       if (scatter) {
1711:         PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1712:         PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1713:       } else {
1714:         PetscBool iscuda, iship;
1715:         PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1716:         PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));

1718:         if (iscuda) {
1719: #if defined(PETSC_HAVE_CUDA)
1720:           PetscOffloadMask ymask = (*Y)->offloadmask;

1722:           /* The offloadmask of X dictates where to move memory
1723:               If X GPU data is valid, then move Y data on GPU if needed
1724:               Otherwise, move back to the CPU */
1725:           switch (X->offloadmask) {
1726:           case PETSC_OFFLOAD_BOTH:
1727:             if (ymask == PETSC_OFFLOAD_CPU) {
1728:               PetscCall(VecCUDAResetArray(*Y));
1729:             } else if (ymask == PETSC_OFFLOAD_GPU) {
1730:               X->offloadmask = PETSC_OFFLOAD_GPU;
1731:             }
1732:             break;
1733:           case PETSC_OFFLOAD_GPU:
1734:             if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1735:             break;
1736:           case PETSC_OFFLOAD_CPU:
1737:             if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1738:             break;
1739:           case PETSC_OFFLOAD_UNALLOCATED:
1740:           case PETSC_OFFLOAD_KOKKOS:
1741:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1742:           }
1743: #endif
1744:         } else if (iship) {
1745: #if defined(PETSC_HAVE_HIP)
1746:           PetscOffloadMask ymask = (*Y)->offloadmask;

1748:           /* The offloadmask of X dictates where to move memory
1749:               If X GPU data is valid, then move Y data on GPU if needed
1750:               Otherwise, move back to the CPU */
1751:           switch (X->offloadmask) {
1752:           case PETSC_OFFLOAD_BOTH:
1753:             if (ymask == PETSC_OFFLOAD_CPU) {
1754:               PetscCall(VecHIPResetArray(*Y));
1755:             } else if (ymask == PETSC_OFFLOAD_GPU) {
1756:               X->offloadmask = PETSC_OFFLOAD_GPU;
1757:             }
1758:             break;
1759:           case PETSC_OFFLOAD_GPU:
1760:             if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1761:             break;
1762:           case PETSC_OFFLOAD_CPU:
1763:             if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1764:             break;
1765:           case PETSC_OFFLOAD_UNALLOCATED:
1766:           case PETSC_OFFLOAD_KOKKOS:
1767:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1768:           }
1769: #endif
1770:         } else {
1771:           /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1772:           PetscCall(VecResetArray(*Y));
1773:         }
1774:         PetscCall(PetscObjectStateIncrease((PetscObject)X));
1775:       }
1776:     }
1777:   }
1778:   PetscCall(VecDestroy(Y));
1779:   PetscFunctionReturn(PETSC_SUCCESS);
1780: }

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

1786:   Not Collective.

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

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

1794:   Level: beginner

1796: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1797: @*/
1798: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1799: {
1800:   VecType  roottype;
1801:   PetscInt n;

1803:   PetscFunctionBegin;
1805:   PetscAssertPointer(w, 2);
1806:   if (v->ops->createlocalvector) {
1807:     PetscUseTypeMethod(v, createlocalvector, w);
1808:     PetscFunctionReturn(PETSC_SUCCESS);
1809:   }
1810:   PetscCall(VecGetRootType_Private(v, &roottype));
1811:   PetscCall(VecCreate(PETSC_COMM_SELF, w));
1812:   PetscCall(VecGetLocalSize(v, &n));
1813:   PetscCall(VecSetSizes(*w, n, n));
1814:   PetscCall(VecGetBlockSize(v, &n));
1815:   PetscCall(VecSetBlockSize(*w, n));
1816:   PetscCall(VecSetType(*w, roottype));
1817:   PetscFunctionReturn(PETSC_SUCCESS);
1818: }

1820: /*@
1821:   VecGetLocalVectorRead - Maps the local portion of a vector into a
1822:   vector.

1824:   Not Collective.

1826:   Input Parameter:
1827: . v - The vector for which the local vector is desired.

1829:   Output Parameter:
1830: . w - Upon exit this contains the local vector.

1832:   Level: beginner

1834:   Notes:
1835:   You must call `VecRestoreLocalVectorRead()` when the local
1836:   vector is no longer needed.

1838:   This function is similar to `VecGetArrayRead()` which maps the local
1839:   portion into a raw pointer.  `VecGetLocalVectorRead()` is usually
1840:   almost as efficient as `VecGetArrayRead()` but in certain circumstances
1841:   `VecGetLocalVectorRead()` can be much more efficient than
1842:   `VecGetArrayRead()`.  This is because the construction of a contiguous
1843:   array representing the vector data required by `VecGetArrayRead()` can
1844:   be an expensive operation for certain vector types.  For example, for
1845:   GPU vectors `VecGetArrayRead()` requires that the data between device
1846:   and host is synchronized.

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

1851: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1852: @*/
1853: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1854: {
1855:   PetscFunctionBegin;
1858:   VecCheckSameLocalSize(v, 1, w, 2);
1859:   if (v->ops->getlocalvectorread) {
1860:     PetscUseTypeMethod(v, getlocalvectorread, w);
1861:   } else {
1862:     PetscScalar *a;

1864:     PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1865:     PetscCall(VecPlaceArray(w, a));
1866:   }
1867:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
1868:   PetscCall(VecLockReadPush(v));
1869:   PetscCall(VecLockReadPush(w));
1870:   PetscFunctionReturn(PETSC_SUCCESS);
1871: }

1873: /*@
1874:   VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1875:   previously mapped into a vector using `VecGetLocalVectorRead()`.

1877:   Not Collective.

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

1883:   Level: beginner

1885: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1886: @*/
1887: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1888: {
1889:   PetscFunctionBegin;
1892:   if (v->ops->restorelocalvectorread) {
1893:     PetscUseTypeMethod(v, restorelocalvectorread, w);
1894:   } else {
1895:     const PetscScalar *a;

1897:     PetscCall(VecGetArrayRead(w, &a));
1898:     PetscCall(VecRestoreArrayRead(v, &a));
1899:     PetscCall(VecResetArray(w));
1900:   }
1901:   PetscCall(VecLockReadPop(v));
1902:   PetscCall(VecLockReadPop(w));
1903:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
1904:   PetscFunctionReturn(PETSC_SUCCESS);
1905: }

1907: /*@
1908:   VecGetLocalVector - Maps the local portion of a vector into a
1909:   vector.

1911:   Collective

1913:   Input Parameter:
1914: . v - The vector for which the local vector is desired.

1916:   Output Parameter:
1917: . w - Upon exit this contains the local vector.

1919:   Level: beginner

1921:   Notes:
1922:   You must call `VecRestoreLocalVector()` when the local
1923:   vector is no longer needed.

1925:   This function is similar to `VecGetArray()` which maps the local
1926:   portion into a raw pointer.  `VecGetLocalVector()` is usually about as
1927:   efficient as `VecGetArray()` but in certain circumstances
1928:   `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1929:   This is because the construction of a contiguous array representing
1930:   the vector data required by `VecGetArray()` can be an expensive
1931:   operation for certain vector types.  For example, for GPU vectors
1932:   `VecGetArray()` requires that the data between device and host is
1933:   synchronized.

1935: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1936: @*/
1937: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1938: {
1939:   PetscFunctionBegin;
1942:   VecCheckSameLocalSize(v, 1, w, 2);
1943:   if (v->ops->getlocalvector) {
1944:     PetscUseTypeMethod(v, getlocalvector, w);
1945:   } else {
1946:     PetscScalar *a;

1948:     PetscCall(VecGetArray(v, &a));
1949:     PetscCall(VecPlaceArray(w, a));
1950:   }
1951:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
1952:   PetscFunctionReturn(PETSC_SUCCESS);
1953: }

1955: /*@
1956:   VecRestoreLocalVector - Unmaps the local portion of a vector
1957:   previously mapped into a vector using `VecGetLocalVector()`.

1959:   Logically Collective.

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

1965:   Level: beginner

1967: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1968: @*/
1969: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1970: {
1971:   PetscFunctionBegin;
1974:   if (v->ops->restorelocalvector) {
1975:     PetscUseTypeMethod(v, restorelocalvector, w);
1976:   } else {
1977:     PetscScalar *a;
1978:     PetscCall(VecGetArray(w, &a));
1979:     PetscCall(VecRestoreArray(v, &a));
1980:     PetscCall(VecResetArray(w));
1981:   }
1982:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
1983:   PetscCall(PetscObjectStateIncrease((PetscObject)v));
1984:   PetscFunctionReturn(PETSC_SUCCESS);
1985: }

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

1991:   Logically Collective

1993:   Input Parameter:
1994: . x - the vector

1996:   Output Parameter:
1997: . a - location to put pointer to the array

1999:   Level: beginner

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

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

2010: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
2011:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2012: @*/
2013: PetscErrorCode VecGetArray(Vec x, PetscScalar **a)
2014: {
2015:   PetscFunctionBegin;
2017:   PetscCall(VecSetErrorIfLocked(x, 1));
2018:   if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
2019:     PetscUseTypeMethod(x, getarray, a);
2020:   } else if (x->petscnative) { /* VECSTANDARD */
2021:     *a = *((PetscScalar **)x->data);
2022:   } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
2023:   PetscFunctionReturn(PETSC_SUCCESS);
2024: }

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

2029:   Logically Collective

2031:   Input Parameters:
2032: + x - the vector
2033: - a - location of pointer to array obtained from `VecGetArray()`

2035:   Level: beginner

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

2040: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2041:           `VecGetArrayPair()`, `VecRestoreArrayPair()`
2042: @*/
2043: PetscErrorCode VecRestoreArray(Vec x, PetscScalar **a)
2044: {
2045:   PetscFunctionBegin;
2047:   if (a) PetscAssertPointer(a, 2);
2048:   if (x->ops->restorearray) {
2049:     PetscUseTypeMethod(x, restorearray, a);
2050:   } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
2051:   if (a) *a = NULL;
2052:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
2053:   PetscFunctionReturn(PETSC_SUCCESS);
2054: }
2055: /*@C
2056:   VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.

2058:   Not Collective

2060:   Input Parameter:
2061: . x - the vector

2063:   Output Parameter:
2064: . a - the array

2066:   Level: beginner

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

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

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

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

2080: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2081: @*/
2082: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar **a)
2083: {
2084:   PetscFunctionBegin;
2086:   PetscAssertPointer(a, 2);
2087:   if (x->ops->getarrayread) {
2088:     PetscUseTypeMethod(x, getarrayread, a);
2089:   } else if (x->ops->getarray) {
2090:     PetscObjectState state;

2092:     /* VECNEST, VECCUDA, VECKOKKOS etc */
2093:     // x->ops->getarray may bump the object state, but since we know this is a read-only get
2094:     // we can just undo that
2095:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2096:     PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
2097:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
2098:   } else if (x->petscnative) {
2099:     /* VECSTANDARD */
2100:     *a = *((PetscScalar **)x->data);
2101:   } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
2102:   PetscFunctionReturn(PETSC_SUCCESS);
2103: }

2105: /*@C
2106:   VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`

2108:   Not Collective

2110:   Input Parameters:
2111: + x - the vector
2112: - a - the array

2114:   Level: beginner

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

2119: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2120: @*/
2121: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar **a)
2122: {
2123:   PetscFunctionBegin;
2125:   if (a) PetscAssertPointer(a, 2);
2126:   if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
2127:     /* nothing */
2128:   } else if (x->ops->restorearrayread) { /* VECNEST */
2129:     PetscUseTypeMethod(x, restorearrayread, a);
2130:   } else { /* No one? */
2131:     PetscObjectState state;

2133:     // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2134:     // we can just undo that
2135:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2136:     PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2137:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
2138:   }
2139:   if (a) *a = NULL;
2140:   PetscFunctionReturn(PETSC_SUCCESS);
2141: }

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

2147:   Logically Collective

2149:   Input Parameter:
2150: . x - the vector

2152:   Output Parameter:
2153: . a - location to put pointer to the array

2155:   Level: intermediate

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

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

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

2166:   Fortran Notes:
2167:   `VecGetArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayF90()`

2169: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
2170:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
2171: @*/
2172: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar **a)
2173: {
2174:   PetscFunctionBegin;
2176:   PetscAssertPointer(a, 2);
2177:   PetscCall(VecSetErrorIfLocked(x, 1));
2178:   if (x->ops->getarraywrite) {
2179:     PetscUseTypeMethod(x, getarraywrite, a);
2180:   } else {
2181:     PetscCall(VecGetArray(x, a));
2182:   }
2183:   PetscFunctionReturn(PETSC_SUCCESS);
2184: }

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

2189:   Logically Collective

2191:   Input Parameters:
2192: + x - the vector
2193: - a - location of pointer to array obtained from `VecGetArray()`

2195:   Level: beginner

2197:   Fortran Notes:
2198:   `VecRestoreArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayF90()`

2200: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2201:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2202: @*/
2203: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar **a)
2204: {
2205:   PetscFunctionBegin;
2207:   if (a) PetscAssertPointer(a, 2);
2208:   if (x->ops->restorearraywrite) {
2209:     PetscUseTypeMethod(x, restorearraywrite, a);
2210:   } else if (x->ops->restorearray) {
2211:     PetscUseTypeMethod(x, restorearray, a);
2212:   }
2213:   if (a) *a = NULL;
2214:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
2215:   PetscFunctionReturn(PETSC_SUCCESS);
2216: }

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

2222:   Logically Collective; No Fortran Support

2224:   Input Parameters:
2225: + x - the vectors
2226: - n - the number of vectors

2228:   Output Parameter:
2229: . a - location to put pointer to the array

2231:   Level: intermediate

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

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

2243:   PetscFunctionBegin;
2244:   PetscAssertPointer(x, 1);
2246:   PetscAssertPointer(a, 3);
2247:   PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2248:   PetscCall(PetscMalloc1(n, &q));
2249:   for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2250:   *a = q;
2251:   PetscFunctionReturn(PETSC_SUCCESS);
2252: }

2254: /*@C
2255:   VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2256:   has been called.

2258:   Logically Collective; No Fortran Support

2260:   Input Parameters:
2261: + x - the vector
2262: . n - the number of vectors
2263: - a - location of pointer to arrays obtained from `VecGetArrays()`

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

2271:   Level: intermediate

2273: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2274: @*/
2275: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2276: {
2277:   PetscInt      i;
2278:   PetscScalar **q = *a;

2280:   PetscFunctionBegin;
2281:   PetscAssertPointer(x, 1);
2283:   PetscAssertPointer(a, 3);

2285:   for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2286:   PetscCall(PetscFree(q));
2287:   PetscFunctionReturn(PETSC_SUCCESS);
2288: }

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

2295:   Logically Collective; No Fortran Support

2297:   Input Parameter:
2298: . x - the vector

2300:   Output Parameters:
2301: + a     - location to put pointer to the array
2302: - mtype - memory type of the array

2304:   Level: beginner

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

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

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

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

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

2342:   Logically Collective; No Fortran Support

2344:   Input Parameters:
2345: + x - the vector
2346: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`

2348:   Level: beginner

2350: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`,
2351:           `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2352: @*/
2353: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar **a)
2354: {
2355:   PetscFunctionBegin;
2358:   if (a) PetscAssertPointer(a, 2);
2359:   if (x->ops->restorearrayandmemtype) {
2360:     /* VECCUDA, VECKOKKOS etc */
2361:     PetscUseTypeMethod(x, restorearrayandmemtype, a);
2362:   } else {
2363:     /* VECNEST, VECVIENNACL */
2364:     PetscCall(VecRestoreArray(x, a));
2365:   } /* VECSTANDARD does nothing */
2366:   if (a) *a = NULL;
2367:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
2368:   PetscFunctionReturn(PETSC_SUCCESS);
2369: }

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

2375:   Not Collective; No Fortran Support

2377:   Input Parameter:
2378: . x - the vector

2380:   Output Parameters:
2381: + a     - the array
2382: - mtype - memory type of the array

2384:   Level: beginner

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

2389: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2390: @*/
2391: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar **a, PetscMemType *mtype)
2392: {
2393:   PetscFunctionBegin;
2396:   PetscAssertPointer(a, 2);
2397:   if (mtype) PetscAssertPointer(mtype, 3);
2398:   if (x->ops->getarrayreadandmemtype) {
2399:     /* VECCUDA/VECHIP though they are also petscnative */
2400:     PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2401:   } else if (x->ops->getarrayandmemtype) {
2402:     /* VECKOKKOS */
2403:     PetscObjectState state;

2405:     // see VecGetArrayRead() for why
2406:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2407:     PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2408:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
2409:   } else {
2410:     PetscCall(VecGetArrayRead(x, a));
2411:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2412:   }
2413:   PetscFunctionReturn(PETSC_SUCCESS);
2414: }

2416: /*@C
2417:   VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`

2419:   Not Collective; No Fortran Support

2421:   Input Parameters:
2422: + x - the vector
2423: - a - the array

2425:   Level: beginner

2427: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2428: @*/
2429: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar **a)
2430: {
2431:   PetscFunctionBegin;
2434:   if (a) PetscAssertPointer(a, 2);
2435:   if (x->ops->restorearrayreadandmemtype) {
2436:     /* VECCUDA/VECHIP */
2437:     PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2438:   } else if (!x->petscnative) {
2439:     /* VECNEST */
2440:     PetscCall(VecRestoreArrayRead(x, a));
2441:   }
2442:   if (a) *a = NULL;
2443:   PetscFunctionReturn(PETSC_SUCCESS);
2444: }

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

2450:   Logically Collective; No Fortran Support

2452:   Input Parameter:
2453: . x - the vector

2455:   Output Parameters:
2456: + a     - the array
2457: - mtype - memory type of the array

2459:   Level: beginner

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

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

2487: /*@C
2488:   VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`

2490:   Logically Collective; No Fortran Support

2492:   Input Parameters:
2493: + x - the vector
2494: - a - the array

2496:   Level: beginner

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

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

2525:   Logically Collective; No Fortran Support

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

2531:   Level: developer

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

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

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

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

2546: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2547: @*/
2548: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2549: {
2550:   PetscFunctionBegin;
2553:   if (array) PetscAssertPointer(array, 2);
2554:   PetscUseTypeMethod(vec, placearray, array);
2555:   PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2556:   PetscFunctionReturn(PETSC_SUCCESS);
2557: }

2559: /*@C
2560:   VecReplaceArray - Allows one to replace the array in a vector with an
2561:   array provided by the user. This is useful to avoid copying an array
2562:   into a vector.

2564:   Logically Collective; No Fortran Support

2566:   Input Parameters:
2567: + vec   - the vector
2568: - array - the array

2570:   Level: developer

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

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

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

2582: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2583: @*/
2584: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2585: {
2586:   PetscFunctionBegin;
2589:   PetscUseTypeMethod(vec, replacearray, array);
2590:   PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2591:   PetscFunctionReturn(PETSC_SUCCESS);
2592: }

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

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

2601:     Collective

2603:     Input Parameters:
2604: +   x - a vector to mimic
2605: -   n - the number of vectors to obtain

2607:     Output Parameters:
2608: +   y - Fortran pointer to the array of vectors
2609: -   ierr - error code

2611:     Example of Usage:
2612: .vb
2613: #include <petsc/finclude/petscvec.h>
2614:     use petscvec

2616:     Vec x
2617:     Vec, pointer :: y(:)
2618:     ....
2619:     call VecDuplicateVecsF90(x,2,y,ierr)
2620:     call VecSet(y(2),alpha,ierr)
2621:     call VecSet(y(2),alpha,ierr)
2622:     ....
2623:     call VecDestroyVecsF90(2,y,ierr)
2624: .ve

2626:     Level: beginner

2628:     Note:
2629:     Use `VecDestroyVecsF90()` to free the space.

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

2634: /*MC
2635:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2636:     `VecGetArrayF90()`.

2638:     Synopsis:
2639:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2641:     Logically Collective

2643:     Input Parameters:
2644: +   x - vector
2645: -   xx_v - the Fortran pointer to the array

2647:     Output Parameter:
2648: .   ierr - error code

2650:     Example of Usage:
2651: .vb
2652: #include <petsc/finclude/petscvec.h>
2653:     use petscvec

2655:     PetscScalar, pointer :: xx_v(:)
2656:     ....
2657:     call VecGetArrayF90(x,xx_v,ierr)
2658:     xx_v(3) = a
2659:     call VecRestoreArrayF90(x,xx_v,ierr)
2660: .ve

2662:     Level: beginner

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

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

2670:     Synopsis:
2671:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

2673:     Collective

2675:     Input Parameters:
2676: +   n - the number of vectors previously obtained
2677: -   x - pointer to array of vector pointers

2679:     Output Parameter:
2680: .   ierr - error code

2682:     Level: beginner

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

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

2693:     Synopsis:
2694:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2696:     Logically Collective

2698:     Input Parameter:
2699: .   x - vector

2701:     Output Parameters:
2702: +   xx_v - the Fortran pointer to the array
2703: -   ierr - error code

2705:     Example of Usage:
2706: .vb
2707: #include <petsc/finclude/petscvec.h>
2708:     use petscvec

2710:     PetscScalar, pointer :: xx_v(:)
2711:     ....
2712:     call VecGetArrayF90(x,xx_v,ierr)
2713:     xx_v(3) = a
2714:     call VecRestoreArrayF90(x,xx_v,ierr)
2715: .ve

2717:      Level: beginner

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

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

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

2731:     Synopsis:
2732:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2734:     Logically Collective

2736:     Input Parameter:
2737: .   x - vector

2739:     Output Parameters:
2740: +   xx_v - the Fortran pointer to the array
2741: -   ierr - error code

2743:     Example of Usage:
2744: .vb
2745: #include <petsc/finclude/petscvec.h>
2746:     use petscvec

2748:     PetscScalar, pointer :: xx_v(:)
2749:     ....
2750:     call VecGetArrayReadF90(x,xx_v,ierr)
2751:     a = xx_v(3)
2752:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2753: .ve

2755:     Level: beginner

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

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

2763: /*MC
2764:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2765:     `VecGetArrayReadF90()`.

2767:     Synopsis:
2768:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2770:     Logically Collective

2772:     Input Parameters:
2773: +   x - vector
2774: -   xx_v - the Fortran pointer to the array

2776:     Output Parameter:
2777: .   ierr - error code

2779:     Example of Usage:
2780: .vb
2781: #include <petsc/finclude/petscvec.h>
2782:     use petscvec

2784:     PetscScalar, pointer :: xx_v(:)
2785:     ....
2786:     call VecGetArrayReadF90(x,xx_v,ierr)
2787:     a = xx_v(3)
2788:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2789: .ve

2791:     Level: beginner

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

2796: /*@C
2797:   VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2798:   processor's portion of the vector data.  You MUST call `VecRestoreArray2d()`
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 VecGetArray2d(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(VecGetArray(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:   VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2848:   processor's portion of the vector data.  You MUST call `VecRestoreArray2dWrite()`
2849:   when you no longer need access to the array.

2851:   Logically Collective

2853:   Input Parameters:
2854: + x      - the vector
2855: . m      - first dimension of two dimensional array
2856: . n      - second dimension of two dimensional array
2857: . mstart - first index you will use in first coordinate direction (often 0)
2858: - nstart - first index in the second coordinate direction (often 0)

2860:   Output Parameter:
2861: . a - location to put pointer to the array

2863:   Level: developer

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

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

2873: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2874:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2875:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2876: @*/
2877: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2878: {
2879:   PetscInt     i, N;
2880:   PetscScalar *aa;

2882:   PetscFunctionBegin;
2884:   PetscAssertPointer(a, 6);
2886:   PetscCall(VecGetLocalSize(x, &N));
2887:   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);
2888:   PetscCall(VecGetArrayWrite(x, &aa));

2890:   PetscCall(PetscMalloc1(m, a));
2891:   for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2892:   *a -= mstart;
2893:   PetscFunctionReturn(PETSC_SUCCESS);
2894: }

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

2899:   Logically Collective

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

2909:   Level: developer

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

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

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

2927:   PetscFunctionBegin;
2929:   PetscAssertPointer(a, 6);
2931:   dummy = (void *)(*a + mstart);
2932:   PetscCall(PetscFree(dummy));
2933:   PetscCall(VecRestoreArray(x, NULL));
2934:   *a = NULL;
2935:   PetscFunctionReturn(PETSC_SUCCESS);
2936: }

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

2941:   Logically Collective

2943:   Input Parameters:
2944: + x      - the vector
2945: . m      - first dimension of two dimensional array
2946: . n      - second dimension of the two dimensional array
2947: . mstart - first index you will use in first coordinate direction (often 0)
2948: . nstart - first index in the second coordinate direction (often 0)
2949: - a      - location of pointer to array obtained from `VecGetArray2d()`

2951:   Level: developer

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

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

2961: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2962:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2963:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2964: @*/
2965: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2966: {
2967:   void *dummy;

2969:   PetscFunctionBegin;
2971:   PetscAssertPointer(a, 6);
2973:   dummy = (void *)(*a + mstart);
2974:   PetscCall(PetscFree(dummy));
2975:   PetscCall(VecRestoreArrayWrite(x, NULL));
2976:   PetscFunctionReturn(PETSC_SUCCESS);
2977: }

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

2984:   Logically Collective

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

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

2994:   Level: developer

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

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

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

3011:   PetscFunctionBegin;
3013:   PetscAssertPointer(a, 4);
3015:   PetscCall(VecGetLocalSize(x, &N));
3016:   PetscCheck(m == N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local array size %" PetscInt_FMT " does not match 1d array dimensions %" PetscInt_FMT, N, m);
3017:   PetscCall(VecGetArray(x, a));
3018:   *a -= mstart;
3019:   PetscFunctionReturn(PETSC_SUCCESS);
3020: }

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

3027:   Logically Collective

3029:   Input Parameters:
3030: + x      - the vector
3031: . m      - first dimension of two dimensional array
3032: - mstart - first index you will use in first coordinate direction (often 0)

3034:   Output Parameter:
3035: . a - location to put pointer to the array

3037:   Level: developer

3039:   Notes:
3040:   For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3041:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3042:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

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

3046: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3047:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3048:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3049: @*/
3050: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3051: {
3052:   PetscInt N;

3054:   PetscFunctionBegin;
3056:   PetscAssertPointer(a, 4);
3058:   PetscCall(VecGetLocalSize(x, &N));
3059:   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);
3060:   PetscCall(VecGetArrayWrite(x, a));
3061:   *a -= mstart;
3062:   PetscFunctionReturn(PETSC_SUCCESS);
3063: }

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

3068:   Logically Collective

3070:   Input Parameters:
3071: + x      - the vector
3072: . m      - first dimension of two dimensional array
3073: . mstart - first index you will use in first coordinate direction (often 0)
3074: - a      - location of pointer to array obtained from `VecGetArray1d()`

3076:   Level: developer

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

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

3086: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3087:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3088:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3089: @*/
3090: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3091: {
3092:   PetscFunctionBegin;
3095:   PetscCall(VecRestoreArray(x, NULL));
3096:   *a = NULL;
3097:   PetscFunctionReturn(PETSC_SUCCESS);
3098: }

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

3103:   Logically Collective

3105:   Input Parameters:
3106: + x      - the vector
3107: . m      - first dimension of two dimensional array
3108: . mstart - first index you will use in first coordinate direction (often 0)
3109: - a      - location of pointer to array obtained from `VecGetArray1d()`

3111:   Level: developer

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

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

3121: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3122:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3123:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3124: @*/
3125: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3126: {
3127:   PetscFunctionBegin;
3130:   PetscCall(VecRestoreArrayWrite(x, NULL));
3131:   *a = NULL;
3132:   PetscFunctionReturn(PETSC_SUCCESS);
3133: }

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

3140:   Logically Collective

3142:   Input Parameters:
3143: + x      - the vector
3144: . m      - first dimension of three dimensional array
3145: . n      - second dimension of three dimensional array
3146: . p      - third dimension of three dimensional array
3147: . mstart - first index you will use in first coordinate direction (often 0)
3148: . nstart - first index in the second coordinate direction (often 0)
3149: - pstart - first index in the third coordinate direction (often 0)

3151:   Output Parameter:
3152: . a - location to put pointer to the array

3154:   Level: developer

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

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

3164: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3165:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
3166:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3167: @*/
3168: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3169: {
3170:   PetscInt     i, N, j;
3171:   PetscScalar *aa, **b;

3173:   PetscFunctionBegin;
3175:   PetscAssertPointer(a, 8);
3177:   PetscCall(VecGetLocalSize(x, &N));
3178:   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);
3179:   PetscCall(VecGetArray(x, &aa));

3181:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3182:   b = (PetscScalar **)((*a) + m);
3183:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3184:   for (i = 0; i < m; i++)
3185:     for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset(aa, i * n * p + j * p - pstart);
3186:   *a -= mstart;
3187:   PetscFunctionReturn(PETSC_SUCCESS);
3188: }

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

3195:   Logically Collective

3197:   Input Parameters:
3198: + x      - the vector
3199: . m      - first dimension of three dimensional array
3200: . n      - second dimension of three dimensional array
3201: . p      - third dimension of three dimensional array
3202: . mstart - first index you will use in first coordinate direction (often 0)
3203: . nstart - first index in the second coordinate direction (often 0)
3204: - pstart - first index in the third coordinate direction (often 0)

3206:   Output Parameter:
3207: . a - location to put pointer to the array

3209:   Level: developer

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

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

3219: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3220:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3221:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3222: @*/
3223: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3224: {
3225:   PetscInt     i, N, j;
3226:   PetscScalar *aa, **b;

3228:   PetscFunctionBegin;
3230:   PetscAssertPointer(a, 8);
3232:   PetscCall(VecGetLocalSize(x, &N));
3233:   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);
3234:   PetscCall(VecGetArrayWrite(x, &aa));

3236:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3237:   b = (PetscScalar **)((*a) + m);
3238:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3239:   for (i = 0; i < m; i++)
3240:     for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;

3242:   *a -= mstart;
3243:   PetscFunctionReturn(PETSC_SUCCESS);
3244: }

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

3249:   Logically Collective

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

3261:   Level: developer

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

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

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

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

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

3293:   Logically Collective

3295:   Input Parameters:
3296: + x      - the vector
3297: . m      - first dimension of three dimensional array
3298: . n      - second dimension of the three dimensional array
3299: . p      - third dimension of the three dimensional array
3300: . mstart - first index you will use in first coordinate direction (often 0)
3301: . nstart - first index in the second coordinate direction (often 0)
3302: . pstart - first index in the third coordinate direction (often 0)
3303: - a      - location of pointer to array obtained from VecGetArray3d()

3305:   Level: developer

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

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

3315: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3316:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3317:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3318: @*/
3319: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3320: {
3321:   void *dummy;

3323:   PetscFunctionBegin;
3325:   PetscAssertPointer(a, 8);
3327:   dummy = (void *)(*a + mstart);
3328:   PetscCall(PetscFree(dummy));
3329:   PetscCall(VecRestoreArrayWrite(x, NULL));
3330:   *a = NULL;
3331:   PetscFunctionReturn(PETSC_SUCCESS);
3332: }

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

3339:   Logically Collective

3341:   Input Parameters:
3342: + x      - the vector
3343: . m      - first dimension of four dimensional array
3344: . n      - second dimension of four dimensional array
3345: . p      - third dimension of four dimensional array
3346: . q      - fourth dimension of four dimensional array
3347: . mstart - first index you will use in first coordinate direction (often 0)
3348: . nstart - first index in the second coordinate direction (often 0)
3349: . pstart - first index in the third coordinate direction (often 0)
3350: - qstart - first index in the fourth coordinate direction (often 0)

3352:   Output Parameter:
3353: . a - location to put pointer to the array

3355:   Level: developer

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

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

3365: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3366:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3367:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3368: @*/
3369: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3370: {
3371:   PetscInt     i, N, j, k;
3372:   PetscScalar *aa, ***b, **c;

3374:   PetscFunctionBegin;
3376:   PetscAssertPointer(a, 10);
3378:   PetscCall(VecGetLocalSize(x, &N));
3379:   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);
3380:   PetscCall(VecGetArray(x, &aa));

3382:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3383:   b = (PetscScalar ***)((*a) + m);
3384:   c = (PetscScalar **)(b + m * n);
3385:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3386:   for (i = 0; i < m; i++)
3387:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3388:   for (i = 0; i < m; i++)
3389:     for (j = 0; j < n; j++)
3390:       for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3391:   *a -= mstart;
3392:   PetscFunctionReturn(PETSC_SUCCESS);
3393: }

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

3400:   Logically Collective

3402:   Input Parameters:
3403: + x      - the vector
3404: . m      - first dimension of four dimensional array
3405: . n      - second dimension of four dimensional array
3406: . p      - third dimension of four dimensional array
3407: . q      - fourth dimension of four dimensional array
3408: . mstart - first index you will use in first coordinate direction (often 0)
3409: . nstart - first index in the second coordinate direction (often 0)
3410: . pstart - first index in the third coordinate direction (often 0)
3411: - qstart - first index in the fourth coordinate direction (often 0)

3413:   Output Parameter:
3414: . a - location to put pointer to the array

3416:   Level: developer

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

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

3426: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3427:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3428:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3429: @*/
3430: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3431: {
3432:   PetscInt     i, N, j, k;
3433:   PetscScalar *aa, ***b, **c;

3435:   PetscFunctionBegin;
3437:   PetscAssertPointer(a, 10);
3439:   PetscCall(VecGetLocalSize(x, &N));
3440:   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);
3441:   PetscCall(VecGetArrayWrite(x, &aa));

3443:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3444:   b = (PetscScalar ***)((*a) + m);
3445:   c = (PetscScalar **)(b + m * n);
3446:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3447:   for (i = 0; i < m; i++)
3448:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3449:   for (i = 0; i < m; i++)
3450:     for (j = 0; j < n; j++)
3451:       for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3452:   *a -= mstart;
3453:   PetscFunctionReturn(PETSC_SUCCESS);
3454: }

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

3459:   Logically Collective

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

3473:   Level: developer

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

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

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

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

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

3505:   Logically Collective

3507:   Input Parameters:
3508: + x      - the vector
3509: . m      - first dimension of four dimensional array
3510: . n      - second dimension of the four dimensional array
3511: . p      - third dimension of the four dimensional array
3512: . q      - fourth dimension of the four dimensional array
3513: . mstart - first index you will use in first coordinate direction (often 0)
3514: . nstart - first index in the second coordinate direction (often 0)
3515: . pstart - first index in the third coordinate direction (often 0)
3516: . qstart - first index in the fourth coordinate direction (often 0)
3517: - a      - location of pointer to array obtained from `VecGetArray4d()`

3519:   Level: developer

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

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

3529: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3530:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3531:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3532: @*/
3533: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3534: {
3535:   void *dummy;

3537:   PetscFunctionBegin;
3539:   PetscAssertPointer(a, 10);
3541:   dummy = (void *)(*a + mstart);
3542:   PetscCall(PetscFree(dummy));
3543:   PetscCall(VecRestoreArrayWrite(x, NULL));
3544:   *a = NULL;
3545:   PetscFunctionReturn(PETSC_SUCCESS);
3546: }

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

3553:   Logically Collective

3555:   Input Parameters:
3556: + x      - the vector
3557: . m      - first dimension of two dimensional array
3558: . n      - second dimension of two dimensional array
3559: . mstart - first index you will use in first coordinate direction (often 0)
3560: - nstart - first index in the second coordinate direction (often 0)

3562:   Output Parameter:
3563: . a - location to put pointer to the array

3565:   Level: developer

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

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

3575: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3576:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3577:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3578: @*/
3579: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3580: {
3581:   PetscInt           i, N;
3582:   const PetscScalar *aa;

3584:   PetscFunctionBegin;
3586:   PetscAssertPointer(a, 6);
3588:   PetscCall(VecGetLocalSize(x, &N));
3589:   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);
3590:   PetscCall(VecGetArrayRead(x, &aa));

3592:   PetscCall(PetscMalloc1(m, a));
3593:   for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3594:   *a -= mstart;
3595:   PetscFunctionReturn(PETSC_SUCCESS);
3596: }

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

3601:   Logically Collective

3603:   Input Parameters:
3604: + x      - the vector
3605: . m      - first dimension of two dimensional array
3606: . n      - second dimension of the two dimensional array
3607: . mstart - first index you will use in first coordinate direction (often 0)
3608: . nstart - first index in the second coordinate direction (often 0)
3609: - a      - location of pointer to array obtained from VecGetArray2d()

3611:   Level: developer

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

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

3621: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3622:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3623:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3624: @*/
3625: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3626: {
3627:   void *dummy;

3629:   PetscFunctionBegin;
3631:   PetscAssertPointer(a, 6);
3633:   dummy = (void *)(*a + mstart);
3634:   PetscCall(PetscFree(dummy));
3635:   PetscCall(VecRestoreArrayRead(x, NULL));
3636:   *a = NULL;
3637:   PetscFunctionReturn(PETSC_SUCCESS);
3638: }

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

3645:   Logically Collective

3647:   Input Parameters:
3648: + x      - the vector
3649: . m      - first dimension of two dimensional array
3650: - mstart - first index you will use in first coordinate direction (often 0)

3652:   Output Parameter:
3653: . a - location to put pointer to the array

3655:   Level: developer

3657:   Notes:
3658:   For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3659:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3660:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

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

3664: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3665:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3666:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3667: @*/
3668: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3669: {
3670:   PetscInt N;

3672:   PetscFunctionBegin;
3674:   PetscAssertPointer(a, 4);
3676:   PetscCall(VecGetLocalSize(x, &N));
3677:   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);
3678:   PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3679:   *a -= mstart;
3680:   PetscFunctionReturn(PETSC_SUCCESS);
3681: }

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

3686:   Logically Collective

3688:   Input Parameters:
3689: + x      - the vector
3690: . m      - first dimension of two dimensional array
3691: . mstart - first index you will use in first coordinate direction (often 0)
3692: - a      - location of pointer to array obtained from `VecGetArray1dRead()`

3694:   Level: developer

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

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

3704: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3705:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3706:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3707: @*/
3708: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3709: {
3710:   PetscFunctionBegin;
3713:   PetscCall(VecRestoreArrayRead(x, NULL));
3714:   *a = NULL;
3715:   PetscFunctionReturn(PETSC_SUCCESS);
3716: }

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

3723:   Logically Collective

3725:   Input Parameters:
3726: + x      - the vector
3727: . m      - first dimension of three dimensional array
3728: . n      - second dimension of three dimensional array
3729: . p      - third dimension of three dimensional array
3730: . mstart - first index you will use in first coordinate direction (often 0)
3731: . nstart - first index in the second coordinate direction (often 0)
3732: - pstart - first index in the third coordinate direction (often 0)

3734:   Output Parameter:
3735: . a - location to put pointer to the array

3737:   Level: developer

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

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

3747: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3748:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3749:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3750: @*/
3751: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3752: {
3753:   PetscInt           i, N, j;
3754:   const PetscScalar *aa;
3755:   PetscScalar      **b;

3757:   PetscFunctionBegin;
3759:   PetscAssertPointer(a, 8);
3761:   PetscCall(VecGetLocalSize(x, &N));
3762:   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);
3763:   PetscCall(VecGetArrayRead(x, &aa));

3765:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3766:   b = (PetscScalar **)((*a) + m);
3767:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3768:   for (i = 0; i < m; i++)
3769:     for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset((PetscScalar *)aa, i * n * p + j * p - pstart);
3770:   *a -= mstart;
3771:   PetscFunctionReturn(PETSC_SUCCESS);
3772: }

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

3777:   Logically Collective

3779:   Input Parameters:
3780: + x      - the vector
3781: . m      - first dimension of three dimensional array
3782: . n      - second dimension of the three dimensional array
3783: . p      - third dimension of the three dimensional array
3784: . mstart - first index you will use in first coordinate direction (often 0)
3785: . nstart - first index in the second coordinate direction (often 0)
3786: . pstart - first index in the third coordinate direction (often 0)
3787: - a      - location of pointer to array obtained from `VecGetArray3dRead()`

3789:   Level: developer

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

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

3799: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3800:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3801:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3802: @*/
3803: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3804: {
3805:   void *dummy;

3807:   PetscFunctionBegin;
3809:   PetscAssertPointer(a, 8);
3811:   dummy = (void *)(*a + mstart);
3812:   PetscCall(PetscFree(dummy));
3813:   PetscCall(VecRestoreArrayRead(x, NULL));
3814:   *a = NULL;
3815:   PetscFunctionReturn(PETSC_SUCCESS);
3816: }

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

3823:   Logically Collective

3825:   Input Parameters:
3826: + x      - the vector
3827: . m      - first dimension of four dimensional array
3828: . n      - second dimension of four dimensional array
3829: . p      - third dimension of four dimensional array
3830: . q      - fourth dimension of four dimensional array
3831: . mstart - first index you will use in first coordinate direction (often 0)
3832: . nstart - first index in the second coordinate direction (often 0)
3833: . pstart - first index in the third coordinate direction (often 0)
3834: - qstart - first index in the fourth coordinate direction (often 0)

3836:   Output Parameter:
3837: . a - location to put pointer to the array

3839:   Level: beginner

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

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

3849: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3850:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3851:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3852: @*/
3853: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3854: {
3855:   PetscInt           i, N, j, k;
3856:   const PetscScalar *aa;
3857:   PetscScalar     ***b, **c;

3859:   PetscFunctionBegin;
3861:   PetscAssertPointer(a, 10);
3863:   PetscCall(VecGetLocalSize(x, &N));
3864:   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);
3865:   PetscCall(VecGetArrayRead(x, &aa));

3867:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3868:   b = (PetscScalar ***)((*a) + m);
3869:   c = (PetscScalar **)(b + m * n);
3870:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3871:   for (i = 0; i < m; i++)
3872:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3873:   for (i = 0; i < m; i++)
3874:     for (j = 0; j < n; j++)
3875:       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;
3876:   *a -= mstart;
3877:   PetscFunctionReturn(PETSC_SUCCESS);
3878: }

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

3883:   Logically Collective

3885:   Input Parameters:
3886: + x      - the vector
3887: . m      - first dimension of four dimensional array
3888: . n      - second dimension of the four dimensional array
3889: . p      - third dimension of the four dimensional array
3890: . q      - fourth dimension of the four dimensional array
3891: . mstart - first index you will use in first coordinate direction (often 0)
3892: . nstart - first index in the second coordinate direction (often 0)
3893: . pstart - first index in the third coordinate direction (often 0)
3894: . qstart - first index in the fourth coordinate direction (often 0)
3895: - a      - location of pointer to array obtained from `VecGetArray4dRead()`

3897:   Level: beginner

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

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

3907: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3908:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3909:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3910: @*/
3911: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3912: {
3913:   void *dummy;

3915:   PetscFunctionBegin;
3917:   PetscAssertPointer(a, 10);
3919:   dummy = (void *)(*a + mstart);
3920:   PetscCall(PetscFree(dummy));
3921:   PetscCall(VecRestoreArrayRead(x, NULL));
3922:   *a = NULL;
3923:   PetscFunctionReturn(PETSC_SUCCESS);
3924: }

3926: /*@
3927:   VecLockGet - Get the current lock status of a vector

3929:   Logically Collective

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

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

3938:   Level: advanced

3940: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3941: @*/
3942: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3943: {
3944:   PetscFunctionBegin;
3946:   PetscAssertPointer(state, 2);
3947:   *state = x->lock;
3948:   PetscFunctionReturn(PETSC_SUCCESS);
3949: }

3951: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3952: {
3953:   PetscFunctionBegin;
3955:   PetscAssertPointer(file, 2);
3956:   PetscAssertPointer(func, 3);
3957:   PetscAssertPointer(line, 4);
3958: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3959:   {
3960:     const int index = x->lockstack.currentsize - 1;

3962:     *file = index < 0 ? NULL : x->lockstack.file[index];
3963:     *func = index < 0 ? NULL : x->lockstack.function[index];
3964:     *line = index < 0 ? 0 : x->lockstack.line[index];
3965:   }
3966: #else
3967:   *file = NULL;
3968:   *func = NULL;
3969:   *line = 0;
3970: #endif
3971:   PetscFunctionReturn(PETSC_SUCCESS);
3972: }

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

3977:   Logically Collective

3979:   Input Parameter:
3980: . x - the vector

3982:   Level: intermediate

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

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

3990: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3991: @*/
3992: PetscErrorCode VecLockReadPush(Vec x)
3993: {
3994:   PetscFunctionBegin;
3996:   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");
3997: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3998:   {
3999:     const char *file, *func;
4000:     int         index, line;

4002:     if ((index = petscstack.currentsize - 2) < 0) {
4003:       // vec was locked "outside" of petsc, either in user-land or main. the error message will
4004:       // now show this function as the culprit, but it will include the stacktrace
4005:       file = "unknown user-file";
4006:       func = "unknown_user_function";
4007:       line = 0;
4008:     } else {
4009:       file = petscstack.file[index];
4010:       func = petscstack.function[index];
4011:       line = petscstack.line[index];
4012:     }
4013:     PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
4014:   }
4015: #endif
4016:   PetscFunctionReturn(PETSC_SUCCESS);
4017: }

4019: /*@
4020:   VecLockReadPop - Pop a read-only lock from a vector

4022:   Logically Collective

4024:   Input Parameter:
4025: . x - the vector

4027:   Level: intermediate

4029: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
4030: @*/
4031: PetscErrorCode VecLockReadPop(Vec x)
4032: {
4033:   PetscFunctionBegin;
4035:   PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
4036: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
4037:   {
4038:     const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];

4040:     PetscStackPop_Private(x->lockstack, previous);
4041:   }
4042: #endif
4043:   PetscFunctionReturn(PETSC_SUCCESS);
4044: }

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

4049:   Logically Collective

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

4055:   Level: intermediate

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

4063: .vb
4064:        VecGetArray(x,&xdata); // begin phase
4065:        VecLockWriteSet(v,PETSC_TRUE);

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

4069:        VecRestoreArray(x,&vdata); // end phase
4070:        VecLockWriteSet(v,PETSC_FALSE);
4071: .ve

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

4076: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
4077: @*/
4078: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
4079: {
4080:   PetscFunctionBegin;
4082:   if (flg) {
4083:     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");
4084:     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");
4085:     x->lock = -1;
4086:   } else {
4087:     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");
4088:     x->lock = 0;
4089:   }
4090:   PetscFunctionReturn(PETSC_SUCCESS);
4091: }