Actual source code: rvector.c

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

  9: PetscInt VecGetSubVectorSavedStateId = -1;

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

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

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

 41:    Logically Collective

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

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

 50:    Level: advanced

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

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

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

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

 82:    Collective

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

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

 91:    Performance Issues:
 92: .vb
 93:     per-processor memory bandwidth
 94:     interprocessor latency
 95:     work load imbalance that causes certain processes to arrive much earlier than others
 96: .ve

 98:    Level: intermediate

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

107:    Use `VecTDot()` for the indefinite form
108: $     val = (x,y) = y^T x,
109:    where y^T denotes the transpose of y.

111: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDotRealPart()`
112: @*/
113: PetscErrorCode VecDot(Vec x, Vec y, PetscScalar *val)
114: {
115:   PetscFunctionBegin;
121:   PetscCheckSameTypeAndComm(x, 1, y, 2);
122:   VecCheckSameSize(x, 1, y, 2);
123:   VecCheckAssembled(x);
124:   VecCheckAssembled(y);

126:   PetscCall(VecLockReadPush(x));
127:   PetscCall(VecLockReadPush(y));
128:   PetscCall(PetscLogEventBegin(VEC_Dot, x, y, 0, 0));
129:   PetscUseTypeMethod(x, dot, y, val);
130:   PetscCall(PetscLogEventEnd(VEC_Dot, x, y, 0, 0));
131:   PetscCall(VecLockReadPop(x));
132:   PetscCall(VecLockReadPop(y));
133:   PetscFunctionReturn(PETSC_SUCCESS);
134: }

136: /*@
137:    VecDotRealPart - Computes the real part of the vector dot product.

139:    Collective

141:    Input Parameters:
142: +  x - first vector
143: -  y - second vector

145:    Output Parameter:
146: .  val - the real part of the dot product;

148:    Level: intermediate

150:    Performance Issues:
151: .vb
152:     per-processor memory bandwidth
153:     interprocessor latency
154:     work load imbalance that causes certain processes to arrive much earlier than others
155: .ve

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

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

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

165:    Developer Note:
166:     This is not currently optimized to compute only the real part of the dot product.

168: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDot()`, `VecDotNorm2()`
169: @*/
170: PetscErrorCode VecDotRealPart(Vec x, Vec y, PetscReal *val)
171: {
172:   PetscScalar fdot;

174:   PetscFunctionBegin;
175:   PetscCall(VecDot(x, y, &fdot));
176:   *val = PetscRealPart(fdot);
177:   PetscFunctionReturn(PETSC_SUCCESS);
178: }

180: /*@
181:    VecNorm  - Computes the vector norm.

183:    Collective

185:    Input Parameters:
186: +  x - the vector
187: -  type - the type of the norm requested

189:    Output Parameter:
190: .  val - the norm

192:    Values of NormType:
193: +     `NORM_1` - sum_i |x[i]|
194: .     `NORM_2` - sqrt(sum_i |x[i]|^2)
195: .     `NORM_INFINITY` - max_i |x[i]|
196: -     `NORM_1_AND_2` - computes efficiently both  `NORM_1` and `NORM_2` and stores them each in an output array

198:     Level: intermediate

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

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

211:    Performance Issues:
212: +    per-processor memory bandwidth - limits the speed of the computation of local portion of the norm
213: .    interprocessor latency - limits the accumulation of the result across ranks, .i.e. MPI_Allreduce() time
214: .    number of ranks - the time for the result will grow with the log base 2 of the number of ranks sharing the vector
215: -    work load imbalance - the rank with the largest number of vector entries will limit the speed up

217: .seealso: [](ch_vectors), `Vec`, `NormType`, `VecDot()`, `VecTDot()`, `VecDotBegin()`, `VecDotEnd()`, `VecNormAvailable()`,
218:           `VecNormBegin()`, `VecNormEnd()`, `NormType()`
219: @*/
220: PetscErrorCode VecNorm(Vec x, NormType type, PetscReal *val)
221: {
222:   PetscBool flg = PETSC_TRUE;

224:   PetscFunctionBegin;
227:   VecCheckAssembled(x);

231:   /* Cached data? */
232:   PetscCall(VecNormAvailable(x, type, &flg, val));
233:   if (flg) PetscFunctionReturn(PETSC_SUCCESS);

235:   PetscCall(VecLockReadPush(x));
236:   PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
237:   PetscUseTypeMethod(x, norm, type, val);
238:   PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
239:   PetscCall(VecLockReadPop(x));

241:   if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
242:   PetscFunctionReturn(PETSC_SUCCESS);
243: }

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

248:    Not Collective

250:    Input Parameters:
251: +  x - the vector
252: -  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
253:           `NORM_1_AND_2`, which computes both norms and stores them
254:           in a two element array.

256:    Output Parameters:
257: +  available - `PETSC_TRUE` if the val returned is valid
258: -  val - the norm

260:    Level: intermediate

262:    Performance Issues:
263: .vb
264:     per-processor memory bandwidth
265:     interprocessor latency
266:     work load imbalance that causes certain processes to arrive much earlier than others
267: .ve

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

274: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
275:           `VecNormBegin()`, `VecNormEnd()`
276: @*/
277: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
278: {
279:   PetscFunctionBegin;

285:   if (type == NORM_1_AND_2) {
286:     *available = PETSC_FALSE;
287:   } else {
288:     PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
289:   }
290:   PetscFunctionReturn(PETSC_SUCCESS);
291: }

293: /*@
294:    VecNormalize - Normalizes a vector by its 2-norm.

296:    Collective

298:    Input Parameter:
299: .  x - the vector

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

304:    Level: intermediate

306: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
307: @*/
308: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
309: {
310:   PetscReal norm;

312:   PetscFunctionBegin;
315:   PetscCall(VecSetErrorIfLocked(x, 1));
317:   PetscCall(PetscLogEventBegin(VEC_Normalize, x, 0, 0, 0));
318:   PetscCall(VecNorm(x, NORM_2, &norm));
319:   if (norm == 0.0) PetscCall(PetscInfo(x, "Vector of zero norm can not be normalized; Returning only the zero norm\n"));
320:   else if (PetscIsInfOrNanReal(norm)) PetscCall(PetscInfo(x, "Vector with Inf or Nan norm can not be normalized; Returning only the norm\n"));
321:   else {
322:     PetscScalar s = 1.0 / norm;
323:     PetscCall(VecScale(x, s));
324:   }
325:   PetscCall(PetscLogEventEnd(VEC_Normalize, x, 0, 0, 0));
326:   if (val) *val = norm;
327:   PetscFunctionReturn(PETSC_SUCCESS);
328: }

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

333:    Collective

335:    Input Parameter:
336: .  x - the vector

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

342:    Level: intermediate

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

347:    Returns the smallest index with the maximum value

349: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
350: @*/
351: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
352: {
353:   PetscFunctionBegin;
356:   VecCheckAssembled(x);
359:   PetscCall(VecLockReadPush(x));
360:   PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
361:   PetscUseTypeMethod(x, max, p, val);
362:   PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
363:   PetscCall(VecLockReadPop(x));
364:   PetscFunctionReturn(PETSC_SUCCESS);
365: }

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

370:    Collective

372:    Input Parameter:
373: .  x - the vector

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

379:    Level: intermediate

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

384:    This returns the smallest index with the minimum value

386: .seealso: [](ch_vectors), `Vec`, `VecMax()`
387: @*/
388: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
389: {
390:   PetscFunctionBegin;
393:   VecCheckAssembled(x);
396:   PetscCall(VecLockReadPush(x));
397:   PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
398:   PetscUseTypeMethod(x, min, p, val);
399:   PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
400:   PetscCall(VecLockReadPop(x));
401:   PetscFunctionReturn(PETSC_SUCCESS);
402: }

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

408:    Collective

410:    Input Parameters:
411: +  x - first vector
412: -  y - second vector

414:    Output Parameter:
415: .  val - the dot product

417:    Level: intermediate

419:    Notes for Users of Complex Numbers:
420:    For complex vectors, `VecTDot()` computes the indefinite form
421: $     val = (x,y) = y^T x,
422:    where y^T denotes the transpose of y.

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

428: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
429: @*/
430: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
431: {
432:   PetscFunctionBegin;
438:   PetscCheckSameTypeAndComm(x, 1, y, 2);
439:   VecCheckSameSize(x, 1, y, 2);
440:   VecCheckAssembled(x);
441:   VecCheckAssembled(y);

443:   PetscCall(VecLockReadPush(x));
444:   PetscCall(VecLockReadPush(y));
445:   PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
446:   PetscUseTypeMethod(x, tdot, y, val);
447:   PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
448:   PetscCall(VecLockReadPop(x));
449:   PetscCall(VecLockReadPop(y));
450:   PetscFunctionReturn(PETSC_SUCCESS);
451: }

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

456:    Not Collective

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

462:    Level: intermediate

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

467: .seealso: [](ch_vectors), `Vec`, `VecSet()`
468: @*/
469: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
470: {
471:   PetscReal   norms[4];
472:   PetscBool   flgs[4];
473:   PetscScalar one = 1.0;

475:   PetscFunctionBegin;
478:   VecCheckAssembled(x);
479:   PetscCall(VecSetErrorIfLocked(x, 1));
480:   if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);

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

485:   PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
486:   PetscUseTypeMethod(x, scale, alpha);
487:   PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));

489:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
490:   /* put the scaled stashed norms back into the Vec */
491:   for (PetscInt i = 0; i < 4; i++) {
492:     PetscReal ar = PetscAbsScalar(alpha);
493:     if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
494:   }
495:   PetscFunctionReturn(PETSC_SUCCESS);
496: }

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

501:    Logically Collective

503:    Input Parameters:
504: +  x  - the vector
505: -  alpha - the scalar

507:    Level: beginner

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

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

518: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
519: @*/
520: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
521: {
522:   PetscFunctionBegin;
525:   VecCheckAssembled(x);
527:   PetscCall(VecSetErrorIfLocked(x, 1));

529:   PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
530:   PetscUseTypeMethod(x, set, alpha);
531:   PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
532:   PetscCall(PetscObjectStateIncrease((PetscObject)x));

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

536:   {
537:     PetscReal      val = PetscAbsScalar(alpha);
538:     const PetscInt N   = x->map->N;

540:     if (N == 0) {
541:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0l));
542:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
543:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
544:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
545:     } else if (val > PETSC_MAX_REAL / N) {
546:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
547:     } else {
548:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
549:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
550:       val *= PetscSqrtReal((PetscReal)N);
551:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
552:       PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
553:     }
554:   }
555:   PetscFunctionReturn(PETSC_SUCCESS);
556: }

