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:       PetscCheck(-b2[0] == b2[1], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Difference in cached %s norms: local %g", NormTypes[type], (double)*val);
221:     }
222:   }
223:   if (flg) PetscFunctionReturn(PETSC_SUCCESS);

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

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

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

238:   Not Collective

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

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

250:   Level: intermediate

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

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

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

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

279:   Collective

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

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

287:   Level: intermediate

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

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

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

316:   Collective

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

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

325:   Level: intermediate

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

330:   Returns the smallest index with the maximum value

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

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

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

356:   Collective

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

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

365:   Level: intermediate

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

370:   This returns the smallest index with the minimum value

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

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

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

397:   Collective

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

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

406:   Level: intermediate

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

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

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

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

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

448:   PetscFunctionBegin;
451:   VecCheckAssembled(x);
452:   PetscCall(VecSetErrorIfLocked(x, 1));
454:   if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);

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

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

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

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

475:   Logically Collective

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

481:   Level: intermediate

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

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

495: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
496: {
497:   PetscFunctionBegin;
500:   VecCheckAssembled(x);
502:   PetscCall(VecSetErrorIfLocked(x, 1));

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

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

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

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

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

542:   Logically Collective

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

548:   Level: beginner

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

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

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

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

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

599:   Logically Collective

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

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

609:   Level: intermediate

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

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

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

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

664:   Logically Collective

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

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

674:   Level: intermediate

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

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

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

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

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

720:   Logically Collective

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

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

731:   Level: intermediate

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

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

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

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

782:   Logically Collective

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

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

795:   Level: intermediate

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

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

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

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

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

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

850:   Logically Collective

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

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

860:   Level: intermediate

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

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

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

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

880:   Not Collective

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

889:   Level: beginner

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

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

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

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

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

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

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

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

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

940:   Not Collective

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

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

950:   Level: beginner

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

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

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

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

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

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

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

983:   Not Collective

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

992:   Level: intermediate

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

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

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

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

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

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

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

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

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

1041:   Not Collective

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

1050:   Level: intermediate

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

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

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

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

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

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

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

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

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

1103:   Not Collective

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

1112:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

1194:   Collective

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

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

1204:   Level: intermediate

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

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

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

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

1228:   Collective

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

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

1238:   Level: intermediate

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

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

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

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

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

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

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

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

1299:   Logically Collective

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

1307:   Level: intermediate

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

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

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

1324:   Logically Collective

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

1333:   Level: intermediate

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

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

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

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

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

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

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

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

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

1396:   Collective

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

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

1406:   Level: advanced

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

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

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

1428:   PetscFunctionBegin;
1432:   PetscAssertPointer(Y, 3);

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

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

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

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

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

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

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

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

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

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

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

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

1548:   Collective

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

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

1557:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

1677:   Collective

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

1684:   Level: advanced

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

1693:   PetscFunctionBegin;
1696:   PetscCheckSameComm(X, 1, is, 2);
1697:   PetscAssertPointer(Y, 3);

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

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

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

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

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

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

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

1787:   Not Collective.

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

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

1795:   Level: beginner

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

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

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

1825:   Not Collective.

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

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

1833:   Level: beginner

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

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

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

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

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

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

1878:   Not Collective.

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

1884:   Level: beginner

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

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

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

1912:   Collective

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

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

1920:   Level: beginner

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

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

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

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

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

1960:   Logically Collective.

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

1966:   Level: beginner

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

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

1992:   Logically Collective

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

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

2000:   Level: beginner

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

2008:   Fortran Note:
2009: .vb
2010:   PetscScalar, pointer :: a(:)
2011: .ve

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

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

2032:   Logically Collective

2034:   Input Parameters:
2035: + x - the vector
2036: - a - location of pointer to array obtained from `VecGetArray()`

2038:   Level: beginner

2040: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `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: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2078: @*/
2079: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar *a[])
2080: {
2081:   PetscFunctionBegin;
2083:   PetscAssertPointer(a, 2);
2084:   if (x->ops->getarrayread) {
2085:     PetscUseTypeMethod(x, getarrayread, a);
2086:   } else if (x->ops->getarray) {
2087:     PetscObjectState state;

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

2102: /*@C
2103:   VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`

2105:   Not Collective

2107:   Input Parameters:
2108: + x - the vector
2109: - a - the array

2111:   Level: beginner

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

2127:     // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2128:     // we can just undo that
2129:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2130:     PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2131:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
2132:   }
2133:   if (a) *a = NULL;
2134:   PetscFunctionReturn(PETSC_SUCCESS);
2135: }

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

2141:   Logically Collective

2143:   Input Parameter:
2144: . x - the vector

2146:   Output Parameter:
2147: . a - location to put pointer to the array

2149:   Level: intermediate

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

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

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