558: /*@
559:    VecAXPY - Computes `y = alpha x + y`.

561:    Logically Collective

563:    Input Parameters:
564: +  alpha - the scalar
565: .  x - vector scale by `alpha`
566: -  y - vector accumulated into

568:    Output Parameter:
569: .  y - output vector

571:    Level: intermediate

573:    Notes:
574:     This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
575: .vb
576:     VecAXPY(y,alpha,x)                   y = alpha x           +      y
577:     VecAYPX(y,beta,x)                    y =       x           + beta y
578:     VecAXPBY(y,alpha,beta,x)             y = alpha x           + beta y
579:     VecWAXPY(w,alpha,x,y)                w = alpha x           +      y
580:     VecAXPBYPCZ(w,alpha,beta,gamma,x,y)  z = alpha x           + beta y + gamma z
581:     VecMAXPY(y,nv,alpha[],x[])           y = sum alpha[i] x[i] +      y
582: .ve

584: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
585: @*/
586: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
587: {
588:   PetscFunctionBegin;
593:   PetscCheckSameTypeAndComm(x, 3, y, 1);
594:   VecCheckSameSize(x, 3, y, 1);
595:   VecCheckAssembled(x);
596:   VecCheckAssembled(y);
598:   if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
599:   PetscCall(VecSetErrorIfLocked(y, 1));
600:   if (x == y) {
601:     PetscCall(VecScale(y, alpha + 1.0));
602:     PetscFunctionReturn(PETSC_SUCCESS);
603:   }
604:   PetscCall(VecLockReadPush(x));
605:   PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
606:   PetscUseTypeMethod(y, axpy, alpha, x);
607:   PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
608:   PetscCall(VecLockReadPop(x));
609:   PetscCall(PetscObjectStateIncrease((PetscObject)y));
610:   PetscFunctionReturn(PETSC_SUCCESS);
611: }

613: /*@
614:    VecAYPX - Computes `y = x + beta y`.

616:    Logically Collective

618:    Input Parameters:
619: +  beta - the scalar
620: .  x - the unscaled vector
621: -  y - the vector to be scaled

623:    Output Parameter:
624: .  y - output vector

626:    Level: intermediate

628:    Developer Note:
629:     The implementation is optimized for `beta` of -1.0, 0.0, and 1.0

631: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
632: @*/
633: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
634: {
635:   PetscFunctionBegin;
640:   PetscCheckSameTypeAndComm(x, 3, y, 1);
641:   VecCheckSameSize(x, 1, y, 3);
642:   VecCheckAssembled(x);
643:   VecCheckAssembled(y);
645:   PetscCall(VecSetErrorIfLocked(y, 1));
646:   if (x == y) {
647:     PetscCall(VecScale(y, beta + 1.0));
648:     PetscFunctionReturn(PETSC_SUCCESS);
649:   }
650:   PetscCall(VecLockReadPush(x));
651:   if (beta == (PetscScalar)0.0) {
652:     PetscCall(VecCopy(x, y));
653:   } else {
654:     PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
655:     PetscUseTypeMethod(y, aypx, beta, x);
656:     PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
657:     PetscCall(PetscObjectStateIncrease((PetscObject)y));
658:   }
659:   PetscCall(VecLockReadPop(x));
660:   PetscFunctionReturn(PETSC_SUCCESS);
661: }

663: /*@
664:    VecAXPBY - Computes `y = alpha x + beta y`.

666:    Logically Collective

668:    Input Parameters:
669: +  alpha - first scalar
670: .  beta - second scalar
671: .  x - the first scaled vector
672: -  y - the second scaled vector

674:    Output Parameter:
675: .  y - output vector

677:    Level: intermediate

679:    Developer Note:
680:    The implementation is optimized for `alpha` and/or `beta` values of 0.0 and 1.0

682: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
683: @*/
684: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
685: {
686:   PetscFunctionBegin;
691:   PetscCheckSameTypeAndComm(x, 4, y, 1);
692:   VecCheckSameSize(y, 1, x, 4);
693:   VecCheckAssembled(x);
694:   VecCheckAssembled(y);
697:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
698:   if (x == y) {
699:     PetscCall(VecScale(y, alpha + beta));
700:     PetscFunctionReturn(PETSC_SUCCESS);
701:   }

703:   PetscCall(VecSetErrorIfLocked(y, 1));
704:   PetscCall(VecLockReadPush(x));
705:   PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
706:   PetscUseTypeMethod(y, axpby, alpha, beta, x);
707:   PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
708:   PetscCall(PetscObjectStateIncrease((PetscObject)y));
709:   PetscCall(VecLockReadPop(x));
710:   PetscFunctionReturn(PETSC_SUCCESS);
711: }

713: /*@
714:    VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`

716:    Logically Collective

718:    Input Parameters:
719: +  alpha - first scalar
720: .  beta - second scalar
721: .  gamma - third scalar
722: .  x  - first vector
723: .  y  - second vector
724: -  z  - third vector

726:    Output Parameter:
727: .  z - output vector

729:    Level: intermediate

731:    Note:
732:    `x`, `y` and `z` must be different vectors

734:    Developer Note:
735:     The implementation is optimized for `alpha` of 1.0 and `gamma` of 1.0 or 0.0

737: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
738: @*/
739: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
740: {
741:   PetscFunctionBegin;
748:   PetscCheckSameTypeAndComm(x, 5, y, 6);
749:   PetscCheckSameTypeAndComm(x, 5, z, 1);
750:   VecCheckSameSize(x, 1, y, 5);
751:   VecCheckSameSize(x, 1, z, 6);
752:   PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
753:   PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
754:   VecCheckAssembled(x);
755:   VecCheckAssembled(y);
756:   VecCheckAssembled(z);
760:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);

762:   PetscCall(VecSetErrorIfLocked(z, 1));
763:   PetscCall(VecLockReadPush(x));
764:   PetscCall(VecLockReadPush(y));
765:   PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
766:   PetscUseTypeMethod(z, axpbypcz, alpha, beta, gamma, x, y);
767:   PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
768:   PetscCall(PetscObjectStateIncrease((PetscObject)z));
769:   PetscCall(VecLockReadPop(x));
770:   PetscCall(VecLockReadPop(y));
771:   PetscFunctionReturn(PETSC_SUCCESS);
772: }

774: /*@
775:    VecWAXPY - Computes `w = alpha x + y`.

777:    Logically Collective

779:    Input Parameters:
780: +  alpha - the scalar
781: .  x  - first vector, multiplied by `alpha`
782: -  y  - second vector

784:    Output Parameter:
785: .  w - the result

787:    Level: intermediate

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

792:    Developer Note:
793:     The implementation is optimized for alpha of -1.0, 0.0, and 1.0

795: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
796: @*/
797: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
798: {
799:   PetscFunctionBegin;
806:   PetscCheckSameTypeAndComm(x, 3, y, 4);
807:   PetscCheckSameTypeAndComm(y, 4, w, 1);
808:   VecCheckSameSize(x, 3, y, 4);
809:   VecCheckSameSize(x, 3, w, 1);
810:   PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
811:   PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
812:   VecCheckAssembled(x);
813:   VecCheckAssembled(y);
815:   PetscCall(VecSetErrorIfLocked(w, 1));

817:   PetscCall(VecLockReadPush(x));
818:   PetscCall(VecLockReadPush(y));
819:   if (alpha == (PetscScalar)0.0) {
820:     PetscCall(VecCopy(y, w));
821:   } else {
822:     PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
823:     PetscUseTypeMethod(w, waxpy, alpha, x, y);
824:     PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
825:     PetscCall(PetscObjectStateIncrease((PetscObject)w));
826:   }
827:   PetscCall(VecLockReadPop(x));
828:   PetscCall(VecLockReadPop(y));
829:   PetscFunctionReturn(PETSC_SUCCESS);
830: }

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

835:    Not Collective

837:    Input Parameters:
838: +  x - vector to insert in
839: .  ni - number of elements to add
840: .  ix - indices where to add
841: .  y - array of values
842: -  iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries

844:    Level: beginner

846:    Notes:
847: .vb
848:    `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
849: .ve

851:    Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
852:    options cannot be mixed without intervening calls to the assembly
853:    routines.

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

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

860:    If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
861:    negative indices may be passed in ix. These rows are
862:    simply ignored. This allows easily inserting element load matrices
863:    with homogeneous Dirchlet boundary conditions that you don't want represented
864:    in the vector.

866: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
867:           `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
868: @*/
869: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
870: {
871:   PetscFunctionBeginHot;
873:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);

878:   PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
879:   PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
880:   PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
881:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
882:   PetscFunctionReturn(PETSC_SUCCESS);
883: }

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

889:     Not Collective

891:    Input Parameters:
892: +  x - vector to get values from
893: .  ni - number of elements to get
894: -  ix - indices where to get them from (in global 1d numbering)

896:    Output Parameter:
897: .   y - array of values

899:    Level: beginner

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

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

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

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

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

914: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
915: @*/
916: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
917: {
918:   PetscFunctionBegin;
920:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
924:   VecCheckAssembled(x);
925:   PetscUseTypeMethod(x, getvalues, ni, ix, y);
926:   PetscFunctionReturn(PETSC_SUCCESS);
927: }

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

932:    Not Collective

934:    Input Parameters:
935: +  x - vector to insert in
936: .  ni - number of blocks to add
937: .  ix - indices where to add in block count, rather than element count
938: .  y - array of values
939: -  iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries

941:    Level: intermediate

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

947:    Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
948:    options cannot be mixed without intervening calls to the assembly
949:    routines.

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

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

956:    Negative indices may be passed in ix, these rows are
957:    simply ignored. This allows easily inserting element load matrices
958:    with homogeneous Dirchlet boundary conditions that you don't want represented
959:    in the vector.

961: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
962:           `VecSetValues()`
963: @*/
964: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
965: {
966:   PetscFunctionBeginHot;
968:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);

973:   PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
974:   PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
975:   PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
976:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
977:   PetscFunctionReturn(PETSC_SUCCESS);
978: }

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

984:    Not Collective

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

993:    Level: intermediate

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

998:    Calls to `VecSetValues()` 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 `VecSetValuesLocal()` have been completed.

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

1007: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1008:           `VecSetValuesBlockedLocal()`
1009: @*/
1010: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1011: {
1012:   PetscInt lixp[128], *lix = lixp;

1014:   PetscFunctionBeginHot;
1016:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);

1021:   PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1022:   if (!x->ops->setvalueslocal) {
1023:     if (x->map->mapping) {
1024:       if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1025:       PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1026:       PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1027:       if (ni > 128) PetscCall(PetscFree(lix));
1028:     } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1029:   } else PetscUseTypeMethod(x, setvalueslocal, ni, ix, y, iora);
1030:   PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1031:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
1032:   PetscFunctionReturn(PETSC_SUCCESS);
1033: }

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