2160: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecPlaceArray()`, `VecGetArray2d()`,
2161:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
2162: @*/
2163: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar *a[])
2164: {
2165:   PetscFunctionBegin;
2167:   PetscAssertPointer(a, 2);
2168:   PetscCall(VecSetErrorIfLocked(x, 1));
2169:   if (x->ops->getarraywrite) {
2170:     PetscUseTypeMethod(x, getarraywrite, a);
2171:   } else {
2172:     PetscCall(VecGetArray(x, a));
2173:   }
2174:   PetscFunctionReturn(PETSC_SUCCESS);
2175: }

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

2180:   Logically Collective

2182:   Input Parameters:
2183: + x - the vector
2184: - a - location of pointer to array obtained from `VecGetArray()`

2186:   Level: beginner

2188: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2189:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2190: @*/
2191: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar *a[])
2192: {
2193:   PetscFunctionBegin;
2195:   if (a) PetscAssertPointer(a, 2);
2196:   if (x->ops->restorearraywrite) {
2197:     PetscUseTypeMethod(x, restorearraywrite, a);
2198:   } else if (x->ops->restorearray) {
2199:     PetscUseTypeMethod(x, restorearray, a);
2200:   }
2201:   if (a) *a = NULL;
2202:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
2203:   PetscFunctionReturn(PETSC_SUCCESS);
2204: }

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

2210:   Logically Collective; No Fortran Support

2212:   Input Parameters:
2213: + x - the vectors
2214: - n - the number of vectors

2216:   Output Parameter:
2217: . a - location to put pointer to the array

2219:   Level: intermediate

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

2224: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2225: @*/
2226: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2227: {
2228:   PetscInt      i;
2229:   PetscScalar **q;

2231:   PetscFunctionBegin;
2232:   PetscAssertPointer(x, 1);
2234:   PetscAssertPointer(a, 3);
2235:   PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2236:   PetscCall(PetscMalloc1(n, &q));
2237:   for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2238:   *a = q;
2239:   PetscFunctionReturn(PETSC_SUCCESS);
2240: }

2242: /*@C
2243:   VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2244:   has been called.

2246:   Logically Collective; No Fortran Support

2248:   Input Parameters:
2249: + x - the vector
2250: . n - the number of vectors
2251: - a - location of pointer to arrays obtained from `VecGetArrays()`

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

2259:   Level: intermediate

2261: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2262: @*/
2263: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2264: {
2265:   PetscInt      i;
2266:   PetscScalar **q = *a;

2268:   PetscFunctionBegin;
2269:   PetscAssertPointer(x, 1);
2271:   PetscAssertPointer(a, 3);

2273:   for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2274:   PetscCall(PetscFree(q));
2275:   PetscFunctionReturn(PETSC_SUCCESS);
2276: }

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

2283:   Logically Collective; No Fortran Support

2285:   Input Parameter:
2286: . x - the vector

2288:   Output Parameters:
2289: + a     - location to put pointer to the array
2290: - mtype - memory type of the array

2292:   Level: beginner

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

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

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

2305: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`,
2306:           `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2307: @*/
2308: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar *a[], PetscMemType *mtype)
2309: {
2310:   PetscFunctionBegin;
2313:   PetscAssertPointer(a, 2);
2314:   if (mtype) PetscAssertPointer(mtype, 3);
2315:   PetscCall(VecSetErrorIfLocked(x, 1));
2316:   if (x->ops->getarrayandmemtype) {
2317:     /* VECCUDA, VECKOKKOS etc */
2318:     PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2319:   } else {
2320:     /* VECSTANDARD, VECNEST, VECVIENNACL */
2321:     PetscCall(VecGetArray(x, a));
2322:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2323:   }
2324:   PetscFunctionReturn(PETSC_SUCCESS);
2325: }

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

2330:   Logically Collective; No Fortran Support

2332:   Input Parameters:
2333: + x - the vector
2334: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`

2336:   Level: beginner

2338: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`,
2339:           `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2340: @*/
2341: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar *a[])
2342: {
2343:   PetscFunctionBegin;
2346:   if (a) PetscAssertPointer(a, 2);
2347:   if (x->ops->restorearrayandmemtype) {
2348:     /* VECCUDA, VECKOKKOS etc */
2349:     PetscUseTypeMethod(x, restorearrayandmemtype, a);
2350:   } else {
2351:     /* VECNEST, VECVIENNACL */
2352:     PetscCall(VecRestoreArray(x, a));
2353:   } /* VECSTANDARD does nothing */
2354:   if (a) *a = NULL;
2355:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
2356:   PetscFunctionReturn(PETSC_SUCCESS);
2357: }

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

2363:   Not Collective; No Fortran Support

2365:   Input Parameter:
2366: . x - the vector

2368:   Output Parameters:
2369: + a     - the array
2370: - mtype - memory type of the array

2372:   Level: beginner

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

2377: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2378: @*/
2379: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar *a[], PetscMemType *mtype)
2380: {
2381:   PetscFunctionBegin;
2384:   PetscAssertPointer(a, 2);
2385:   if (mtype) PetscAssertPointer(mtype, 3);
2386:   if (x->ops->getarrayreadandmemtype) {
2387:     /* VECCUDA/VECHIP though they are also petscnative */
2388:     PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2389:   } else if (x->ops->getarrayandmemtype) {
2390:     /* VECKOKKOS */
2391:     PetscObjectState state;

2393:     // see VecGetArrayRead() for why
2394:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2395:     PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2396:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
2397:   } else {
2398:     PetscCall(VecGetArrayRead(x, a));
2399:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2400:   }
2401:   PetscFunctionReturn(PETSC_SUCCESS);
2402: }

2404: /*@C
2405:   VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`

2407:   Not Collective; No Fortran Support

2409:   Input Parameters:
2410: + x - the vector
2411: - a - the array

2413:   Level: beginner

2415: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2416: @*/
2417: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar *a[])
2418: {
2419:   PetscFunctionBegin;
2422:   if (a) PetscAssertPointer(a, 2);
2423:   if (x->ops->restorearrayreadandmemtype) {
2424:     /* VECCUDA/VECHIP */
2425:     PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2426:   } else if (!x->petscnative) {
2427:     /* VECNEST */
2428:     PetscCall(VecRestoreArrayRead(x, a));
2429:   }
2430:   if (a) *a = NULL;
2431:   PetscFunctionReturn(PETSC_SUCCESS);
2432: }

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

2438:   Logically Collective; No Fortran Support

2440:   Input Parameter:
2441: . x - the vector

2443:   Output Parameters:
2444: + a     - the array
2445: - mtype - memory type of the array

2447:   Level: beginner

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

2452: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2453: @*/
2454: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar *a[], PetscMemType *mtype)
2455: {
2456:   PetscFunctionBegin;
2459:   PetscCall(VecSetErrorIfLocked(x, 1));
2460:   PetscAssertPointer(a, 2);
2461:   if (mtype) PetscAssertPointer(mtype, 3);
2462:   if (x->ops->getarraywriteandmemtype) {
2463:     /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2464:     PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2465:   } else if (x->ops->getarrayandmemtype) {
2466:     PetscCall(VecGetArrayAndMemType(x, a, mtype));
2467:   } else {
2468:     /* VECNEST, VECVIENNACL */
2469:     PetscCall(VecGetArrayWrite(x, a));
2470:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2471:   }
2472:   PetscFunctionReturn(PETSC_SUCCESS);
2473: }

2475: /*@C
2476:   VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`

2478:   Logically Collective; No Fortran Support

2480:   Input Parameters:
2481: + x - the vector
2482: - a - the array

2484:   Level: beginner

2486: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2487: @*/
2488: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar *a[])
2489: {
2490:   PetscFunctionBegin;
2493:   PetscCall(VecSetErrorIfLocked(x, 1));
2494:   if (a) PetscAssertPointer(a, 2);
2495:   if (x->ops->restorearraywriteandmemtype) {
2496:     /* VECCUDA/VECHIP */
2497:     PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2498:     PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2499:   } else if (x->ops->restorearrayandmemtype) {
2500:     PetscCall(VecRestoreArrayAndMemType(x, a));
2501:   } else {
2502:     PetscCall(VecRestoreArray(x, a));
2503:   }
2504:   if (a) *a = NULL;
2505:   PetscFunctionReturn(PETSC_SUCCESS);
2506: }

2508: /*@
2509:   VecPlaceArray - Allows one to replace the array in a vector with an
2510:   array provided by the user. This is useful to avoid copying an array
2511:   into a vector.

2513:   Logically Collective; No Fortran Support

2515:   Input Parameters:
2516: + vec   - the vector
2517: - array - the array

2519:   Level: developer

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

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

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

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

2534: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2535: @*/
2536: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2537: {
2538:   PetscFunctionBegin;
2541:   if (array) PetscAssertPointer(array, 2);
2542:   PetscUseTypeMethod(vec, placearray, array);
2543:   PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2544:   PetscFunctionReturn(PETSC_SUCCESS);
2545: }

2547: /*@C
2548:   VecReplaceArray - Allows one to replace the array in a vector with an
2549:   array provided by the user. This is useful to avoid copying an array
2550:   into a vector.

2552:   Logically Collective; No Fortran Support

2554:   Input Parameters:
2555: + vec   - the vector
2556: - array - the array

2558:   Level: developer

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

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

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

2570: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2571: @*/
2572: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2573: {
2574:   PetscFunctionBegin;
2577:   PetscUseTypeMethod(vec, replacearray, array);
2578:   PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2579:   PetscFunctionReturn(PETSC_SUCCESS);
2580: }

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

2587:   Logically Collective

2589:   Input Parameters:
2590: + x      - the vector
2591: . m      - first dimension of two dimensional array
2592: . n      - second dimension of two dimensional array
2593: . mstart - first index you will use in first coordinate direction (often 0)
2594: - nstart - first index in the second coordinate direction (often 0)

2596:   Output Parameter:
2597: . a - location to put pointer to the array

2599:   Level: developer

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

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