1039:    Not Collective

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

1048:    Level: intermediate

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

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

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

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

1063: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1064:           `VecSetLocalToGlobalMapping()`
1065: @*/
1066: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1067: {
1068:   PetscInt lixp[128], *lix = lixp;

1070:   PetscFunctionBeginHot;
1072:   if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1076:   PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1077:   if (x->map->mapping) {
1078:     if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1079:     PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1080:     PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1081:     if (ni > 128) PetscCall(PetscFree(lix));
1082:   } else {
1083:     PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1084:   }
1085:   PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1086:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
1087:   PetscFunctionReturn(PETSC_SUCCESS);
1088: }

1090: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1091: {
1092:   PetscFunctionBegin;
1095:   VecCheckAssembled(x);
1097:   if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1099:   for (PetscInt i = 0; i < nv; ++i) {
1102:     PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1103:     VecCheckSameSize(x, 1, y[i], 3);
1104:     VecCheckAssembled(y[i]);
1105:     PetscCall(VecLockReadPush(y[i]));
1106:   }

1110:   PetscCall(VecLockReadPush(x));
1111:   PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1112:   PetscCall((*mxdot)(x, nv, y, result));
1113:   PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1114:   PetscCall(VecLockReadPop(x));
1115:   for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1116:   PetscFunctionReturn(PETSC_SUCCESS);
1117: }

1119: /*@
1120:    VecMTDot - Computes indefinite vector multiple dot products.
1121:    That is, it does NOT use the complex conjugate.

1123:    Collective

1125:    Input Parameters:
1126: +  x - one vector
1127: .  nv - number of vectors
1128: -  y - array of vectors.  Note that vectors are pointers

1130:    Output Parameter:
1131: .  val - array of the dot products

1133:    Level: intermediate

1135:    Notes for Users of Complex Numbers:
1136:    For complex vectors, `VecMTDot()` computes the indefinite form
1137: $      val = (x,y) = y^T x,
1138:    where y^T denotes the transpose of y.

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

1144: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1145: @*/
1146: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1147: {
1148:   PetscFunctionBegin;
1150:   PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1151:   PetscFunctionReturn(PETSC_SUCCESS);
1152: }

1154: /*@
1155:    VecMDot - Computes multiple vector dot products.

1157:    Collective

1159:    Input Parameters:
1160: +  x - one vector
1161: .  nv - number of vectors
1162: -  y - array of vectors.

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

1167:    Level: intermediate

1169:    Notes for Users of Complex Numbers:
1170:    For complex vectors, `VecMDot()` computes
1171: $     val = (x,y) = y^H x,
1172:    where y^H denotes the conjugate transpose of y.

1174:    Use `VecMTDot()` for the indefinite form
1175: $     val = (x,y) = y^T x,
1176:    where y^T denotes the transpose of y.

1178: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`
1179: @*/
1180: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1181: {
1182:   PetscFunctionBegin;
1184:   PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1185:   PetscFunctionReturn(PETSC_SUCCESS);
1186: }

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

1191:    Logically Collective

1193:    Input Parameters:
1194: +  nv - number of scalars and x-vectors
1195: .  alpha - array of scalars
1196: .  y - one vector
1197: -  x - array of vectors

1199:    Level: intermediate

1201:    Note:
1202:     `y` cannot be any of the `x` vectors

1204: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1205: @*/
1206: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1207: {
1208:   PetscFunctionBegin;
1210:   VecCheckAssembled(y);
1212:   PetscCall(VecSetErrorIfLocked(y, 1));
1213:   PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1214:   if (nv) {
1215:     PetscInt zeros = 0;

1219:     for (PetscInt i = 0; i < nv; ++i) {
1223:       PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1224:       VecCheckSameSize(y, 1, x[i], 4);
1225:       PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1226:       VecCheckAssembled(x[i]);
1227:       PetscCall(VecLockReadPush(x[i]));
1228:       zeros += alpha[i] == (PetscScalar)0.0;
1229:     }

1231:     if (zeros < nv) {
1232:       PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1233:       PetscUseTypeMethod(y, maxpy, nv, alpha, x);
1234:       PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1235:       PetscCall(PetscObjectStateIncrease((PetscObject)y));
1236:     }

1238:     for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1239:   }
1240:   PetscFunctionReturn(PETSC_SUCCESS);
1241: }

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

1248:    Collective

1250:    Input Parameters:
1251: +  nx   - number of vectors to be concatenated
1252: -  X    - array containing the vectors to be concatenated in the order of concatenation

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

1258:    Level: advanced

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

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

1270: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1271: @*/
1272: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1273: {
1274:   MPI_Comm comm;
1275:   VecType  vec_type;
1276:   Vec      Ytmp, Xtmp;
1277:   IS      *is_tmp;
1278:   PetscInt i, shift = 0, Xnl, Xng, Xbegin;

1280:   PetscFunctionBegin;

1286:   if ((*X)->ops->concatenate) {
1287:     /* use the dedicated concatenation function if available */
1288:     PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1289:   } else {
1290:     /* loop over vectors and start creating IS */
1291:     comm = PetscObjectComm((PetscObject)(*X));
1292:     PetscCall(VecGetType(*X, &vec_type));
1293:     PetscCall(PetscMalloc1(nx, &is_tmp));
1294:     for (i = 0; i < nx; i++) {
1295:       PetscCall(VecGetSize(X[i], &Xng));
1296:       PetscCall(VecGetLocalSize(X[i], &Xnl));
1297:       PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1298:       PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1299:       shift += Xng;
1300:     }
1301:     /* create the concatenated vector */
1302:     PetscCall(VecCreate(comm, &Ytmp));
1303:     PetscCall(VecSetType(Ytmp, vec_type));
1304:     PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1305:     PetscCall(VecSetUp(Ytmp));
1306:     /* copy data from X array to Y and return */
1307:     for (i = 0; i < nx; i++) {
1308:       PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1309:       PetscCall(VecCopy(X[i], Xtmp));
1310:       PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1311:     }
1312:     *Y = Ytmp;
1313:     if (x_is) {
1314:       *x_is = is_tmp;
1315:     } else {
1316:       for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1317:       PetscCall(PetscFree(is_tmp));
1318:     }
1319:   }
1320:   PetscFunctionReturn(PETSC_SUCCESS);
1321: }

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

1326:     Input Parameters:
1327: +   X - the original vector
1328: -   is - the index set of the subvector

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

1335: */
1336: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1337: {
1338:   PetscInt  gstart, gend, lstart;
1339:   PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1340:   PetscInt  n, N, ibs, vbs, bs = -1;

1342:   PetscFunctionBegin;
1343:   PetscCall(ISGetLocalSize(is, &n));
1344:   PetscCall(ISGetSize(is, &N));
1345:   PetscCall(ISGetBlockSize(is, &ibs));
1346:   PetscCall(VecGetBlockSize(X, &vbs));
1347:   PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1348:   PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1349:   /* block size is given by IS if ibs > 1; otherwise, check the vector */
1350:   if (ibs > 1) {
1351:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1352:     bs = ibs;
1353:   } else {
1354:     if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1355:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1356:     if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1357:   }

1359:   *contig    = red[0];
1360:   *start     = lstart;
1361:   *blocksize = bs;
1362:   PetscFunctionReturn(PETSC_SUCCESS);
1363: }

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

1367:     Input Parameters:
1368: +   X - the original vector
1369: .   is - the index set of the subvector
1370: -   bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()

1372:     Output Parameter:
1373: .   Z  - the subvector, which will compose the VecScatter context on output
1374: */
1375: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1376: {
1377:   PetscInt   n, N;
1378:   VecScatter vscat;
1379:   Vec        Y;

1381:   PetscFunctionBegin;
1382:   PetscCall(ISGetLocalSize(is, &n));
1383:   PetscCall(ISGetSize(is, &N));
1384:   PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1385:   PetscCall(VecSetSizes(Y, n, N));
1386:   PetscCall(VecSetBlockSize(Y, bs));
1387:   PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1388:   PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1389:   PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1390:   PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1391:   PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1392:   PetscCall(VecScatterDestroy(&vscat));
1393:   *Z = Y;
1394:   PetscFunctionReturn(PETSC_SUCCESS);
1395: }

1397: /*@
1398:    VecGetSubVector - Gets a vector representing part of another vector

1400:    Collective

1402:    Input Parameters:
1403: +  X - vector from which to extract a subvector
1404: -  is - index set representing portion of `X` to extract

1406:    Output Parameter:
1407: .  Y - subvector corresponding to `is`

1409:    Level: advanced

1411:    Notes:
1412:    The subvector `Y` should be returned with `VecRestoreSubVector()`.
1413:    `X` and must be defined on the same communicator

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

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

1420: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1421: @*/
1422: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1423: {
1424:   Vec Z;

1426:   PetscFunctionBegin;
1429:   PetscCheckSameComm(X, 1, is, 2);
1431:   if (X->ops->getsubvector) {
1432:     PetscUseTypeMethod(X, getsubvector, is, &Z);
1433:   } else { /* Default implementation currently does no caching */
1434:     PetscBool contig;
1435:     PetscInt  n, N, start, bs;

1437:     PetscCall(ISGetLocalSize(is, &n));
1438:     PetscCall(ISGetSize(is, &N));
1439:     PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1440:     if (contig) { /* We can do a no-copy implementation */
1441:       const PetscScalar *x;
1442:       PetscInt           state = 0;
1443:       PetscBool          isstd, iscuda, iship;

1445:       PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1446:       PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1447:       PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1448:       if (iscuda) {
1449: #if defined(PETSC_HAVE_CUDA)
1450:         const PetscScalar *x_d;
1451:         PetscMPIInt        size;
1452:         PetscOffloadMask   flg;

1454:         PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1455:         PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1456:         PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1457:         if (x) x += start;
1458:         if (x_d) x_d += start;
1459:         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1460:         if (size == 1) {
1461:           PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1462:         } else {
1463:           PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1464:         }
1465:         Z->offloadmask = flg;
1466: #endif
1467:       } else if (iship) {
1468: #if defined(PETSC_HAVE_HIP)
1469:         const PetscScalar *x_d;
1470:         PetscMPIInt        size;
1471:         PetscOffloadMask   flg;

1473:         PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1474:         PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1475:         PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1476:         if (x) x += start;
1477:         if (x_d) x_d += start;
1478:         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1479:         if (size == 1) {
1480:           PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1481:         } else {
1482:           PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1483:         }
1484:         Z->offloadmask = flg;
1485: #endif
1486:       } else if (isstd) {
1487:         PetscMPIInt size;

1489:         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1490:         PetscCall(VecGetArrayRead(X, &x));
1491:         if (x) x += start;
1492:         if (size == 1) {
1493:           PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1494:         } else {
1495:           PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1496:         }
1497:         PetscCall(VecRestoreArrayRead(X, &x));
1498:       } else { /* default implementation: use place array */
1499:         PetscCall(VecGetArrayRead(X, &x));
1500:         PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1501:         PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1502:         PetscCall(VecSetSizes(Z, n, N));
1503:         PetscCall(VecSetBlockSize(Z, bs));
1504:         PetscCall(VecPlaceArray(Z, x ? x + start : NULL));
1505:         PetscCall(VecRestoreArrayRead(X, &x));
1506:       }

1508:       /* this is relevant only in debug mode */
1509:       PetscCall(VecLockGet(X, &state));
1510:       if (state) PetscCall(VecLockReadPush(Z));
1511:       Z->ops->placearray   = NULL;
1512:       Z->ops->replacearray = NULL;
1513:     } else { /* Have to create a scatter and do a copy */
1514:       PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1515:     }
1516:   }
1517:   /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1518:   if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1519:   PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1520:   *Y = Z;
1521:   PetscFunctionReturn(PETSC_SUCCESS);
1522: }

1524: /*@
1525:    VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`

1527:    Collective

1529:    Input Parameters:
1530: + X - vector from which subvector was obtained
1531: . is - index set representing the subset of `X`
1532: - Y - subvector being restored

1534:    Level: advanced

1536: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1537: @*/
1538: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1539: {
1540:   PETSC_UNUSED PetscObjectState dummystate = 0;
1541:   PetscBool                     unchanged;

1543:   PetscFunctionBegin;
1546:   PetscCheckSameComm(X, 1, is, 2);

1550:   if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1551:   else {
1552:     PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1553:     if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1554:       VecScatter scatter;
1555:       PetscInt   state;

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

1560:       PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1561:       if (scatter) {
1562:         PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1563:         PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1564:       } else {
1565:         PetscBool iscuda, iship;
1566:         PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1567:         PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));

1569:         if (iscuda) {
1570: #if defined(PETSC_HAVE_CUDA)
1571:           PetscOffloadMask ymask = (*Y)->offloadmask;

1573:           /* The offloadmask of X dictates where to move memory
1574:               If X GPU data is valid, then move Y data on GPU if needed
1575:               Otherwise, move back to the CPU */
1576:           switch (X->offloadmask) {
1577:           case PETSC_OFFLOAD_BOTH:
1578:             if (ymask == PETSC_OFFLOAD_CPU) {
1579:               PetscCall(VecCUDAResetArray(*Y));
1580:             } else if (ymask == PETSC_OFFLOAD_GPU) {
1581:               X->offloadmask = PETSC_OFFLOAD_GPU;
1582:             }
1583:             break;
1584:           case PETSC_OFFLOAD_GPU:
1585:             if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1586:             break;
1587:           case PETSC_OFFLOAD_CPU:
1588:             if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1589:             break;
1590:           case PETSC_OFFLOAD_UNALLOCATED:
1591:           case PETSC_OFFLOAD_KOKKOS:
1592:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1593:           }
1594: #endif
1595:         } else if (iship) {
1596: #if defined(PETSC_HAVE_HIP)
1597:           PetscOffloadMask ymask = (*Y)->offloadmask;

1599:           /* The offloadmask of X dictates where to move memory
1600:               If X GPU data is valid, then move Y data on GPU if needed
1601:               Otherwise, move back to the CPU */
1602:           switch (X->offloadmask) {
1603:           case PETSC_OFFLOAD_BOTH:
1604:             if (ymask == PETSC_OFFLOAD_CPU) {
1605:               PetscCall(VecHIPResetArray(*Y));
1606:             } else if (ymask == PETSC_OFFLOAD_GPU) {
1607:               X->offloadmask = PETSC_OFFLOAD_GPU;
1608:             }
1609:             break;
1610:           case PETSC_OFFLOAD_GPU:
1611:             if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1612:             break;
1613:           case PETSC_OFFLOAD_CPU:
1614:             if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1615:             break;
1616:           case PETSC_OFFLOAD_UNALLOCATED:
1617:           case PETSC_OFFLOAD_KOKKOS:
1618:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1619:           }
1620: #endif
1621:         } else {
1622:           /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1623:           PetscCall(VecResetArray(*Y));
1624:         }
1625:         PetscCall(PetscObjectStateIncrease((PetscObject)X));
1626:       }
1627:     }
1628:   }
1629:   PetscCall(VecDestroy(Y));
1630:   PetscFunctionReturn(PETSC_SUCCESS);
1631: }

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

1637:    Not Collective.

1639:    Input parameter:
1640: .  v - The vector for which the local vector is desired.

1642:    Output parameter:
1643: .  w - Upon exit this contains the local vector.

1645:    Level: beginner

1647: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1648: @*/
1649: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1650: {
1651:   PetscMPIInt size;

1653:   PetscFunctionBegin;
1656:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1657:   if (size == 1) PetscCall(VecDuplicate(v, w));
1658:   else if (v->ops->createlocalvector) PetscUseTypeMethod(v, createlocalvector, w);
1659:   else {
1660:     VecType  type;
1661:     PetscInt n;

1663:     PetscCall(VecCreate(PETSC_COMM_SELF, w));
1664:     PetscCall(VecGetLocalSize(v, &n));
1665:     PetscCall(VecSetSizes(*w, n, n));
1666:     PetscCall(VecGetBlockSize(v, &n));
1667:     PetscCall(VecSetBlockSize(*w, n));
1668:     PetscCall(VecGetType(v, &type));
1669:     PetscCall(VecSetType(*w, type));
1670:   }
1671:   PetscFunctionReturn(PETSC_SUCCESS);
1672: }

1674: /*@
1675:    VecGetLocalVectorRead - Maps the local portion of a vector into a
1676:    vector.

1678:    Not Collective.

1680:    Input parameter:
1681: .  v - The vector for which the local vector is desired.

1683:    Output parameter:
1684: .  w - Upon exit this contains the local vector.

1686:    Level: beginner

1688:    Notes:
1689:    You must call `VecRestoreLocalVectorRead()` when the local
1690:    vector is no longer needed.

1692:    This function is similar to `VecGetArrayRead()` which maps the local
1693:    portion into a raw pointer.  `VecGetLocalVectorRead()` is usually
1694:    almost as efficient as `VecGetArrayRead()` but in certain circumstances
1695:    `VecGetLocalVectorRead()` can be much more efficient than
1696:    `VecGetArrayRead()`.  This is because the construction of a contiguous
1697:    array representing the vector data required by `VecGetArrayRead()` can
1698:    be an expensive operation for certain vector types.  For example, for
1699:    GPU vectors `VecGetArrayRead()` requires that the data between device
1700:    and host is synchronized.

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

1705: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1706: @*/
1707: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1708: {
1709:   PetscFunctionBegin;
1712:   VecCheckSameLocalSize(v, 1, w, 2);
1713:   if (v->ops->getlocalvectorread) {
1714:     PetscUseTypeMethod(v, getlocalvectorread, w);
1715:   } else {
1716:     PetscScalar *a;

1718:     PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1719:     PetscCall(VecPlaceArray(w, a));
1720:   }
1721:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
1722:   PetscCall(VecLockReadPush(v));
1723:   PetscCall(VecLockReadPush(w));
1724:   PetscFunctionReturn(PETSC_SUCCESS);
1725: }

1727: /*@
1728:    VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1729:    previously mapped into a vector using `VecGetLocalVectorRead()`.

1731:    Not Collective.

1733:    Input parameter:
1734: +  v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
1735: -  w - The vector into which the local portion of `v` was mapped.

1737:    Level: beginner

1739: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1740: @*/
1741: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1742: {
1743:   PetscFunctionBegin;
1746:   if (v->ops->restorelocalvectorread) {
1747:     PetscUseTypeMethod(v, restorelocalvectorread, w);
1748:   } else {
1749:     const PetscScalar *a;

1751:     PetscCall(VecGetArrayRead(w, &a));
1752:     PetscCall(VecRestoreArrayRead(v, &a));
1753:     PetscCall(VecResetArray(w));
1754:   }
1755:   PetscCall(VecLockReadPop(v));
1756:   PetscCall(VecLockReadPop(w));
1757:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
1758:   PetscFunctionReturn(PETSC_SUCCESS);
1759: }

1761: /*@
1762:    VecGetLocalVector - Maps the local portion of a vector into a
1763:    vector.

1765:    Collective

1767:    Input parameter:
1768: .  v - The vector for which the local vector is desired.

1770:    Output parameter:
1771: .  w - Upon exit this contains the local vector.

1773:    Level: beginner

1775:    Notes:
1776:    You must call `VecRestoreLocalVector()` when the local
1777:    vector is no longer needed.

1779:    This function is similar to `VecGetArray()` which maps the local
1780:    portion into a raw pointer.  `VecGetLocalVector()` is usually about as
1781:    efficient as `VecGetArray()` but in certain circumstances
1782:    `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1783:    This is because the construction of a contiguous array representing
1784:    the vector data required by `VecGetArray()` can be an expensive
1785:    operation for certain vector types.  For example, for GPU vectors
1786:    `VecGetArray()` requires that the data between device and host is
1787:    synchronized.

1789: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1790: @*/
1791: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1792: {
1793:   PetscFunctionBegin;
1796:   VecCheckSameLocalSize(v, 1, w, 2);
1797:   if (v->ops->getlocalvector) {
1798:     PetscUseTypeMethod(v, getlocalvector, w);
1799:   } else {
1800:     PetscScalar *a;

1802:     PetscCall(VecGetArray(v, &a));
1803:     PetscCall(VecPlaceArray(w, a));
1804:   }
1805:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
1806:   PetscFunctionReturn(PETSC_SUCCESS);
1807: }

1809: /*@
1810:    VecRestoreLocalVector - Unmaps the local portion of a vector
1811:    previously mapped into a vector using `VecGetLocalVector()`.

1813:    Logically Collective.

1815:    Input parameter:
1816: +  v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
1817: -  w - The vector into which the local portion of `v` was mapped.

1819:    Level: beginner

1821: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1822: @*/
1823: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1824: {
1825:   PetscFunctionBegin;
1828:   if (v->ops->restorelocalvector) {
1829:     PetscUseTypeMethod(v, restorelocalvector, w);
1830:   } else {
1831:     PetscScalar *a;
1832:     PetscCall(VecGetArray(w, &a));
1833:     PetscCall(VecRestoreArray(v, &a));
1834:     PetscCall(VecResetArray(w));
1835:   }
1836:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
1837:   PetscCall(PetscObjectStateIncrease((PetscObject)v));
1838:   PetscFunctionReturn(PETSC_SUCCESS);
1839: }