2609: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2610:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2611:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2612: @*/
2613: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2614: {
2615:   PetscInt     i, N;
2616:   PetscScalar *aa;

2618:   PetscFunctionBegin;
2620:   PetscAssertPointer(a, 6);
2622:   PetscCall(VecGetLocalSize(x, &N));
2623:   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);
2624:   PetscCall(VecGetArray(x, &aa));

2626:   PetscCall(PetscMalloc1(m, a));
2627:   for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2628:   *a -= mstart;
2629:   PetscFunctionReturn(PETSC_SUCCESS);
2630: }

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

2637:   Logically Collective

2639:   Input Parameters:
2640: + x      - the vector
2641: . m      - first dimension of two dimensional array
2642: . n      - second dimension of two dimensional array
2643: . mstart - first index you will use in first coordinate direction (often 0)
2644: - nstart - first index in the second coordinate direction (often 0)

2646:   Output Parameter:
2647: . a - location to put pointer to the array

2649:   Level: developer

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

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

2659: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2660:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2661:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2662: @*/
2663: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2664: {
2665:   PetscInt     i, N;
2666:   PetscScalar *aa;

2668:   PetscFunctionBegin;
2670:   PetscAssertPointer(a, 6);
2672:   PetscCall(VecGetLocalSize(x, &N));
2673:   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);
2674:   PetscCall(VecGetArrayWrite(x, &aa));

2676:   PetscCall(PetscMalloc1(m, a));
2677:   for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2678:   *a -= mstart;
2679:   PetscFunctionReturn(PETSC_SUCCESS);
2680: }

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

2685:   Logically Collective

2687:   Input Parameters:
2688: + x      - the vector
2689: . m      - first dimension of two dimensional array
2690: . n      - second dimension of the two dimensional array
2691: . mstart - first index you will use in first coordinate direction (often 0)
2692: . nstart - first index in the second coordinate direction (often 0)
2693: - a      - location of pointer to array obtained from `VecGetArray2d()`

2695:   Level: developer

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

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

2705: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2706:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2707:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2708: @*/
2709: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2710: {
2711:   void *dummy;

2713:   PetscFunctionBegin;
2715:   PetscAssertPointer(a, 6);
2717:   dummy = (void *)(*a + mstart);
2718:   PetscCall(PetscFree(dummy));
2719:   PetscCall(VecRestoreArray(x, NULL));
2720:   *a = NULL;
2721:   PetscFunctionReturn(PETSC_SUCCESS);
2722: }

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

2727:   Logically Collective

2729:   Input Parameters:
2730: + x      - the vector
2731: . m      - first dimension of two dimensional array
2732: . n      - second dimension of the two dimensional array
2733: . mstart - first index you will use in first coordinate direction (often 0)
2734: . nstart - first index in the second coordinate direction (often 0)
2735: - a      - location of pointer to array obtained from `VecGetArray2d()`

2737:   Level: developer

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

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

2747: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2748:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2749:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2750: @*/
2751: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2752: {
2753:   void *dummy;

2755:   PetscFunctionBegin;
2757:   PetscAssertPointer(a, 6);
2759:   dummy = (void *)(*a + mstart);
2760:   PetscCall(PetscFree(dummy));
2761:   PetscCall(VecRestoreArrayWrite(x, NULL));
2762:   PetscFunctionReturn(PETSC_SUCCESS);
2763: }

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

2770:   Logically Collective

2772:   Input Parameters:
2773: + x      - the vector
2774: . m      - first dimension of two dimensional array
2775: - mstart - first index you will use in first coordinate direction (often 0)

2777:   Output Parameter:
2778: . a - location to put pointer to the array

2780:   Level: developer

2782:   Notes:
2783:   For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2784:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2785:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

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

2789: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2790:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2791:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2792: @*/
2793: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2794: {
2795:   PetscInt N;

2797:   PetscFunctionBegin;
2799:   PetscAssertPointer(a, 4);
2801:   PetscCall(VecGetLocalSize(x, &N));
2802:   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);
2803:   PetscCall(VecGetArray(x, a));
2804:   *a -= mstart;
2805:   PetscFunctionReturn(PETSC_SUCCESS);
2806: }

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

2813:   Logically Collective

2815:   Input Parameters:
2816: + x      - the vector
2817: . m      - first dimension of two dimensional array
2818: - mstart - first index you will use in first coordinate direction (often 0)

2820:   Output Parameter:
2821: . a - location to put pointer to the array

2823:   Level: developer

2825:   Notes:
2826:   For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2827:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2828:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

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

2832: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2833:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2834:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2835: @*/
2836: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2837: {
2838:   PetscInt N;

2840:   PetscFunctionBegin;
2842:   PetscAssertPointer(a, 4);
2844:   PetscCall(VecGetLocalSize(x, &N));
2845:   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);
2846:   PetscCall(VecGetArrayWrite(x, a));
2847:   *a -= mstart;
2848:   PetscFunctionReturn(PETSC_SUCCESS);
2849: }

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

2854:   Logically Collective

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

2862:   Level: developer

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

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

2872: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2873:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2874:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2875: @*/
2876: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2877: {
2878:   PetscFunctionBegin;
2881:   PetscCall(VecRestoreArray(x, NULL));
2882:   *a = NULL;
2883:   PetscFunctionReturn(PETSC_SUCCESS);
2884: }

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

2889:   Logically Collective

2891:   Input Parameters:
2892: + x      - the vector
2893: . m      - first dimension of two dimensional array
2894: . mstart - first index you will use in first coordinate direction (often 0)
2895: - a      - location of pointer to array obtained from `VecGetArray1d()`

2897:   Level: developer

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

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

2907: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2908:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2909:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2910: @*/
2911: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2912: {
2913:   PetscFunctionBegin;
2916:   PetscCall(VecRestoreArrayWrite(x, NULL));
2917:   *a = NULL;
2918:   PetscFunctionReturn(PETSC_SUCCESS);
2919: }

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

2926:   Logically Collective

2928:   Input Parameters:
2929: + x      - the vector
2930: . m      - first dimension of three dimensional array
2931: . n      - second dimension of three dimensional array
2932: . p      - third dimension of three dimensional array
2933: . mstart - first index you will use in first coordinate direction (often 0)
2934: . nstart - first index in the second coordinate direction (often 0)
2935: - pstart - first index in the third coordinate direction (often 0)

2937:   Output Parameter:
2938: . a - location to put pointer to the array

2940:   Level: developer

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

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

2950: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2951:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
2952:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2953: @*/
2954: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
2955: {
2956:   PetscInt     i, N, j;
2957:   PetscScalar *aa, **b;

2959:   PetscFunctionBegin;
2961:   PetscAssertPointer(a, 8);
2963:   PetscCall(VecGetLocalSize(x, &N));
2964:   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);
2965:   PetscCall(VecGetArray(x, &aa));

2967:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
2968:   b = (PetscScalar **)((*a) + m);
2969:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
2970:   for (i = 0; i < m; i++)
2971:     for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset(aa, i * n * p + j * p - pstart);
2972:   *a -= mstart;
2973:   PetscFunctionReturn(PETSC_SUCCESS);
2974: }

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

2981:   Logically Collective

2983:   Input Parameters:
2984: + x      - the vector
2985: . m      - first dimension of three dimensional array
2986: . n      - second dimension of three dimensional array
2987: . p      - third dimension of three dimensional array
2988: . mstart - first index you will use in first coordinate direction (often 0)
2989: . nstart - first index in the second coordinate direction (often 0)
2990: - pstart - first index in the third coordinate direction (often 0)

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

2995:   Level: developer

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

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

3005: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3006:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3007:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3008: @*/
3009: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3010: {
3011:   PetscInt     i, N, j;
3012:   PetscScalar *aa, **b;

3014:   PetscFunctionBegin;
3016:   PetscAssertPointer(a, 8);
3018:   PetscCall(VecGetLocalSize(x, &N));
3019:   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);
3020:   PetscCall(VecGetArrayWrite(x, &aa));

3022:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3023:   b = (PetscScalar **)((*a) + m);
3024:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3025:   for (i = 0; i < m; i++)
3026:     for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;

3028:   *a -= mstart;
3029:   PetscFunctionReturn(PETSC_SUCCESS);
3030: }

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

3035:   Logically Collective

3037:   Input Parameters:
3038: + x      - the vector
3039: . m      - first dimension of three dimensional array
3040: . n      - second dimension of the three dimensional array
3041: . p      - third dimension of the three dimensional array
3042: . mstart - first index you will use in first coordinate direction (often 0)
3043: . nstart - first index in the second coordinate direction (often 0)
3044: . pstart - first index in the third coordinate direction (often 0)
3045: - a      - location of pointer to array obtained from VecGetArray3d()

3047:   Level: developer

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

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