1841: /*@C
1842:    VecGetArray - Returns a pointer to a contiguous array that contains this
1843:    processor's portion of the vector data. For the standard PETSc
1844:    vectors, `VecGetArray()` returns a pointer to the local data array and
1845:    does not use any copies. If the underlying vector data is not stored
1846:    in a contiguous array this routine will copy the data to a contiguous
1847:    array and return a pointer to that. You MUST call `VecRestoreArray()`
1848:    when you no longer need access to the array.

1850:    Logically Collective

1852:    Input Parameter:
1853: .  x - the vector

1855:    Output Parameter:
1856: .  a - location to put pointer to the array

1858:    Level: beginner

1860:    Fortran Note:
1861:    `VecGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayF90()`

1863: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
1864:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
1865: @*/
1866: PetscErrorCode VecGetArray(Vec x, PetscScalar **a)
1867: {
1868:   PetscFunctionBegin;
1870:   PetscCall(VecSetErrorIfLocked(x, 1));
1871:   if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
1872:     PetscUseTypeMethod(x, getarray, a);
1873:   } else if (x->petscnative) { /* VECSTANDARD */
1874:     *a = *((PetscScalar **)x->data);
1875:   } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
1876:   PetscFunctionReturn(PETSC_SUCCESS);
1877: }

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

1882:    Logically Collective

1884:    Input Parameters:
1885: +  x - the vector
1886: -  a - location of pointer to array obtained from `VecGetArray()`

1888:    Level: beginner

1890:    Fortran Note:
1891:    `VecRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayF90()`

1893: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
1894:           `VecGetArrayPair()`, `VecRestoreArrayPair()`
1895: @*/
1896: PetscErrorCode VecRestoreArray(Vec x, PetscScalar **a)
1897: {
1898:   PetscFunctionBegin;
1901:   if (x->ops->restorearray) {
1902:     PetscUseTypeMethod(x, restorearray, a);
1903:   } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
1904:   if (a) *a = NULL;
1905:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
1906:   PetscFunctionReturn(PETSC_SUCCESS);
1907: }
1908: /*@C
1909:    VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.

1911:    Not Collective

1913:    Input Parameter:
1914: .  x - the vector

1916:    Output Parameter:
1917: .  a - the array

1919:    Level: beginner

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

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

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

1930:    Fortran Note:
1931:    `VecGetArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayReadF90()`

1933: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
1934: @*/
1935: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar **a)
1936: {
1937:   PetscFunctionBegin;
1940:   if (x->ops->getarrayread) {
1941:     PetscUseTypeMethod(x, getarrayread, a);
1942:   } else if (x->ops->getarray) {
1943:     PetscObjectState state;

1945:     /* VECNEST, VECCUDA, VECKOKKOS etc */
1946:     // x->ops->getarray may bump the object state, but since we know this is a read-only get
1947:     // we can just undo that
1948:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
1949:     PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
1950:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
1951:   } else if (x->petscnative) {
1952:     /* VECSTANDARD */
1953:     *a = *((PetscScalar **)x->data);
1954:   } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
1955:   PetscFunctionReturn(PETSC_SUCCESS);
1956: }

1958: /*@C
1959:    VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`

1961:    Not Collective

1963:    Input Parameters:
1964: +  vec - the vector
1965: -  array - the array

1967:    Level: beginner

1969:    Fortran Note:
1970:    `VecRestoreArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayReadF90()`

1972: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
1973: @*/
1974: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar **a)
1975: {
1976:   PetscFunctionBegin;
1979:   if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
1980:     /* nothing */
1981:   } else if (x->ops->restorearrayread) { /* VECNEST */
1982:     PetscUseTypeMethod(x, restorearrayread, a);
1983:   } else { /* No one? */
1984:     PetscObjectState state;

1986:     // x->ops->restorearray may bump the object state, but since we know this is a read-restore
1987:     // we can just undo that
1988:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
1989:     PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
1990:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
1991:   }
1992:   if (a) *a = NULL;
1993:   PetscFunctionReturn(PETSC_SUCCESS);
1994: }

1996: /*@C
1997:    VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
1998:    processor's portion of the vector data. The values in this array are NOT valid, the caller of this
1999:    routine is responsible for putting values into the array; any values it does not set will be invalid

2001:    Logically Collective

2003:    Input Parameter:
2004: .  x - the vector

2006:    Output Parameter:
2007: .  a - location to put pointer to the array

2009:    Level: intermediate

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

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

2017:    Fortran Note:
2018:    `VecGetArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayWriteF90()`

2020: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteF90()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
2021:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
2022: @*/
2023: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar **a)
2024: {
2025:   PetscFunctionBegin;
2028:   PetscCall(VecSetErrorIfLocked(x, 1));
2029:   if (x->ops->getarraywrite) {
2030:     PetscUseTypeMethod(x, getarraywrite, a);
2031:   } else {
2032:     PetscCall(VecGetArray(x, a));
2033:   }
2034:   PetscFunctionReturn(PETSC_SUCCESS);
2035: }

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

2040:    Logically Collective

2042:    Input Parameters:
2043: +  x - the vector
2044: -  a - location of pointer to array obtained from `VecGetArray()`

2046:    Level: beginner

2048:    Fortran Note:
2049:    `VecRestoreArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayWriteF90()`

2051: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteF90()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2052:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2053: @*/
2054: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar **a)
2055: {
2056:   PetscFunctionBegin;
2059:   if (x->ops->restorearraywrite) {
2060:     PetscUseTypeMethod(x, restorearraywrite, a);
2061:   } else if (x->ops->restorearray) {
2062:     PetscUseTypeMethod(x, restorearray, a);
2063:   }
2064:   if (a) *a = NULL;
2065:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
2066:   PetscFunctionReturn(PETSC_SUCCESS);
2067: }

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

2073:    Logically Collective; No Fortran Support

2075:    Input Parameters:
2076: +  x - the vectors
2077: -  n - the number of vectors

2079:    Output Parameter:
2080: .  a - location to put pointer to the array

2082:    Level: intermediate

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

2087: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2088: @*/
2089: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2090: {
2091:   PetscInt      i;
2092:   PetscScalar **q;

2094:   PetscFunctionBegin;
2098:   PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2099:   PetscCall(PetscMalloc1(n, &q));
2100:   for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2101:   *a = q;
2102:   PetscFunctionReturn(PETSC_SUCCESS);
2103: }

2105: /*@C
2106:    VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2107:    has been called.

2109:    Logically Collective; No Fortran Support

2111:    Input Parameters:
2112: +  x - the vector
2113: .  n - the number of vectors
2114: -  a - location of pointer to arrays obtained from `VecGetArrays()`

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

2122:    Level: intermediate

2124: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2125: @*/
2126: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2127: {
2128:   PetscInt      i;
2129:   PetscScalar **q = *a;

2131:   PetscFunctionBegin;

2136:   for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2137:   PetscCall(PetscFree(q));
2138:   PetscFunctionReturn(PETSC_SUCCESS);
2139: }

2141: /*@C
2142:    VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g., `VECCUDA`), the returned pointer will be a device
2143:    pointer to the device memory that contains this processor's portion of the vector data. Device data is guaranteed to have the latest value.
2144:    Otherwise, when this is a host vector (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host pointer.

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

2149:    Logically Collective; No Fortran Support

2151:    Input Parameter:
2152: .  x - the vector

2154:    Output Parameters:
2155: +  a - location to put pointer to the array
2156: -  mtype - memory type of the array

2158:    Level: beginner

2160:    Note:
2161:    Use `VecRestoreArrayAndMemType()` when the array access is no longer needed

2163: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`,
2164:           `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2165: @*/
2166: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2167: {
2168:   PetscFunctionBegin;
2173:   PetscCall(VecSetErrorIfLocked(x, 1));
2174:   if (x->ops->getarrayandmemtype) {
2175:     /* VECCUDA, VECKOKKOS etc */
2176:     PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2177:   } else {
2178:     /* VECSTANDARD, VECNEST, VECVIENNACL */
2179:     PetscCall(VecGetArray(x, a));
2180:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2181:   }
2182:   PetscFunctionReturn(PETSC_SUCCESS);
2183: }

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

2188:    Logically Collective; No Fortran Support

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

2194:    Level: beginner

2196: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`,
2197:           `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2198: @*/
2199: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar **a)
2200: {
2201:   PetscFunctionBegin;
2205:   if (x->ops->restorearrayandmemtype) {
2206:     /* VECCUDA, VECKOKKOS etc */
2207:     PetscUseTypeMethod(x, restorearrayandmemtype, a);
2208:   } else {
2209:     /* VECNEST, VECVIENNACL */
2210:     PetscCall(VecRestoreArray(x, a));
2211:   } /* VECSTANDARD does nothing */
2212:   if (a) *a = NULL;
2213:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
2214:   PetscFunctionReturn(PETSC_SUCCESS);
2215: }

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

2221:    Not Collective; No Fortran Support

2223:    Input Parameter:
2224: .  x - the vector

2226:    Output Parameters:
2227: +  a - the array
2228: -  mtype - memory type of the array

2230:    Level: beginner

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

2235: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayAndMemType()`
2236: @*/
2237: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar **a, PetscMemType *mtype)
2238: {
2239:   PetscFunctionBegin;
2244:   if (x->ops->getarrayreadandmemtype) {
2245:     /* VECCUDA/VECHIP though they are also petscnative */
2246:     PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2247:   } else if (x->ops->getarrayandmemtype) {
2248:     /* VECKOKKOS */
2249:     PetscObjectState state;

2251:     // see VecGetArrayRead() for why
2252:     PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2253:     PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2254:     PetscCall(PetscObjectStateSet((PetscObject)x, state));
2255:   } else {
2256:     PetscCall(VecGetArrayRead(x, a));
2257:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2258:   }
2259:   PetscFunctionReturn(PETSC_SUCCESS);
2260: }

2262: /*@C
2263:    VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`

2265:    Not Collective; No Fortran Support

2267:    Input Parameters:
2268: +  vec - the vector
2269: -  array - the array

2271:    Level: beginner

2273: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2274: @*/
2275: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar **a)
2276: {
2277:   PetscFunctionBegin;
2281:   if (x->ops->restorearrayreadandmemtype) {
2282:     /* VECCUDA/VECHIP */
2283:     PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2284:   } else if (!x->petscnative) {
2285:     /* VECNEST */
2286:     PetscCall(VecRestoreArrayRead(x, a));
2287:   }
2288:   if (a) *a = NULL;
2289:   PetscFunctionReturn(PETSC_SUCCESS);
2290: }

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

2296:    Not Collective; No Fortran Support

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

2301:    Output Parameters:
2302: +  a - the array
2303: -  mtype - memory type of the array

2305:    Level: beginner

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

2310: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2311: @*/
2312: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2313: {
2314:   PetscFunctionBegin;
2317:   PetscCall(VecSetErrorIfLocked(x, 1));
2320:   if (x->ops->getarraywriteandmemtype) {
2321:     /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2322:     PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2323:   } else if (x->ops->getarrayandmemtype) {
2324:     PetscCall(VecGetArrayAndMemType(x, a, mtype));
2325:   } else {
2326:     /* VECNEST, VECVIENNACL */
2327:     PetscCall(VecGetArrayWrite(x, a));
2328:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2329:   }
2330:   PetscFunctionReturn(PETSC_SUCCESS);
2331: }

2333: /*@C
2334:    VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`

2336:    Not Collective; No Fortran Support

2338:    Input Parameters:
2339: +  vec - the vector
2340: -  array - the array

2342:    Level: beginner

2344: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2345: @*/
2346: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar **a)
2347: {
2348:   PetscFunctionBegin;
2351:   PetscCall(VecSetErrorIfLocked(x, 1));
2353:   if (x->ops->restorearraywriteandmemtype) {
2354:     /* VECCUDA/VECHIP */
2355:     PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2356:     PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2357:   } else if (x->ops->restorearrayandmemtype) {
2358:     PetscCall(VecRestoreArrayAndMemType(x, a));
2359:   } else {
2360:     PetscCall(VecRestoreArray(x, a));
2361:   }
2362:   if (a) *a = NULL;
2363:   PetscFunctionReturn(PETSC_SUCCESS);
2364: }

2366: /*@
2367:    VecPlaceArray - Allows one to replace the array in a vector with an
2368:    array provided by the user. This is useful to avoid copying an array
2369:    into a vector.

2371:    Not Collective; No Fortran Support

2373:    Input Parameters:
2374: +  vec - the vector
2375: -  array - the array

2377:    Level: developer

2379:    Notes:
2380:    Use `VecReplaceArray()` instead to permanently replace the array

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

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

2389: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2390: @*/
2391: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2392: {
2393:   PetscFunctionBegin;
2397:   PetscUseTypeMethod(vec, placearray, array);
2398:   PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2399:   PetscFunctionReturn(PETSC_SUCCESS);
2400: }

2402: /*@C
2403:    VecReplaceArray - Allows one to replace the array in a vector with an
2404:    array provided by the user. This is useful to avoid copying an array
2405:    into a vector.

2407:    Not Collective; No Fortran Support

2409:    Input Parameters:
2410: +  vec - the vector
2411: -  array - the array

2413:    Level: developer

2415:    Notes:
2416:    This permanently replaces the array and frees the memory associated
2417:    with the old array. Use `VecPlaceArray()` to temporarily replace the array.

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

2422: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2423: @*/
2424: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2425: {
2426:   PetscFunctionBegin;
2429:   PetscUseTypeMethod(vec, replacearray, array);
2430:   PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2431:   PetscFunctionReturn(PETSC_SUCCESS);
2432: }

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

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

2441:     Collective

2443:     Input Parameters:
2444: +   x - a vector to mimic
2445: -   n - the number of vectors to obtain

2447:     Output Parameters:
2448: +   y - Fortran pointer to the array of vectors
2449: -   ierr - error code

2451:     Example of Usage:
2452: .vb
2453: #include <petsc/finclude/petscvec.h>
2454:     use petscvec

2456:     Vec x
2457:     Vec, pointer :: y(:)
2458:     ....
2459:     call VecDuplicateVecsF90(x,2,y,ierr)
2460:     call VecSet(y(2),alpha,ierr)
2461:     call VecSet(y(2),alpha,ierr)
2462:     ....
2463:     call VecDestroyVecsF90(2,y,ierr)
2464: .ve

2466:     Level: beginner

2468:     Note:
2469:     Use `VecDestroyVecsF90()` to free the space.

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

2474: /*MC
2475:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2476:     `VecGetArrayF90()`.

2478:     Synopsis:
2479:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2481:     Logically Collective

2483:     Input Parameters:
2484: +   x - vector
2485: -   xx_v - the Fortran pointer to the array

2487:     Output Parameter:
2488: .   ierr - error code

2490:     Example of Usage:
2491: .vb
2492: #include <petsc/finclude/petscvec.h>
2493:     use petscvec

2495:     PetscScalar, pointer :: xx_v(:)
2496:     ....
2497:     call VecGetArrayF90(x,xx_v,ierr)
2498:     xx_v(3) = a
2499:     call VecRestoreArrayF90(x,xx_v,ierr)
2500: .ve

2502:     Level: beginner

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

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

2510:     Synopsis:
2511:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

2513:     Collective

2515:     Input Parameters:
2516: +   n - the number of vectors previously obtained
2517: -   x - pointer to array of vector pointers

2519:     Output Parameter:
2520: .   ierr - error code

2522:     Level: beginner

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

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

2533:     Synopsis:
2534:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2536:     Logically Collective

2538:     Input Parameter:
2539: .   x - vector

2541:     Output Parameters:
2542: +   xx_v - the Fortran pointer to the array
2543: -   ierr - error code

2545:     Example of Usage:
2546: .vb
2547: #include <petsc/finclude/petscvec.h>
2548:     use petscvec

2550:     PetscScalar, pointer :: xx_v(:)
2551:     ....
2552:     call VecGetArrayF90(x,xx_v,ierr)
2553:     xx_v(3) = a
2554:     call VecRestoreArrayF90(x,xx_v,ierr)
2555: .ve

2557:      Level: beginner

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

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

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

2571:     Synopsis:
2572:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2574:     Logically Collective

2576:     Input Parameter:
2577: .   x - vector

2579:     Output Parameters:
2580: +   xx_v - the Fortran pointer to the array
2581: -   ierr - error code

2583:     Example of Usage:
2584: .vb
2585: #include <petsc/finclude/petscvec.h>
2586:     use petscvec

2588:     PetscScalar, pointer :: xx_v(:)
2589:     ....
2590:     call VecGetArrayReadF90(x,xx_v,ierr)
2591:     a = xx_v(3)
2592:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2593: .ve

2595:     Level: beginner

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

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

2603: /*MC
2604:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2605:     `VecGetArrayReadF90()`.

2607:     Synopsis:
2608:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2610:     Logically Collective

2612:     Input Parameters:
2613: +   x - vector
2614: -   xx_v - the Fortran pointer to the array

2616:     Output Parameter:
2617: .   ierr - error code

2619:     Example of Usage:
2620: .vb
2621: #include <petsc/finclude/petscvec.h>
2622:     use petscvec

2624:     PetscScalar, pointer :: xx_v(:)
2625:     ....
2626:     call VecGetArrayReadF90(x,xx_v,ierr)
2627:     a = xx_v(3)
2628:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2629: .ve

2631:     Level: beginner

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

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

2641:    Logically Collective

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

2650:    Output Parameter:
2651: .  a - location to put pointer to the array

2653:    Level: developer

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

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

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

2672:   PetscFunctionBegin;
2676:   PetscCall(VecGetLocalSize(x, &N));
2677:   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);
2678:   PetscCall(VecGetArray(x, &aa));

2680:   PetscCall(PetscMalloc1(m, a));
2681:   for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2682:   *a -= mstart;
2683:   PetscFunctionReturn(PETSC_SUCCESS);
2684: }

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

2691:    Logically Collective

2693:    Input Parameters:
2694: +  x - the vector
2695: .  m - first dimension of two dimensional array
2696: .  n - second dimension of two dimensional array
2697: .  mstart - first index you will use in first coordinate direction (often 0)
2698: -  nstart - first index in the second coordinate direction (often 0)

2700:    Output Parameter:
2701: .  a - location to put pointer to the array

2703:    Level: developer

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

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

2713: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2714:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2715:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2716: @*/
2717: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2718: {
2719:   PetscInt     i, N;
2720:   PetscScalar *aa;

2722:   PetscFunctionBegin;
2726:   PetscCall(VecGetLocalSize(x, &N));
2727:   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);
2728:   PetscCall(VecGetArrayWrite(x, &aa));

2730:   PetscCall(PetscMalloc1(m, a));
2731:   for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2732:   *a -= mstart;
2733:   PetscFunctionReturn(PETSC_SUCCESS);
2734: }

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

2739:    Logically Collective

2741:    Input Parameters:
2742: +  x - the vector
2743: .  m - first dimension of two dimensional array
2744: .  n - second dimension of the two dimensional array
2745: .  mstart - first index you will use in first coordinate direction (often 0)
2746: .  nstart - first index in the second coordinate direction (often 0)
2747: -  a - location of pointer to array obtained from `VecGetArray2d()`

2749:    Level: developer

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

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

2759: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2760:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2761:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2762: @*/
2763: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2764: {
2765:   void *dummy;

2767:   PetscFunctionBegin;
2771:   dummy = (void *)(*a + mstart);
2772:   PetscCall(PetscFree(dummy));
2773:   PetscCall(VecRestoreArray(x, NULL));
2774:   *a = NULL;
2775:   PetscFunctionReturn(PETSC_SUCCESS);
2776: }

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

2781:    Logically Collective

2783:    Input Parameters:
2784: +  x - the vector
2785: .  m - first dimension of two dimensional array
2786: .  n - second dimension of the two dimensional array
2787: .  mstart - first index you will use in first coordinate direction (often 0)
2788: .  nstart - first index in the second coordinate direction (often 0)
2789: -  a - location of pointer to array obtained from `VecGetArray2d()`

2791:    Level: developer

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

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

2801: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2802:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2803:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2804: @*/
2805: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2806: {
2807:   void *dummy;

2809:   PetscFunctionBegin;
2813:   dummy = (void *)(*a + mstart);
2814:   PetscCall(PetscFree(dummy));
2815:   PetscCall(VecRestoreArrayWrite(x, NULL));
2816:   PetscFunctionReturn(PETSC_SUCCESS);
2817: }

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

2824:    Logically Collective

2826:    Input Parameters:
2827: +  x - the vector
2828: .  m - first dimension of two dimensional array
2829: -  mstart - first index you will use in first coordinate direction (often 0)

2831:    Output Parameter:
2832: .  a - location to put pointer to the array

2834:    Level: developer

2836:   Notes:
2837:    For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2838:    obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2839:    `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

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

2843: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2844:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2845:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2846: @*/
2847: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2848: {
2849:   PetscInt N;

2851:   PetscFunctionBegin;
2855:   PetscCall(VecGetLocalSize(x, &N));
2856:   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);
2857:   PetscCall(VecGetArray(x, a));
2858:   *a -= mstart;
2859:   PetscFunctionReturn(PETSC_SUCCESS);
2860: }

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

2867:    Logically Collective

2869:    Input Parameters:
2870: +  x - the vector
2871: .  m - first dimension of two dimensional array
2872: -  mstart - first index you will use in first coordinate direction (often 0)

2874:    Output Parameter:
2875: .  a - location to put pointer to the array

2877:    Level: developer

2879:   Notes:
2880:    For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2881:    obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2882:    `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

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