3057: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3058:           `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3059:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3060: @*/
3061: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3062: {
3063:   void *dummy;

3065:   PetscFunctionBegin;
3067:   PetscAssertPointer(a, 8);
3069:   dummy = (void *)(*a + mstart);
3070:   PetscCall(PetscFree(dummy));
3071:   PetscCall(VecRestoreArray(x, NULL));
3072:   *a = NULL;
3073:   PetscFunctionReturn(PETSC_SUCCESS);
3074: }

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

3079:   Logically Collective

3081:   Input Parameters:
3082: + x      - the vector
3083: . m      - first dimension of three dimensional array
3084: . n      - second dimension of the three dimensional array
3085: . p      - third dimension of the three dimensional array
3086: . mstart - first index you will use in first coordinate direction (often 0)
3087: . nstart - first index in the second coordinate direction (often 0)
3088: . pstart - first index in the third coordinate direction (often 0)
3089: - a      - location of pointer to array obtained from VecGetArray3d()

3091:   Level: developer

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

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

3101: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3102:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3103:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3104: @*/
3105: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3106: {
3107:   void *dummy;

3109:   PetscFunctionBegin;
3111:   PetscAssertPointer(a, 8);
3113:   dummy = (void *)(*a + mstart);
3114:   PetscCall(PetscFree(dummy));
3115:   PetscCall(VecRestoreArrayWrite(x, NULL));
3116:   *a = NULL;
3117:   PetscFunctionReturn(PETSC_SUCCESS);
3118: }

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

3125:   Logically Collective

3127:   Input Parameters:
3128: + x      - the vector
3129: . m      - first dimension of four dimensional array
3130: . n      - second dimension of four dimensional array
3131: . p      - third dimension of four dimensional array
3132: . q      - fourth dimension of four dimensional array
3133: . mstart - first index you will use in first coordinate direction (often 0)
3134: . nstart - first index in the second coordinate direction (often 0)
3135: . pstart - first index in the third coordinate direction (often 0)
3136: - qstart - first index in the fourth coordinate direction (often 0)

3138:   Output Parameter:
3139: . a - location to put pointer to the array

3141:   Level: developer

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

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

3151: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3152:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3153:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3154: @*/
3155: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3156: {
3157:   PetscInt     i, N, j, k;
3158:   PetscScalar *aa, ***b, **c;

3160:   PetscFunctionBegin;
3162:   PetscAssertPointer(a, 10);
3164:   PetscCall(VecGetLocalSize(x, &N));
3165:   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);
3166:   PetscCall(VecGetArray(x, &aa));

3168:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3169:   b = (PetscScalar ***)((*a) + m);
3170:   c = (PetscScalar **)(b + m * n);
3171:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3172:   for (i = 0; i < m; i++)
3173:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3174:   for (i = 0; i < m; i++)
3175:     for (j = 0; j < n; j++)
3176:       for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3177:   *a -= mstart;
3178:   PetscFunctionReturn(PETSC_SUCCESS);
3179: }

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

3186:   Logically Collective

3188:   Input Parameters:
3189: + x      - the vector
3190: . m      - first dimension of four dimensional array
3191: . n      - second dimension of four dimensional array
3192: . p      - third dimension of four dimensional array
3193: . q      - fourth dimension of four dimensional array
3194: . mstart - first index you will use in first coordinate direction (often 0)
3195: . nstart - first index in the second coordinate direction (often 0)
3196: . pstart - first index in the third coordinate direction (often 0)
3197: - qstart - first index in the fourth coordinate direction (often 0)

3199:   Output Parameter:
3200: . a - location to put pointer to the array

3202:   Level: developer

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

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

3212: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3213:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3214:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3215: @*/
3216: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3217: {
3218:   PetscInt     i, N, j, k;
3219:   PetscScalar *aa, ***b, **c;

3221:   PetscFunctionBegin;
3223:   PetscAssertPointer(a, 10);
3225:   PetscCall(VecGetLocalSize(x, &N));
3226:   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);
3227:   PetscCall(VecGetArrayWrite(x, &aa));

3229:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3230:   b = (PetscScalar ***)((*a) + m);
3231:   c = (PetscScalar **)(b + m * n);
3232:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3233:   for (i = 0; i < m; i++)
3234:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3235:   for (i = 0; i < m; i++)
3236:     for (j = 0; j < n; j++)
3237:       for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3238:   *a -= mstart;
3239:   PetscFunctionReturn(PETSC_SUCCESS);
3240: }

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

3245:   Logically Collective

3247:   Input Parameters:
3248: + x      - the vector
3249: . m      - first dimension of four dimensional array
3250: . n      - second dimension of the four dimensional array
3251: . p      - third dimension of the four dimensional array
3252: . q      - fourth dimension of the four dimensional array
3253: . mstart - first index you will use in first coordinate direction (often 0)
3254: . nstart - first index in the second coordinate direction (often 0)
3255: . pstart - first index in the third coordinate direction (often 0)
3256: . qstart - first index in the fourth coordinate direction (often 0)
3257: - a      - location of pointer to array obtained from VecGetArray4d()

3259:   Level: developer

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

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

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

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

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

3291:   Logically Collective

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

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()`, `VecPlaceArray()`,
3316:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3317:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3318: @*/
3319: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3320: {
3321:   void *dummy;

3323:   PetscFunctionBegin;
3325:   PetscAssertPointer(a, 10);
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:   VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3336:   processor's portion of the vector data.  You MUST call `VecRestoreArray2dRead()`
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 two dimensional array
3344: . n      - second dimension of two dimensional array
3345: . mstart - first index you will use in first coordinate direction (often 0)
3346: - nstart - first index in the second coordinate direction (often 0)

3348:   Output Parameter:
3349: . a - location to put pointer to the array

3351:   Level: developer

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

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

3361: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3362:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3363:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3364: @*/
3365: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3366: {
3367:   PetscInt           i, N;
3368:   const PetscScalar *aa;

3370:   PetscFunctionBegin;
3372:   PetscAssertPointer(a, 6);
3374:   PetscCall(VecGetLocalSize(x, &N));
3375:   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);
3376:   PetscCall(VecGetArrayRead(x, &aa));

3378:   PetscCall(PetscMalloc1(m, a));
3379:   for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3380:   *a -= mstart;
3381:   PetscFunctionReturn(PETSC_SUCCESS);
3382: }

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

3387:   Logically Collective

3389:   Input Parameters:
3390: + x      - the vector
3391: . m      - first dimension of two dimensional array
3392: . n      - second dimension of the two dimensional array
3393: . mstart - first index you will use in first coordinate direction (often 0)
3394: . nstart - first index in the second coordinate direction (often 0)
3395: - a      - location of pointer to array obtained from VecGetArray2d()

3397:   Level: developer

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

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

3407: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3408:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3409:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3410: @*/
3411: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3412: {
3413:   void *dummy;

3415:   PetscFunctionBegin;
3417:   PetscAssertPointer(a, 6);
3419:   dummy = (void *)(*a + mstart);
3420:   PetscCall(PetscFree(dummy));
3421:   PetscCall(VecRestoreArrayRead(x, NULL));
3422:   *a = NULL;
3423:   PetscFunctionReturn(PETSC_SUCCESS);
3424: }

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

3431:   Logically Collective

3433:   Input Parameters:
3434: + x      - the vector
3435: . m      - first dimension of two dimensional array
3436: - mstart - first index you will use in first coordinate direction (often 0)

3438:   Output Parameter:
3439: . a - location to put pointer to the array

3441:   Level: developer

3443:   Notes:
3444:   For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3445:   obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3446:   `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

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

3450: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3451:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3452:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3453: @*/
3454: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3455: {
3456:   PetscInt N;

3458:   PetscFunctionBegin;
3460:   PetscAssertPointer(a, 4);
3462:   PetscCall(VecGetLocalSize(x, &N));
3463:   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);
3464:   PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3465:   *a -= mstart;
3466:   PetscFunctionReturn(PETSC_SUCCESS);
3467: }

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

3472:   Logically Collective

3474:   Input Parameters:
3475: + x      - the vector
3476: . m      - first dimension of two dimensional array
3477: . mstart - first index you will use in first coordinate direction (often 0)
3478: - a      - location of pointer to array obtained from `VecGetArray1dRead()`

3480:   Level: developer

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

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

3490: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3491:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3492:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3493: @*/
3494: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3495: {
3496:   PetscFunctionBegin;
3499:   PetscCall(VecRestoreArrayRead(x, NULL));
3500:   *a = NULL;
3501:   PetscFunctionReturn(PETSC_SUCCESS);
3502: }

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

3509:   Logically Collective

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

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

3523:   Level: developer

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

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

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

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

3551:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3552:   b = (PetscScalar **)((*a) + m);
3553:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3554:   for (i = 0; i < m; i++)
3555:     for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset((PetscScalar *)aa, i * n * p + j * p - pstart);
3556:   *a -= mstart;
3557:   PetscFunctionReturn(PETSC_SUCCESS);
3558: }

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

3563:   Logically Collective

3565:   Input Parameters:
3566: + x      - the vector
3567: . m      - first dimension of three dimensional array
3568: . n      - second dimension of the three dimensional array
3569: . p      - third dimension of the three dimensional array
3570: . mstart - first index you will use in first coordinate direction (often 0)
3571: . nstart - first index in the second coordinate direction (often 0)
3572: . pstart - first index in the third coordinate direction (often 0)
3573: - a      - location of pointer to array obtained from `VecGetArray3dRead()`

3575:   Level: developer

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

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

3585: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3586:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3587:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3588: @*/
3589: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3590: {
3591:   void *dummy;

3593:   PetscFunctionBegin;
3595:   PetscAssertPointer(a, 8);
3597:   dummy = (void *)(*a + mstart);
3598:   PetscCall(PetscFree(dummy));
3599:   PetscCall(VecRestoreArrayRead(x, NULL));
3600:   *a = NULL;
3601:   PetscFunctionReturn(PETSC_SUCCESS);
3602: }

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

3609:   Logically Collective

3611:   Input Parameters:
3612: + x      - the vector
3613: . m      - first dimension of four dimensional array
3614: . n      - second dimension of four dimensional array
3615: . p      - third dimension of four dimensional array
3616: . q      - fourth dimension of four dimensional array
3617: . mstart - first index you will use in first coordinate direction (often 0)
3618: . nstart - first index in the second coordinate direction (often 0)
3619: . pstart - first index in the third coordinate direction (often 0)
3620: - qstart - first index in the fourth coordinate direction (often 0)

3622:   Output Parameter:
3623: . a - location to put pointer to the array

3625:   Level: beginner

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

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

3635: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3636:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3637:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3638: @*/
3639: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3640: {
3641:   PetscInt           i, N, j, k;
3642:   const PetscScalar *aa;
3643:   PetscScalar     ***b, **c;

3645:   PetscFunctionBegin;
3647:   PetscAssertPointer(a, 10);
3649:   PetscCall(VecGetLocalSize(x, &N));
3650:   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);
3651:   PetscCall(VecGetArrayRead(x, &aa));

3653:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3654:   b = (PetscScalar ***)((*a) + m);
3655:   c = (PetscScalar **)(b + m * n);
3656:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3657:   for (i = 0; i < m; i++)
3658:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3659:   for (i = 0; i < m; i++)
3660:     for (j = 0; j < n; j++)
3661:       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;
3662:   *a -= mstart;
3663:   PetscFunctionReturn(PETSC_SUCCESS);
3664: }

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

3669:   Logically Collective

3671:   Input Parameters:
3672: + x      - the vector
3673: . m      - first dimension of four dimensional array
3674: . n      - second dimension of the four dimensional array
3675: . p      - third dimension of the four dimensional array
3676: . q      - fourth dimension of the four dimensional array
3677: . mstart - first index you will use in first coordinate direction (often 0)
3678: . nstart - first index in the second coordinate direction (often 0)
3679: . pstart - first index in the third coordinate direction (often 0)
3680: . qstart - first index in the fourth coordinate direction (often 0)
3681: - a      - location of pointer to array obtained from `VecGetArray4dRead()`

3683:   Level: beginner

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

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

3693: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3694:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3695:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3696: @*/
3697: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3698: {
3699:   void *dummy;

3701:   PetscFunctionBegin;
3703:   PetscAssertPointer(a, 10);
3705:   dummy = (void *)(*a + mstart);
3706:   PetscCall(PetscFree(dummy));
3707:   PetscCall(VecRestoreArrayRead(x, NULL));
3708:   *a = NULL;
3709:   PetscFunctionReturn(PETSC_SUCCESS);
3710: }

3712: /*@
3713:   VecLockGet - Get the current lock status of a vector

3715:   Logically Collective

3717:   Input Parameter:
3718: . x - the vector

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

3724:   Level: advanced

3726: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3727: @*/
3728: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3729: {
3730:   PetscFunctionBegin;
3732:   PetscAssertPointer(state, 2);
3733:   *state = x->lock;
3734:   PetscFunctionReturn(PETSC_SUCCESS);
3735: }

3737: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3738: {
3739:   PetscFunctionBegin;
3741:   PetscAssertPointer(file, 2);
3742:   PetscAssertPointer(func, 3);
3743:   PetscAssertPointer(line, 4);
3744: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3745:   {
3746:     const int index = x->lockstack.currentsize - 1;

3748:     *file = index < 0 ? NULL : x->lockstack.file[index];
3749:     *func = index < 0 ? NULL : x->lockstack.function[index];
3750:     *line = index < 0 ? 0 : x->lockstack.line[index];
3751:   }
3752: #else
3753:   *file = NULL;
3754:   *func = NULL;
3755:   *line = 0;
3756: #endif
3757:   PetscFunctionReturn(PETSC_SUCCESS);
3758: }

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

3763:   Logically Collective

3765:   Input Parameter:
3766: . x - the vector

3768:   Level: intermediate

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

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

3776: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3777: @*/
3778: PetscErrorCode VecLockReadPush(Vec x)
3779: {
3780:   PetscFunctionBegin;
3782:   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");
3783: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3784:   {
3785:     const char *file, *func;
3786:     int         index, line;

3788:     if ((index = petscstack.currentsize - 2) < 0) {
3789:       // vec was locked "outside" of petsc, either in user-land or main. the error message will
3790:       // now show this function as the culprit, but it will include the stacktrace
3791:       file = "unknown user-file";
3792:       func = "unknown_user_function";
3793:       line = 0;
3794:     } else {
3795:       file = petscstack.file[index];
3796:       func = petscstack.function[index];
3797:       line = petscstack.line[index];
3798:     }
3799:     PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
3800:   }
3801: #endif
3802:   PetscFunctionReturn(PETSC_SUCCESS);
3803: }

3805: /*@
3806:   VecLockReadPop - Pop a read-only lock from a vector

3808:   Logically Collective

3810:   Input Parameter:
3811: . x - the vector

3813:   Level: intermediate

3815: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
3816: @*/
3817: PetscErrorCode VecLockReadPop(Vec x)
3818: {
3819:   PetscFunctionBegin;
3821:   PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
3822: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3823:   {
3824:     const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];

3826:     PetscStackPop_Private(x->lockstack, previous);
3827:   }
3828: #endif
3829:   PetscFunctionReturn(PETSC_SUCCESS);
3830: }

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

3835:   Logically Collective

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

3841:   Level: intermediate

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

3849: .vb
3850:        VecGetArray(x,&xdata); // begin phase
3851:        VecLockWriteSet(v,PETSC_TRUE);

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

3855:        VecRestoreArray(x,&vdata); // end phase
3856:        VecLockWriteSet(v,PETSC_FALSE);
3857: .ve

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

3862: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
3863: @*/
3864: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
3865: {
3866:   PetscFunctionBegin;
3868:   if (flg) {
3869:     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");
3870:     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");
3871:     x->lock = -1;
3872:   } else {
3873:     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");
3874:     x->lock = 0;
3875:   }
3876:   PetscFunctionReturn(PETSC_SUCCESS);
3877: }