2886: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2887:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2888:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2889: @*/
2890: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2891: {
2892:   PetscInt N;

2894:   PetscFunctionBegin;
2898:   PetscCall(VecGetLocalSize(x, &N));
2899:   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);
2900:   PetscCall(VecGetArrayWrite(x, a));
2901:   *a -= mstart;
2902:   PetscFunctionReturn(PETSC_SUCCESS);
2903: }

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

2908:    Logically Collective

2910:    Input Parameters:
2911: +  x - the vector
2912: .  m - first dimension of two dimensional array
2913: .  mstart - first index you will use in first coordinate direction (often 0)
2914: -  a - location of pointer to array obtained from `VecGetArray1d()`

2916:    Level: developer

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

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

2926: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2927:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2928:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2929: @*/
2930: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2931: {
2932:   PetscFunctionBegin;
2935:   PetscCall(VecRestoreArray(x, NULL));
2936:   *a = NULL;
2937:   PetscFunctionReturn(PETSC_SUCCESS);
2938: }

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

2943:    Logically Collective

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

2951:    Level: developer

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

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

2961: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2962:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2963:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2964: @*/
2965: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2966: {
2967:   PetscFunctionBegin;
2970:   PetscCall(VecRestoreArrayWrite(x, NULL));
2971:   *a = NULL;
2972:   PetscFunctionReturn(PETSC_SUCCESS);
2973: }

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

2980:    Logically Collective

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

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

2994:    Level: developer

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

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

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

3013:   PetscFunctionBegin;
3017:   PetscCall(VecGetLocalSize(x, &N));
3018:   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);
3019:   PetscCall(VecGetArray(x, &aa));

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

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

3035:    Logically Collective

3037:    Input Parameters:
3038: +  x - the vector
3039: .  m - first dimension of three dimensional array
3040: .  n - second dimension of three dimensional array
3041: .  p - third dimension of 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)

3046:    Output Parameter:
3047: .  a - location to put pointer to the array

3049:    Level: developer

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

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

3059: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3060:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3061:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3062: @*/
3063: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3064: {
3065:   PetscInt     i, N, j;
3066:   PetscScalar *aa, **b;

3068:   PetscFunctionBegin;
3072:   PetscCall(VecGetLocalSize(x, &N));
3073:   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);
3074:   PetscCall(VecGetArrayWrite(x, &aa));

3076:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3077:   b = (PetscScalar **)((*a) + m);
3078:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3079:   for (i = 0; i < m; i++)
3080:     for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;

3082:   *a -= mstart;
3083:   PetscFunctionReturn(PETSC_SUCCESS);
3084: }

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

3089:    Logically Collective

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

3101:    Level: developer

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

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

3111: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3112:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3113:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3114: @*/
3115: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3116: {
3117:   void *dummy;

3119:   PetscFunctionBegin;
3123:   dummy = (void *)(*a + mstart);
3124:   PetscCall(PetscFree(dummy));
3125:   PetscCall(VecRestoreArray(x, NULL));
3126:   *a = NULL;
3127:   PetscFunctionReturn(PETSC_SUCCESS);
3128: }

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

3133:    Logically Collective

3135:    Input Parameters:
3136: +  x - the vector
3137: .  m - first dimension of three dimensional array
3138: .  n - second dimension of the three dimensional array
3139: .  p - third dimension of the three dimensional array
3140: .  mstart - first index you will use in first coordinate direction (often 0)
3141: .  nstart - first index in the second coordinate direction (often 0)
3142: .  pstart - first index in the third coordinate direction (often 0)
3143: -  a - location of pointer to array obtained from VecGetArray3d()

3145:    Level: developer

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

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

3155: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3156:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3157:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3158: @*/
3159: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3160: {
3161:   void *dummy;

3163:   PetscFunctionBegin;
3167:   dummy = (void *)(*a + mstart);
3168:   PetscCall(PetscFree(dummy));
3169:   PetscCall(VecRestoreArrayWrite(x, NULL));
3170:   *a = NULL;
3171:   PetscFunctionReturn(PETSC_SUCCESS);
3172: }

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

3179:    Logically Collective

3181:    Input Parameters:
3182: +  x - the vector
3183: .  m - first dimension of four dimensional array
3184: .  n - second dimension of four dimensional array
3185: .  p - third dimension of four dimensional array
3186: .  q - fourth dimension of four dimensional array
3187: .  mstart - first index you will use in first coordinate direction (often 0)
3188: .  nstart - first index in the second coordinate direction (often 0)
3189: .  pstart - first index in the third coordinate direction (often 0)
3190: -  qstart - first index in the fourth coordinate direction (often 0)

3192:    Output Parameter:
3193: .  a - location to put pointer to the array

3195:    Level: developer

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

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

3205: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3206:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3207:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3208: @*/
3209: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3210: {
3211:   PetscInt     i, N, j, k;
3212:   PetscScalar *aa, ***b, **c;

3214:   PetscFunctionBegin;
3218:   PetscCall(VecGetLocalSize(x, &N));
3219:   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);
3220:   PetscCall(VecGetArray(x, &aa));

3222:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3223:   b = (PetscScalar ***)((*a) + m);
3224:   c = (PetscScalar **)(b + m * n);
3225:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3226:   for (i = 0; i < m; i++)
3227:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3228:   for (i = 0; i < m; i++)
3229:     for (j = 0; j < n; j++)
3230:       for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3231:   *a -= mstart;
3232:   PetscFunctionReturn(PETSC_SUCCESS);
3233: }

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

3240:    Logically Collective

3242:    Input Parameters:
3243: +  x - the vector
3244: .  m - first dimension of four dimensional array
3245: .  n - second dimension of four dimensional array
3246: .  p - third dimension of four dimensional array
3247: .  q - fourth dimension of four dimensional array
3248: .  mstart - first index you will use in first coordinate direction (often 0)
3249: .  nstart - first index in the second coordinate direction (often 0)
3250: .  pstart - first index in the third coordinate direction (often 0)
3251: -  qstart - first index in the fourth coordinate direction (often 0)

3253:    Output Parameter:
3254: .  a - location to put pointer to the array

3256:    Level: developer

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

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

3266: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3267:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3268:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3269: @*/
3270: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3271: {
3272:   PetscInt     i, N, j, k;
3273:   PetscScalar *aa, ***b, **c;

3275:   PetscFunctionBegin;
3279:   PetscCall(VecGetLocalSize(x, &N));
3280:   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);
3281:   PetscCall(VecGetArrayWrite(x, &aa));

3283:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3284:   b = (PetscScalar ***)((*a) + m);
3285:   c = (PetscScalar **)(b + m * n);
3286:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3287:   for (i = 0; i < m; i++)
3288:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3289:   for (i = 0; i < m; i++)
3290:     for (j = 0; j < n; j++)
3291:       for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3292:   *a -= mstart;
3293:   PetscFunctionReturn(PETSC_SUCCESS);
3294: }

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

3299:    Logically Collective

3301:    Input Parameters:
3302: +  x - the vector
3303: .  m - first dimension of four dimensional array
3304: .  n - second dimension of the four dimensional array
3305: .  p - third dimension of the four dimensional array
3306: .  q - fourth dimension of the four dimensional array
3307: .  mstart - first index you will use in first coordinate direction (often 0)
3308: .  nstart - first index in the second coordinate direction (often 0)
3309: .  pstart - first index in the third coordinate direction (often 0)
3310: .  qstart - first index in the fourth coordinate direction (often 0)
3311: -  a - location of pointer to array obtained from VecGetArray4d()

3313:    Level: developer

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

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

3323: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3324:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3325:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3326: @*/
3327: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3328: {
3329:   void *dummy;

3331:   PetscFunctionBegin;
3335:   dummy = (void *)(*a + mstart);
3336:   PetscCall(PetscFree(dummy));
3337:   PetscCall(VecRestoreArray(x, NULL));
3338:   *a = NULL;
3339:   PetscFunctionReturn(PETSC_SUCCESS);
3340: }

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

3345:    Logically Collective

3347:    Input Parameters:
3348: +  x - the vector
3349: .  m - first dimension of four dimensional array
3350: .  n - second dimension of the four dimensional array
3351: .  p - third dimension of the four dimensional array
3352: .  q - fourth dimension of the four dimensional array
3353: .  mstart - first index you will use in first coordinate direction (often 0)
3354: .  nstart - first index in the second coordinate direction (often 0)
3355: .  pstart - first index in the third coordinate direction (often 0)
3356: .  qstart - first index in the fourth coordinate direction (often 0)
3357: -  a - location of pointer to array obtained from `VecGetArray4d()`

3359:    Level: developer

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

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

3369: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3370:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3371:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3372: @*/
3373: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3374: {
3375:   void *dummy;

3377:   PetscFunctionBegin;
3381:   dummy = (void *)(*a + mstart);
3382:   PetscCall(PetscFree(dummy));
3383:   PetscCall(VecRestoreArrayWrite(x, NULL));
3384:   *a = NULL;
3385:   PetscFunctionReturn(PETSC_SUCCESS);
3386: }

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

3393:    Logically Collective

3395:    Input Parameters:
3396: +  x - the vector
3397: .  m - first dimension of two dimensional array
3398: .  n - second dimension of two dimensional array
3399: .  mstart - first index you will use in first coordinate direction (often 0)
3400: -  nstart - first index in the second coordinate direction (often 0)

3402:    Output Parameter:
3403: .  a - location to put pointer to the array

3405:    Level: developer

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

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

3415: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3416:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3417:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3418: @*/
3419: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3420: {
3421:   PetscInt           i, N;
3422:   const PetscScalar *aa;

3424:   PetscFunctionBegin;
3428:   PetscCall(VecGetLocalSize(x, &N));
3429:   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);
3430:   PetscCall(VecGetArrayRead(x, &aa));

3432:   PetscCall(PetscMalloc1(m, a));
3433:   for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3434:   *a -= mstart;
3435:   PetscFunctionReturn(PETSC_SUCCESS);
3436: }

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

3441:    Logically Collective

3443:    Input Parameters:
3444: +  x - the vector
3445: .  m - first dimension of two dimensional array
3446: .  n - second dimension of the two dimensional array
3447: .  mstart - first index you will use in first coordinate direction (often 0)
3448: .  nstart - first index in the second coordinate direction (often 0)
3449: -  a - location of pointer to array obtained from VecGetArray2d()

3451:    Level: developer

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

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

3461: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3462:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3463:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3464: @*/
3465: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3466: {
3467:   void *dummy;

3469:   PetscFunctionBegin;
3473:   dummy = (void *)(*a + mstart);
3474:   PetscCall(PetscFree(dummy));
3475:   PetscCall(VecRestoreArrayRead(x, NULL));
3476:   *a = NULL;
3477:   PetscFunctionReturn(PETSC_SUCCESS);
3478: }

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

3485:    Logically Collective

3487:    Input Parameters:
3488: +  x - the vector
3489: .  m - first dimension of two dimensional array
3490: -  mstart - first index you will use in first coordinate direction (often 0)

3492:    Output Parameter:
3493: .  a - location to put pointer to the array

3495:    Level: developer

3497:   Notes:
3498:    For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3499:    obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3500:    `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.

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

3504: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3505:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3506:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3507: @*/
3508: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3509: {
3510:   PetscInt N;

3512:   PetscFunctionBegin;
3516:   PetscCall(VecGetLocalSize(x, &N));
3517:   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);
3518:   PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3519:   *a -= mstart;
3520:   PetscFunctionReturn(PETSC_SUCCESS);
3521: }

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

3526:    Logically Collective

3528:    Input Parameters:
3529: +  x - the vector
3530: .  m - first dimension of two dimensional array
3531: .  mstart - first index you will use in first coordinate direction (often 0)
3532: -  a - location of pointer to array obtained from `VecGetArray1dRead()`

3534:    Level: developer

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

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

3544: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3545:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3546:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3547: @*/
3548: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3549: {
3550:   PetscFunctionBegin;
3553:   PetscCall(VecRestoreArrayRead(x, NULL));
3554:   *a = NULL;
3555:   PetscFunctionReturn(PETSC_SUCCESS);
3556: }

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

3563:    Logically Collective

3565:    Input Parameters:
3566: +  x - the vector
3567: .  m - first dimension of three dimensional array
3568: .  n - second dimension of three dimensional array
3569: .  p - third dimension of 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)

3574:    Output Parameter:
3575: .  a - location to put pointer to the array

3577:    Level: developer

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

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

3587: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3588:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3589:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3590: @*/
3591: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3592: {
3593:   PetscInt           i, N, j;
3594:   const PetscScalar *aa;
3595:   PetscScalar      **b;

3597:   PetscFunctionBegin;
3601:   PetscCall(VecGetLocalSize(x, &N));
3602:   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);
3603:   PetscCall(VecGetArrayRead(x, &aa));

3605:   PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3606:   b = (PetscScalar **)((*a) + m);
3607:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3608:   for (i = 0; i < m; i++)
3609:     for (j = 0; j < n; j++) b[i * n + j] = (PetscScalar *)aa + i * n * p + j * p - pstart;
3610:   *a -= mstart;
3611:   PetscFunctionReturn(PETSC_SUCCESS);
3612: }

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

3617:    Logically Collective

3619:    Input Parameters:
3620: +  x - the vector
3621: .  m - first dimension of three dimensional array
3622: .  n - second dimension of the three dimensional array
3623: .  p - third dimension of the three dimensional array
3624: .  mstart - first index you will use in first coordinate direction (often 0)
3625: .  nstart - first index in the second coordinate direction (often 0)
3626: .  pstart - first index in the third coordinate direction (often 0)
3627: -  a - location of pointer to array obtained from `VecGetArray3dRead()`

3629:    Level: developer

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

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

3639: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3640:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3641:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3642: @*/
3643: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3644: {
3645:   void *dummy;

3647:   PetscFunctionBegin;
3651:   dummy = (void *)(*a + mstart);
3652:   PetscCall(PetscFree(dummy));
3653:   PetscCall(VecRestoreArrayRead(x, NULL));
3654:   *a = NULL;
3655:   PetscFunctionReturn(PETSC_SUCCESS);
3656: }

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

3663:    Logically Collective

3665:    Input Parameters:
3666: +  x - the vector
3667: .  m - first dimension of four dimensional array
3668: .  n - second dimension of four dimensional array
3669: .  p - third dimension of four dimensional array
3670: .  q - fourth dimension of four dimensional array
3671: .  mstart - first index you will use in first coordinate direction (often 0)
3672: .  nstart - first index in the second coordinate direction (often 0)
3673: .  pstart - first index in the third coordinate direction (often 0)
3674: -  qstart - first index in the fourth coordinate direction (often 0)

3676:    Output Parameter:
3677: .  a - location to put pointer to the array

3679:    Level: beginner

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

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

3689: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3690:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3691:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3692: @*/
3693: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3694: {
3695:   PetscInt           i, N, j, k;
3696:   const PetscScalar *aa;
3697:   PetscScalar     ***b, **c;

3699:   PetscFunctionBegin;
3703:   PetscCall(VecGetLocalSize(x, &N));
3704:   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);
3705:   PetscCall(VecGetArrayRead(x, &aa));

3707:   PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3708:   b = (PetscScalar ***)((*a) + m);
3709:   c = (PetscScalar **)(b + m * n);
3710:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3711:   for (i = 0; i < m; i++)
3712:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3713:   for (i = 0; i < m; i++)
3714:     for (j = 0; j < n; j++)
3715:       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;
3716:   *a -= mstart;
3717:   PetscFunctionReturn(PETSC_SUCCESS);
3718: }

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

3723:    Logically Collective

3725:    Input Parameters:
3726: +  x - the vector
3727: .  m - first dimension of four dimensional array
3728: .  n - second dimension of the four dimensional array
3729: .  p - third dimension of the four dimensional array
3730: .  q - fourth dimension of the four dimensional array
3731: .  mstart - first index you will use in first coordinate direction (often 0)
3732: .  nstart - first index in the second coordinate direction (often 0)
3733: .  pstart - first index in the third coordinate direction (often 0)
3734: .  qstart - first index in the fourth coordinate direction (often 0)
3735: -  a - location of pointer to array obtained from `VecGetArray4dRead()`

3737:    Level: beginner

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

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

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

3755:   PetscFunctionBegin;
3759:   dummy = (void *)(*a + mstart);
3760:   PetscCall(PetscFree(dummy));
3761:   PetscCall(VecRestoreArrayRead(x, NULL));
3762:   *a = NULL;
3763:   PetscFunctionReturn(PETSC_SUCCESS);
3764: }

3766: #if defined(PETSC_USE_DEBUG)

3768: /*@
3769:    VecLockGet  - Gets the current lock status of a vector

3771:    Logically Collective

3773:    Input Parameter:
3774: .  x - the vector

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

3780:    Level: advanced

3782: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3783: @*/
3784: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3785: {
3786:   PetscFunctionBegin;
3788:   *state = x->lock;
3789:   PetscFunctionReturn(PETSC_SUCCESS);
3790: }

3792: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3793: {
3794:   PetscFunctionBegin;
3799:   #if !PetscDefined(HAVE_THREADSAFETY)
3800:   {
3801:     const int index = x->lockstack.currentsize - 1;

3803:     PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupted vec lock stack, have negative index %d", index);
3804:     *file = x->lockstack.file[index];
3805:     *func = x->lockstack.function[index];
3806:     *line = x->lockstack.line[index];
3807:   }
3808:   #else
3809:   *file = NULL;
3810:   *func = NULL;
3811:   *line = 0;
3812:   #endif
3813:   PetscFunctionReturn(PETSC_SUCCESS);
3814: }

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

3819:    Logically Collective

3821:    Input Parameter:
3822: .  x - the vector

3824:    Level: intermediate

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

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

3832: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3833: @*/
3834: PetscErrorCode VecLockReadPush(Vec x)
3835: {
3836:   PetscFunctionBegin;
3838:   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");
3839:   #if !PetscDefined(HAVE_THREADSAFETY)
3840:   {
3841:     const char *file, *func;
3842:     int         index, line;

3844:     if ((index = petscstack.currentsize - 2) == -1) {
3845:       // vec was locked "outside" of petsc, either in user-land or main. the error message will
3846:       // now show this function as the culprit, but it will include the stacktrace
3847:       file = "unknown user-file";
3848:       func = "unknown_user_function";
3849:       line = 0;
3850:     } else {
3851:       PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected petscstack, have negative index %d", index);
3852:       file = petscstack.file[index];
3853:       func = petscstack.function[index];
3854:       line = petscstack.line[index];
3855:     }
3856:     PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
3857:   }
3858:   #endif
3859:   PetscFunctionReturn(PETSC_SUCCESS);
3860: }

3862: /*@
3863:    VecLockReadPop  - Pops a read-only lock from a vector

3865:    Logically Collective

3867:    Input Parameter:
3868: .  x - the vector

3870:    Level: intermediate

3872: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
3873: @*/
3874: PetscErrorCode VecLockReadPop(Vec x)
3875: {
3876:   PetscFunctionBegin;
3878:   PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
3879:   #if !PetscDefined(HAVE_THREADSAFETY)
3880:   {
3881:     const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];

3883:     PetscStackPop_Private(x->lockstack, previous);
3884:   }
3885:   #endif
3886:   PetscFunctionReturn(PETSC_SUCCESS);
3887: }

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

3892:    Logically Collective

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

3898:    Level: intermediate

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

3906: .vb
3907:        VecGetArray(x,&xdata); // begin phase
3908:        VecLockWriteSet(v,PETSC_TRUE);

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

3912:        VecRestoreArray(x,&vdata); // end phase
3913:        VecLockWriteSet(v,PETSC_FALSE);
3914: .ve

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

3919: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
3920: @*/
3921: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
3922: {
3923:   PetscFunctionBegin;
3925:   if (flg) {
3926:     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");
3927:     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");
3928:     x->lock = -1;
3929:   } else {
3930:     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");
3931:     x->lock = 0;
3932:   }
3933:   PetscFunctionReturn(PETSC_SUCCESS);
3934: }

3936: /*@
3937:    VecLockPush  - Pushes a read-only lock on a vector to prevent it from being written to

3939:    Level: deprecated

3941: .seealso: [](ch_vectors), `Vec`, `VecLockReadPush()`
3942: @*/
3943: PetscErrorCode VecLockPush(Vec x)
3944: {
3945:   PetscFunctionBegin;
3946:   PetscCall(VecLockReadPush(x));
3947:   PetscFunctionReturn(PETSC_SUCCESS);
3948: }

3950: /*@
3951:    VecLockPop  - Pops a read-only lock from a vector

3953:    Level: deprecated

3955: .seealso: [](ch_vectors), `Vec`, `VecLockReadPop()`
3956: @*/
3957: PetscErrorCode VecLockPop(Vec x)
3958: {
3959:   PetscFunctionBegin;
3960:   PetscCall(VecLockReadPop(x));
3961:   PetscFunctionReturn(PETSC_SUCCESS);
3962: }

3964: #endif