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>
  8: #if defined(PETSC_HAVE_CUDA)
  9: #include <../src/vec/vec/impls/dvecimpl.h>
 10: #include <petsc/private/cudavecimpl.h>
 11: #endif
 12: #if defined(PETSC_HAVE_HIP)
 13: #include <../src/vec/vec/impls/dvecimpl.h>
 14: #include <petsc/private/hipvecimpl.h>
 15: #endif
 16: PetscInt VecGetSubVectorSavedStateId = -1;

 18: #if PetscDefined(USE_DEBUG)
 19: // this is a no-op '0' macro in optimized builds
 20: PetscErrorCode VecValidValues_Internal(Vec vec, PetscInt argnum, PetscBool begin)
 21: {
 22:   if (vec->petscnative || vec->ops->getarray) {
 23:     PetscInt           n;
 24:     const PetscScalar *x;
 25:     PetscOffloadMask   mask;

 27:     VecGetOffloadMask(vec, &mask);
 28:     if (!PetscOffloadHost(mask)) return 0;
 29:     VecGetLocalSize(vec, &n);
 30:     VecGetArrayRead(vec, &x);
 31:     for (PetscInt i = 0; i < n; i++) {
 32:       if (begin) {
 34:       } else {
 36:       }
 37:     }
 38:     VecRestoreArrayRead(vec, &x);
 39:   }
 40:   return 0;
 41: }
 42: #endif

 44: /*@
 45:    VecMaxPointwiseDivide - Computes the maximum of the componentwise division max = max_i abs(x_i/y_i).

 47:    Logically Collective on x

 49:    Input Parameters:
 50: .  x, y  - the vectors

 52:    Output Parameter:
 53: .  max - the result

 55:    Level: advanced

 57:    Notes:
 58:    x and y may be the same vector

 60:   if a particular y_i is zero, it is treated as 1 in the above formula

 62: .seealso: `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`
 63: @*/
 64: PetscErrorCode VecMaxPointwiseDivide(Vec x, Vec y, PetscReal *max)
 65: {
 72:   VecCheckSameSize(x, 1, y, 2);
 73:   VecLockReadPush(x);
 74:   VecLockReadPush(y);
 75:   PetscUseTypeMethod(x, maxpointwisedivide, y, max);
 76:   VecLockReadPop(x);
 77:   VecLockReadPop(y);
 78:   return 0;
 79: }

 81: /*@
 82:    VecDot - Computes the vector dot product.

 84:    Collective on x

 86:    Input Parameters:
 87: .  x, y - the vectors

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

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

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

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

110:    Level: intermediate

112: .seealso: `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDotRealPart()`
113: @*/
114: PetscErrorCode VecDot(Vec x, Vec y, PetscScalar *val)
115: {
122:   VecCheckSameSize(x, 1, y, 2);

124:   VecLockReadPush(x);
125:   VecLockReadPush(y);
126:   PetscLogEventBegin(VEC_Dot, x, y, 0, 0);
127:   PetscUseTypeMethod(x, dot, y, val);
128:   PetscLogEventEnd(VEC_Dot, x, y, 0, 0);
129:   VecLockReadPop(x);
130:   VecLockReadPop(y);
131:   return 0;
132: }

134: /*@
135:    VecDotRealPart - Computes the real part of the vector dot product.

137:    Collective on x

139:    Input Parameters:
140: .  x, y - the vectors

142:    Output Parameter:
143: .  val - the real part of the dot product;

145:    Level: intermediate

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

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

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

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

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

165: .seealso: `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDot()`, `VecDotNorm2()`
166: @*/
167: PetscErrorCode VecDotRealPart(Vec x, Vec y, PetscReal *val)
168: {
169:   PetscScalar fdot;

171:   VecDot(x, y, &fdot);
172:   *val = PetscRealPart(fdot);
173:   return 0;
174: }

176: /*@
177:    VecNorm  - Computes the vector norm.

179:    Collective on Vec

181:    Input Parameters:
182: +  x - the vector
183: -  type - the type of the norm requested

185:    Output Parameter:
186: .  val - the norm

188:    Values of NormType:
189: +     NORM_1 - sum_i |x_i|
190: .     NORM_2 - sqrt(sum_i |x_i|^2)
191: .     NORM_INFINITY - max_i |x_i|
192: -     NORM_1_AND_2 - computes efficiently both  NORM_1 and NORM_2 and stores them each in an output array

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

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

205:    Level: intermediate

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

213: .seealso: `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecNormAvailable()`,
214:           `VecNormBegin()`, `VecNormEnd()`, `NormType()`

216: @*/
217: PetscErrorCode VecNorm(Vec x, NormType type, PetscReal *val)
218: {

223:   /* Cached data? */
224:   if (type != NORM_1_AND_2) {
225:     PetscBool flg;

227:     PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, flg);
228:     if (flg) return 0;
229:   }

231:   VecLockReadPush(x);
232:   PetscLogEventBegin(VEC_Norm, x, 0, 0, 0);
233:   PetscUseTypeMethod(x, norm, type, val);
234:   PetscLogEventEnd(VEC_Norm, x, 0, 0, 0);
235:   VecLockReadPop(x);

237:   if (type != NORM_1_AND_2) PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val);
238:   return 0;
239: }

241: /*@
242:    VecNormAvailable  - Returns the vector norm if it is already known.

244:    Not Collective

246:    Input Parameters:
247: +  x - the vector
248: -  type - one of NORM_1, NORM_2, NORM_INFINITY.  Also available
249:           NORM_1_AND_2, which computes both norms and stores them
250:           in a two element array.

252:    Output Parameters:
253: +  available - PETSC_TRUE if the val returned is valid
254: -  val - the norm

256:    Notes:
257: $     NORM_1 denotes sum_i |x_i|
258: $     NORM_2 denotes sqrt(sum_i (x_i)^2)
259: $     NORM_INFINITY denotes max_i |x_i|

261:    Level: intermediate

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

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

273: .seealso: `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecNorm()`
274:           `VecNormBegin()`, `VecNormEnd()`

276: @*/
277: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
278: {

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

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

295:    Collective on Vec

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

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

303:    Level: intermediate

305: @*/
306: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
307: {
308:   PetscReal norm;

312:   VecSetErrorIfLocked(x, 1);
314:   PetscLogEventBegin(VEC_Normalize, x, 0, 0, 0);
315:   VecNorm(x, NORM_2, &norm);
316:   if (norm == 0.0) {
317:     PetscInfo(x, "Vector of zero norm can not be normalized; Returning only the zero norm\n");
318:   } else if (norm != 1.0) {
319:     VecScale(x, 1.0 / norm);
320:   }
321:   PetscLogEventEnd(VEC_Normalize, x, 0, 0, 0);
322:   if (val) *val = norm;
323:   return 0;
324: }

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

329:    Collective on Vec

331:    Input Parameter:
332: .  x - the vector

334:    Output Parameters:
335: +  p - the location of val (pass NULL if you don't want this)
336: -  val - the maximum component

338:    Notes:
339:    Returns the value PETSC_MIN_REAL and negative p if the vector is of length 0.

341:    Returns the smallest index with the maximum value
342:    Level: intermediate

344: .seealso: `VecNorm()`, `VecMin()`
345: @*/
346: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
347: {
352:   VecLockReadPush(x);
353:   PetscLogEventBegin(VEC_Max, x, 0, 0, 0);
354:   PetscUseTypeMethod(x, max, p, val);
355:   PetscLogEventEnd(VEC_Max, x, 0, 0, 0);
356:   VecLockReadPop(x);
357:   return 0;
358: }

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

363:    Collective on Vec

365:    Input Parameter:
366: .  x - the vector

368:    Output Parameters:
369: +  p - the location of val (pass NULL if you don't want this location)
370: -  val - the minimum component

372:    Level: intermediate

374:    Notes:
375:    Returns the value PETSC_MAX_REAL and negative p if the vector is of length 0.

377:    This returns the smallest index with the minumum value

379: .seealso: `VecMax()`
380: @*/
381: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
382: {
387:   VecLockReadPush(x);
388:   PetscLogEventBegin(VEC_Min, x, 0, 0, 0);
389:   PetscUseTypeMethod(x, min, p, val);
390:   PetscLogEventEnd(VEC_Min, x, 0, 0, 0);
391:   VecLockReadPop(x);
392:   return 0;
393: }

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

399:    Collective on Vec

401:    Input Parameters:
402: .  x, y - the vectors

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

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

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

416:    Level: intermediate

418: .seealso: `VecDot()`, `VecMTDot()`
419: @*/
420: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
421: {
428:   VecCheckSameSize(x, 1, y, 2);

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

440: /*@
441:    VecScale - Scales a vector.

443:    Not collective on Vec

445:    Input Parameters:
446: +  x - the vector
447: -  alpha - the scalar

449:    Note:
450:    For a vector with n components, VecScale() computes
451: $      x[i] = alpha * x[i], for i=1,...,n.

453:    Level: intermediate

455: @*/
456: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
457: {
458:   PetscReal norms[4];
459:   PetscBool flgs[4];

464:   VecSetErrorIfLocked(x, 1);
465:   if (alpha == (PetscScalar)1.0) return 0;

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

470:   PetscLogEventBegin(VEC_Scale, x, 0, 0, 0);
471:   PetscUseTypeMethod(x, scale, alpha);
472:   PetscLogEventEnd(VEC_Scale, x, 0, 0, 0);

474:   PetscObjectStateIncrease((PetscObject)x);
475:   /* put the scaled stashed norms back into the Vec */
476:   for (PetscInt i = 0; i < 4; i++) {
477:     if (flgs[i]) PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], PetscAbsScalar(alpha) * norms[i]);
478:   }
479:   return 0;
480: }

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

485:    Logically Collective on x

487:    Input Parameters:
488: +  x  - the vector
489: -  alpha - the scalar

491:    Output Parameter:
492: .  x  - the vector

494:    Level: beginner

496:    Note:
497:    For a vector of dimension n, `VecSet()` computes
498: $     x[i] = alpha, for i=1,...,n,
499:    so that all vector entries then equal the identical
500:    scalar value, alpha.  Use the more general routine
501:    `VecSetValues()` to set different vector entries.

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

506: .seealso: `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
507: @*/
508: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
509: {
514:   VecSetErrorIfLocked(x, 1);

516:   PetscLogEventBegin(VEC_Set, x, 0, 0, 0);
517:   PetscUseTypeMethod(x, set, alpha);
518:   PetscLogEventEnd(VEC_Set, x, 0, 0, 0);
519:   PetscObjectStateIncrease((PetscObject)x);

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

523:   {
524:     PetscReal      val = PetscAbsScalar(alpha);
525:     const PetscInt N   = x->map->N;

527:     if (N == 0) {
528:       PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0l);
529:       PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0);
530:       PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0);
531:       PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0);
532:     } else if (val > PETSC_MAX_REAL / N) {
533:       PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val);
534:     } else {
535:       PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val);
536:       PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val);
537:       val = PetscSqrtReal((PetscReal)N) * val;
538:       PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val);
539:       PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val);
540:     }
541:   }
542:   return 0;
543: }

545: /*@
546:    VecAXPY - Computes y = alpha x + y.

548:    Logically Collective on Vec

550:    Input Parameters:
551: +  alpha - the scalar
552: -  x, y  - the vectors

554:    Output Parameter:
555: .  y - output vector

557:    Level: intermediate

559:    Notes:
560:     x and y MUST be different vectors
561:     This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine

563: $    VecAXPY(y,alpha,x)                   y = alpha x           +      y
564: $    VecAYPX(y,beta,x)                    y =       x           + beta y
565: $    VecAXPBY(y,alpha,beta,x)             y = alpha x           + beta y
566: $    VecWAXPY(w,alpha,x,y)                w = alpha x           +      y
567: $    VecAXPBYPCZ(w,alpha,beta,gamma,x,y)  z = alpha x           + beta y + gamma z
568: $    VecMAXPY(y,nv,alpha[],x[])           y = sum alpha[i] x[i] +      y

570: .seealso: `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
571: @*/
572: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
573: {
579:   VecCheckSameSize(x, 3, y, 1);
582:   if (alpha == (PetscScalar)0.0) return 0;

584:   VecSetErrorIfLocked(y, 1);
585:   VecLockReadPush(x);
586:   PetscLogEventBegin(VEC_AXPY, x, y, 0, 0);
587:   PetscUseTypeMethod(y, axpy, alpha, x);
588:   PetscLogEventEnd(VEC_AXPY, x, y, 0, 0);
589:   VecLockReadPop(x);
590:   PetscObjectStateIncrease((PetscObject)y);
591:   return 0;
592: }

594: /*@
595:    VecAYPX - Computes y = x + beta y.

597:    Logically Collective on Vec

599:    Input Parameters:
600: +  beta - the scalar
601: -  x, y  - the vectors

603:    Output Parameter:
604: .  y - output vector

606:    Level: intermediate

608:    Notes:
609:     x and y MUST be different vectors
610:     The implementation is optimized for beta of -1.0, 0.0, and 1.0

612: .seealso: `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
613: @*/
614: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
615: {
621:   VecCheckSameSize(x, 1, y, 3);
624:   VecSetErrorIfLocked(y, 1);
625:   VecLockReadPush(x);
626:   if (beta == (PetscScalar)0.0) {
627:     VecCopy(x, y);
628:   } else {
629:     PetscLogEventBegin(VEC_AYPX, x, y, 0, 0);
630:     PetscUseTypeMethod(y, aypx, beta, x);
631:     PetscLogEventEnd(VEC_AYPX, x, y, 0, 0);
632:     PetscObjectStateIncrease((PetscObject)y);
633:   }
634:   VecLockReadPop(x);
635:   return 0;
636: }

638: /*@
639:    VecAXPBY - Computes y = alpha x + beta y.

641:    Logically Collective on Vec

643:    Input Parameters:
644: +  alpha,beta - the scalars
645: -  x, y  - the vectors

647:    Output Parameter:
648: .  y - output vector

650:    Level: intermediate

652:    Notes:
653:     x and y MUST be different vectors
654:     The implementation is optimized for alpha and/or beta values of 0.0 and 1.0

656: .seealso: `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
657: @*/
658: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
659: {
665:   VecCheckSameSize(y, 1, x, 4);
669:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) return 0;

671:   VecSetErrorIfLocked(y, 1);
672:   VecLockReadPush(x);
673:   PetscLogEventBegin(VEC_AXPY, y, x, 0, 0);
674:   PetscUseTypeMethod(y, axpby, alpha, beta, x);
675:   PetscLogEventEnd(VEC_AXPY, y, x, 0, 0);
676:   PetscObjectStateIncrease((PetscObject)y);
677:   VecLockReadPop(x);
678:   return 0;
679: }

681: /*@
682:    VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z

684:    Logically Collective on Vec

686:    Input Parameters:
687: +  alpha,beta, gamma - the scalars
688: -  x, y, z  - the vectors

690:    Output Parameter:
691: .  z - output vector

693:    Level: intermediate

695:    Notes:
696:     x, y and z must be different vectors
697:     The implementation is optimized for alpha of 1.0 and gamma of 1.0 or 0.0

699: .seealso: `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
700: @*/
701: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
702: {
711:   VecCheckSameSize(x, 1, y, 5);
712:   VecCheckSameSize(x, 1, z, 6);
718:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) return 0;

720:   VecSetErrorIfLocked(z, 1);
721:   VecLockReadPush(x);
722:   VecLockReadPush(y);
723:   PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0);
724:   PetscUseTypeMethod(z, axpbypcz, alpha, beta, gamma, x, y);
725:   PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0);
726:   PetscObjectStateIncrease((PetscObject)z);
727:   VecLockReadPop(x);
728:   VecLockReadPop(y);
729:   return 0;
730: }

732: /*@
733:    VecWAXPY - Computes w = alpha x + y.

735:    Logically Collective on Vec

737:    Input Parameters:
738: +  alpha - the scalar
739: -  x, y  - the vectors

741:    Output Parameter:
742: .  w - the result

744:    Level: intermediate

746:    Notes:
747:     w cannot be either x or y, but x and y can be the same
748:     The implementation is optimzed for alpha of -1.0, 0.0, and 1.0

750: .seealso: `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
751: @*/
752: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
753: {
762:   VecCheckSameSize(x, 3, y, 4);
763:   VecCheckSameSize(x, 3, w, 1);
767:   VecSetErrorIfLocked(w, 1);

769:   VecLockReadPush(x);
770:   VecLockReadPush(y);
771:   if (alpha == (PetscScalar)0.0) {
772:     VecCopy(y, w);
773:   } else {
774:     PetscLogEventBegin(VEC_WAXPY, x, y, w, 0);
775:     PetscUseTypeMethod(w, waxpy, alpha, x, y);
776:     PetscLogEventEnd(VEC_WAXPY, x, y, w, 0);
777:     PetscObjectStateIncrease((PetscObject)w);
778:   }
779:   VecLockReadPop(x);
780:   VecLockReadPop(y);
781:   return 0;
782: }

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

787:    Not Collective

789:    Input Parameters:
790: +  x - vector to insert in
791: .  ni - number of elements to add
792: .  ix - indices where to add
793: .  y - array of values
794: -  iora - either INSERT_VALUES or ADD_VALUES, where
795:    ADD_VALUES adds values to any existing entries, and
796:    INSERT_VALUES replaces existing entries with new values

798:    Notes:
799:    VecSetValues() sets x[ix[i]] = y[i], for i=0,...,ni-1.

801:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
802:    options cannot be mixed without intervening calls to the assembly
803:    routines.

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

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

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

816:    Level: beginner

818: .seealso: `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
819:           `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
820: @*/
821: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
822: {
825:   if (!ni) return 0;

830:   PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0);
831:   PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
832:   PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0);
833:   PetscObjectStateIncrease((PetscObject)x);
834:   return 0;
835: }

837: /*@C
838:    VecGetValues - Gets values from certain locations of a vector. Currently
839:           can only get values on the same processor

841:     Not Collective

843:    Input Parameters:
844: +  x - vector to get values from
845: .  ni - number of elements to get
846: -  ix - indices where to get them from (in global 1d numbering)

848:    Output Parameter:
849: .   y - array of values

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

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

856:    VecAssemblyBegin() and VecAssemblyEnd()  MUST be called before calling this

858:    VecGetValues() 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.

864:    Level: beginner

866: .seealso: `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
867: @*/
868: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
869: {
871:   if (!ni) return 0;
875:   PetscUseTypeMethod(x, getvalues, ni, ix, y);
876:   return 0;
877: }

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

882:    Not Collective

884:    Input Parameters:
885: +  x - vector to insert in
886: .  ni - number of blocks to add
887: .  ix - indices where to add in block count, rather than element count
888: .  y - array of values
889: -  iora - either INSERT_VALUES or ADD_VALUES, where
890:    ADD_VALUES adds values to any existing entries, and
891:    INSERT_VALUES replaces existing entries with new values

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

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

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

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

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

911:    Level: intermediate

913: .seealso: `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
914:           `VecSetValues()`
915: @*/
916: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
917: {
920:   if (!ni) return 0;

925:   PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0);
926:   PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
927:   PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0);
928:   PetscObjectStateIncrease((PetscObject)x);
929:   return 0;
930: }

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

936:    Not Collective

938:    Input Parameters:
939: +  x - vector to insert in
940: .  ni - number of elements to add
941: .  ix - indices where to add
942: .  y - array of values
943: -  iora - either INSERT_VALUES or ADD_VALUES, where
944:    ADD_VALUES adds values to any existing entries, and
945:    INSERT_VALUES replaces existing entries with new values

947:    Level: intermediate

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

952:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
953:    options cannot be mixed without intervening calls to the assembly
954:    routines.

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

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

961: .seealso: `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
962:           `VecSetValuesBlockedLocal()`
963: @*/
964: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
965: {
966:   PetscInt lixp[128], *lix = lixp;

970:   if (!ni) return 0;

975:   PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0);
976:   if (!x->ops->setvalueslocal) {
977:     if (x->map->mapping) {
978:       if (ni > 128) PetscMalloc1(ni, &lix);
979:       ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix);
980:       PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
981:       if (ni > 128) PetscFree(lix);
982:     } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
983:   } else PetscUseTypeMethod(x, setvalueslocal, ni, ix, y, iora);
984:   PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0);
985:   PetscObjectStateIncrease((PetscObject)x);
986:   return 0;
987: }

989: /*@
990:    VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
991:    using a local ordering of the nodes.

993:    Not Collective

995:    Input Parameters:
996: +  x - vector to insert in
997: .  ni - number of blocks to add
998: .  ix - indices where to add in block count, not element count
999: .  y - array of values
1000: -  iora - either INSERT_VALUES or ADD_VALUES, where
1001:    ADD_VALUES adds values to any existing entries, and
1002:    INSERT_VALUES replaces existing entries with new values

1004:    Level: intermediate

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

1010:    Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1011:    options cannot be mixed without intervening calls to the assembly
1012:    routines.

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

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

1019: .seealso: `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1020:           `VecSetLocalToGlobalMapping()`
1021: @*/
1022: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1023: {
1024:   PetscInt lixp[128], *lix = lixp;

1028:   if (!ni) return 0;
1032:   PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0);
1033:   if (x->map->mapping) {
1034:     if (ni > 128) PetscMalloc1(ni, &lix);
1035:     ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix);
1036:     PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1037:     if (ni > 128) PetscFree(lix);
1038:   } else {
1039:     PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1040:   }
1041:   PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0);
1042:   PetscObjectStateIncrease((PetscObject)x);
1043:   return 0;
1044: }

1046: /*@
1047:    VecMTDot - Computes indefinite vector multiple dot products.
1048:    That is, it does NOT use the complex conjugate.

1050:    Collective on Vec

1052:    Input Parameters:
1053: +  x - one vector
1054: .  nv - number of vectors
1055: -  y - array of vectors.  Note that vectors are pointers

1057:    Output Parameter:
1058: .  val - array of the dot products

1060:    Notes for Users of Complex Numbers:
1061:    For complex vectors, VecMTDot() computes the indefinite form
1062: $      val = (x,y) = y^T x,
1063:    where y^T denotes the transpose of y.

1065:    Use VecMDot() for the inner product
1066: $      val = (x,y) = y^H x,
1067:    where y^H denotes the conjugate transpose of y.

1069:    Level: intermediate

1071: .seealso: `VecMDot()`, `VecTDot()`
1072: @*/
1073: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1074: {
1078:   if (!nv) return 0;
1080:   for (PetscInt i = 0; i < nv; ++i) {
1084:     VecCheckSameSize(x, 1, y[i], 3);
1085:     VecLockReadPush(y[i]);
1086:   }

1089:   VecLockReadPush(x);
1090:   PetscLogEventBegin(VEC_MTDot, x, *y, 0, 0);
1091:   PetscUseTypeMethod(x, mtdot, nv, y, val);
1092:   PetscLogEventEnd(VEC_MTDot, x, *y, 0, 0);
1093:   VecLockReadPop(x);
1094:   for (PetscInt i = 0; i < nv; ++i) VecLockReadPop(y[i]);
1095:   return 0;
1096: }

1098: /*@
1099:    VecMDot - Computes vector multiple dot products.

1101:    Collective on Vec

1103:    Input Parameters:
1104: +  x - one vector
1105: .  nv - number of vectors
1106: -  y - array of vectors.

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

1111:    Notes for Users of Complex Numbers:
1112:    For complex vectors, VecMDot() computes
1113: $     val = (x,y) = y^H x,
1114:    where y^H denotes the conjugate transpose of y.

1116:    Use VecMTDot() for the indefinite form
1117: $     val = (x,y) = y^T x,
1118:    where y^T denotes the transpose of y.

1120:    Level: intermediate

1122: .seealso: `VecMTDot()`, `VecDot()`
1123: @*/
1124: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1125: {
1129:   if (!nv) return 0;
1131:   for (PetscInt i = 0; i < nv; ++i) {
1135:     VecCheckSameSize(x, 1, y[i], 3);
1136:     VecLockReadPush(y[i]);
1137:   }

1140:   VecLockReadPush(x);
1141:   PetscLogEventBegin(VEC_MDot, x, *y, 0, 0);
1142:   PetscUseTypeMethod(x, mdot, nv, y, val);
1143:   PetscLogEventEnd(VEC_MDot, x, *y, 0, 0);
1144:   VecLockReadPop(x);
1145:   for (PetscInt i = 0; i < nv; ++i) VecLockReadPop(y[i]);
1146:   return 0;
1147: }

1149: /*@
1150:    VecMAXPY - Computes y = y + sum alpha[i] x[i]

1152:    Logically Collective on Vec

1154:    Input Parameters:
1155: +  nv - number of scalars and x-vectors
1156: .  alpha - array of scalars
1157: .  y - one vector
1158: -  x - array of vectors

1160:    Level: intermediate

1162:    Notes:
1163:     y cannot be any of the x vectors

1165: .seealso: `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1166: @*/
1167: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1168: {
1171:   VecSetErrorIfLocked(y, 1);
1173:   if (nv) {
1174:     PetscInt zeros = 0;

1178:     for (PetscInt i = 0; i < nv; ++i) {
1183:       VecCheckSameSize(y, 1, x[i], 4);
1185:       VecLockReadPush(x[i]);
1186:       zeros += alpha[i] == (PetscScalar)0.0;
1187:     }

1189:     if (zeros < nv) {
1190:       PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0);
1191:       PetscUseTypeMethod(y, maxpy, nv, alpha, x);
1192:       PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0);
1193:       PetscObjectStateIncrease((PetscObject)y);
1194:     }

1196:     for (PetscInt i = 0; i < nv; ++i) VecLockReadPop(x[i]);
1197:   }
1198:   return 0;
1199: }

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

1206:    Collective on X

1208:    Input Parameters:
1209: +  nx   - number of vectors to be concatenated
1210: -  X    - array containing the vectors to be concatenated in the order of concatenation

1212:    Output Parameters:
1213: +  Y    - concatenated vector
1214: -  x_is - array of index sets corresponding to the concatenated components of Y (NULL if not needed)

1216:    Notes:
1217:    Concatenation is similar to the functionality of a VecNest object; they both represent combination of
1218:    different vector spaces. However, concatenated vectors do not store any information about their
1219:    sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1220:    manipulation of data in the concatenated vector that corresponds to the original components at creation.

1222:    This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1223:    has to operate on combined vector spaces and cannot utilize VecNest objects due to incompatibility with
1224:    bound projections.

1226:    Level: advanced

1228: .seealso: `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1229: @*/
1230: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1231: {
1232:   MPI_Comm comm;
1233:   VecType  vec_type;
1234:   Vec      Ytmp, Xtmp;
1235:   IS      *is_tmp;
1236:   PetscInt i, shift = 0, Xnl, Xng, Xbegin;


1243:   if ((*X)->ops->concatenate) {
1244:     /* use the dedicated concatenation function if available */
1245:     (*(*X)->ops->concatenate)(nx, X, Y, x_is);
1246:   } else {
1247:     /* loop over vectors and start creating IS */
1248:     comm = PetscObjectComm((PetscObject)(*X));
1249:     VecGetType(*X, &vec_type);
1250:     PetscMalloc1(nx, &is_tmp);
1251:     for (i = 0; i < nx; i++) {
1252:       VecGetSize(X[i], &Xng);
1253:       VecGetLocalSize(X[i], &Xnl);
1254:       VecGetOwnershipRange(X[i], &Xbegin, NULL);
1255:       ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]);
1256:       shift += Xng;
1257:     }
1258:     /* create the concatenated vector */
1259:     VecCreate(comm, &Ytmp);
1260:     VecSetType(Ytmp, vec_type);
1261:     VecSetSizes(Ytmp, PETSC_DECIDE, shift);
1262:     VecSetUp(Ytmp);
1263:     /* copy data from X array to Y and return */
1264:     for (i = 0; i < nx; i++) {
1265:       VecGetSubVector(Ytmp, is_tmp[i], &Xtmp);
1266:       VecCopy(X[i], Xtmp);
1267:       VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp);
1268:     }
1269:     *Y = Ytmp;
1270:     if (x_is) {
1271:       *x_is = is_tmp;
1272:     } else {
1273:       for (i = 0; i < nx; i++) ISDestroy(&is_tmp[i]);
1274:       PetscFree(is_tmp);
1275:     }
1276:   }
1277:   return 0;
1278: }

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

1283:     Input Parameters:
1284: +   X - the original vector
1285: -   is - the index set of the subvector

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

1292: */
1293: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1294: {
1295:   PetscInt  gstart, gend, lstart;
1296:   PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1297:   PetscInt  n, N, ibs, vbs, bs = -1;

1299:   ISGetLocalSize(is, &n);
1300:   ISGetSize(is, &N);
1301:   ISGetBlockSize(is, &ibs);
1302:   VecGetBlockSize(X, &vbs);
1303:   VecGetOwnershipRange(X, &gstart, &gend);
1304:   ISContiguousLocal(is, gstart, gend, &lstart, &red[0]);
1305:   /* block size is given by IS if ibs > 1; otherwise, check the vector */
1306:   if (ibs > 1) {
1307:     MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is));
1308:     bs = ibs;
1309:   } else {
1310:     if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1311:     MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is));
1312:     if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1313:   }

1315:   *contig    = red[0];
1316:   *start     = lstart;
1317:   *blocksize = bs;
1318:   return 0;
1319: }

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

1323:     Input Parameters:
1324: +   X - the original vector
1325: .   is - the index set of the subvector
1326: -   bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()

1328:     Output Parameters:
1329: .   Z  - the subvector, which will compose the VecScatter context on output
1330: */
1331: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1332: {
1333:   PetscInt   n, N;
1334:   VecScatter vscat;
1335:   Vec        Y;

1337:   ISGetLocalSize(is, &n);
1338:   ISGetSize(is, &N);
1339:   VecCreate(PetscObjectComm((PetscObject)is), &Y);
1340:   VecSetSizes(Y, n, N);
1341:   VecSetBlockSize(Y, bs);
1342:   VecSetType(Y, ((PetscObject)X)->type_name);
1343:   VecScatterCreate(X, is, Y, NULL, &vscat);
1344:   VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD);
1345:   VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD);
1346:   PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat);
1347:   VecScatterDestroy(&vscat);
1348:   *Z = Y;
1349:   return 0;
1350: }

1352: /*@
1353:    VecGetSubVector - Gets a vector representing part of another vector

1355:    Collective on X and IS

1357:    Input Parameters:
1358: + X - vector from which to extract a subvector
1359: - is - index set representing portion of X to extract

1361:    Output Parameter:
1362: . Y - subvector corresponding to is

1364:    Level: advanced

1366:    Notes:
1367:    The subvector Y should be returned with VecRestoreSubVector().
1368:    X and is must be defined on the same communicator

1370:    This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1371:    modifying the subvector.  Other non-overlapping subvectors can still be obtained from X using this function.
1372:    The resulting subvector inherits the block size from the IS if greater than one. Otherwise, the block size is guessed from the block size of the original vec.

1374: .seealso: `MatCreateSubMatrix()`
1375: @*/
1376: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1377: {
1378:   Vec Z;

1384:   if (X->ops->getsubvector) {
1385:     PetscUseTypeMethod(X, getsubvector, is, &Z);
1386:   } else { /* Default implementation currently does no caching */
1387:     PetscBool contig;
1388:     PetscInt  n, N, start, bs;

1390:     ISGetLocalSize(is, &n);
1391:     ISGetSize(is, &N);
1392:     VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs);
1393:     if (contig) { /* We can do a no-copy implementation */
1394:       const PetscScalar *x;
1395:       PetscInt           state = 0;
1396:       PetscBool          isstd, iscuda, iship;

1398:       PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, "");
1399:       PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, "");
1400:       PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, "");
1401:       if (iscuda) {
1402: #if defined(PETSC_HAVE_CUDA)
1403:         const PetscScalar *x_d;
1404:         PetscMPIInt        size;
1405:         PetscOffloadMask   flg;

1407:         VecCUDAGetArrays_Private(X, &x, &x_d, &flg);
1410:         if (x) x += start;
1411:         if (x_d) x_d += start;
1412:         MPI_Comm_size(PetscObjectComm((PetscObject)X), &size);
1413:         if (size == 1) {
1414:           VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z);
1415:         } else {
1416:           VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z);
1417:         }
1418:         Z->offloadmask = flg;
1419: #endif
1420:       } else if (iship) {
1421: #if defined(PETSC_HAVE_HIP)
1422:         const PetscScalar *x_d;
1423:         PetscMPIInt        size;
1424:         PetscOffloadMask   flg;

1426:         VecHIPGetArrays_Private(X, &x, &x_d, &flg);
1429:         if (x) x += start;
1430:         if (x_d) x_d += start;
1431:         MPI_Comm_size(PetscObjectComm((PetscObject)X), &size);
1432:         if (size == 1) {
1433:           VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z);
1434:         } else {
1435:           VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z);
1436:         }
1437:         Z->offloadmask = flg;
1438: #endif
1439:       } else if (isstd) {
1440:         PetscMPIInt size;

1442:         MPI_Comm_size(PetscObjectComm((PetscObject)X), &size);
1443:         VecGetArrayRead(X, &x);
1444:         if (x) x += start;
1445:         if (size == 1) {
1446:           VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z);
1447:         } else {
1448:           VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z);
1449:         }
1450:         VecRestoreArrayRead(X, &x);
1451:       } else { /* default implementation: use place array */
1452:         VecGetArrayRead(X, &x);
1453:         VecCreate(PetscObjectComm((PetscObject)X), &Z);
1454:         VecSetType(Z, ((PetscObject)X)->type_name);
1455:         VecSetSizes(Z, n, N);
1456:         VecSetBlockSize(Z, bs);
1457:         VecPlaceArray(Z, x ? x + start : NULL);
1458:         VecRestoreArrayRead(X, &x);
1459:       }

1461:       /* this is relevant only in debug mode */
1462:       VecLockGet(X, &state);
1463:       if (state) VecLockReadPush(Z);
1464:       Z->ops->placearray   = NULL;
1465:       Z->ops->replacearray = NULL;
1466:     } else { /* Have to create a scatter and do a copy */
1467:       VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z);
1468:     }
1469:   }
1470:   /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1471:   if (VecGetSubVectorSavedStateId < 0) PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);
1472:   PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1);
1473:   *Y = Z;
1474:   return 0;
1475: }

1477: /*@
1478:    VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()

1480:    Collective on IS

1482:    Input Parameters:
1483: + X - vector from which subvector was obtained
1484: . is - index set representing the subset of X
1485: - Y - subvector being restored

1487:    Level: advanced

1489: .seealso: `VecGetSubVector()`
1490: @*/
1491: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1492: {
1493:   PETSC_UNUSED PetscObjectState dummystate = 0;
1494:   PetscBool                     unchanged;


1502:   if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1503:   else {
1504:     PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged);
1505:     if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1506:       VecScatter scatter;
1507:       PetscInt   state;

1509:       VecLockGet(X, &state);

1512:       PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter);
1513:       if (scatter) {
1514:         VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE);
1515:         VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE);
1516:       } else {
1517:         PetscBool iscuda, iship;
1518:         PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, "");
1519:         PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, "");

1521:         if (iscuda) {
1522: #if defined(PETSC_HAVE_CUDA)
1523:           PetscOffloadMask ymask = (*Y)->offloadmask;

1525:           /* The offloadmask of X dictates where to move memory
1526:               If X GPU data is valid, then move Y data on GPU if needed
1527:               Otherwise, move back to the CPU */
1528:           switch (X->offloadmask) {
1529:           case PETSC_OFFLOAD_BOTH:
1530:             if (ymask == PETSC_OFFLOAD_CPU) {
1531:               VecCUDAResetArray(*Y);
1532:             } else if (ymask == PETSC_OFFLOAD_GPU) {
1533:               X->offloadmask = PETSC_OFFLOAD_GPU;
1534:             }
1535:             break;
1536:           case PETSC_OFFLOAD_GPU:
1537:             if (ymask == PETSC_OFFLOAD_CPU) VecCUDAResetArray(*Y);
1538:             break;
1539:           case PETSC_OFFLOAD_CPU:
1540:             if (ymask == PETSC_OFFLOAD_GPU) VecResetArray(*Y);
1541:             break;
1542:           case PETSC_OFFLOAD_UNALLOCATED:
1543:           case PETSC_OFFLOAD_KOKKOS:
1544:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1545:           }
1546: #endif
1547:         } else if (iship) {
1548: #if defined(PETSC_HAVE_HIP)
1549:           PetscOffloadMask ymask = (*Y)->offloadmask;

1551:           /* The offloadmask of X dictates where to move memory
1552:               If X GPU data is valid, then move Y data on GPU if needed
1553:               Otherwise, move back to the CPU */
1554:           switch (X->offloadmask) {
1555:           case PETSC_OFFLOAD_BOTH:
1556:             if (ymask == PETSC_OFFLOAD_CPU) {
1557:               VecHIPResetArray(*Y);
1558:             } else if (ymask == PETSC_OFFLOAD_GPU) {
1559:               X->offloadmask = PETSC_OFFLOAD_GPU;
1560:             }
1561:             break;
1562:           case PETSC_OFFLOAD_GPU:
1563:             if (ymask == PETSC_OFFLOAD_CPU) VecHIPResetArray(*Y);
1564:             break;
1565:           case PETSC_OFFLOAD_CPU:
1566:             if (ymask == PETSC_OFFLOAD_GPU) VecResetArray(*Y);
1567:             break;
1568:           case PETSC_OFFLOAD_UNALLOCATED:
1569:           case PETSC_OFFLOAD_KOKKOS:
1570:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1571:           }
1572: #endif
1573:         } else {
1574:           /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1575:           VecResetArray(*Y);
1576:         }
1577:         PetscObjectStateIncrease((PetscObject)X);
1578:       }
1579:     }
1580:   }
1581:   VecDestroy(Y);
1582:   return 0;
1583: }

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

1589:    Not collective.

1591:    Input parameter:
1592: .  v - The vector for which the local vector is desired.

1594:    Output parameter:
1595: .  w - Upon exit this contains the local vector.

1597:    Level: beginner

1599: .seealso: `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1600: @*/
1601: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1602: {
1603:   PetscMPIInt size;

1607:   MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
1608:   if (size == 1) VecDuplicate(v, w);
1609:   else if (v->ops->createlocalvector) PetscUseTypeMethod(v, createlocalvector, w);
1610:   else {
1611:     VecType  type;
1612:     PetscInt n;

1614:     VecCreate(PETSC_COMM_SELF, w);
1615:     VecGetLocalSize(v, &n);
1616:     VecSetSizes(*w, n, n);
1617:     VecGetBlockSize(v, &n);
1618:     VecSetBlockSize(*w, n);
1619:     VecGetType(v, &type);
1620:     VecSetType(*w, type);
1621:   }
1622:   return 0;
1623: }

1625: /*@
1626:    VecGetLocalVectorRead - Maps the local portion of a vector into a
1627:    vector.  You must call VecRestoreLocalVectorRead() when the local
1628:    vector is no longer needed.

1630:    Not collective.

1632:    Input parameter:
1633: .  v - The vector for which the local vector is desired.

1635:    Output parameter:
1636: .  w - Upon exit this contains the local vector.

1638:    Level: beginner

1640:    Notes:
1641:    This function is similar to VecGetArrayRead() which maps the local
1642:    portion into a raw pointer.  VecGetLocalVectorRead() is usually
1643:    almost as efficient as VecGetArrayRead() but in certain circumstances
1644:    VecGetLocalVectorRead() can be much more efficient than
1645:    VecGetArrayRead().  This is because the construction of a contiguous
1646:    array representing the vector data required by VecGetArrayRead() can
1647:    be an expensive operation for certain vector types.  For example, for
1648:    GPU vectors VecGetArrayRead() requires that the data between device
1649:    and host is synchronized.

1651:    Unlike VecGetLocalVector(), this routine is not collective and
1652:    preserves cached information.

1654: .seealso: `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1655: @*/
1656: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1657: {
1660:   VecCheckSameLocalSize(v, 1, w, 2);
1661:   if (v->ops->getlocalvectorread) {
1662:     PetscUseTypeMethod(v, getlocalvectorread, w);
1663:   } else {
1664:     PetscScalar *a;

1666:     VecGetArrayRead(v, (const PetscScalar **)&a);
1667:     VecPlaceArray(w, a);
1668:   }
1669:   PetscObjectStateIncrease((PetscObject)w);
1670:   VecLockReadPush(v);
1671:   VecLockReadPush(w);
1672:   return 0;
1673: }

1675: /*@
1676:    VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1677:    previously mapped into a vector using VecGetLocalVectorRead().

1679:    Not collective.

1681:    Input parameter:
1682: +  v - The local portion of this vector was previously mapped into w using VecGetLocalVectorRead().
1683: -  w - The vector into which the local portion of v was mapped.

1685:    Level: beginner

1687: .seealso: `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1688: @*/
1689: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1690: {
1693:   if (v->ops->restorelocalvectorread) {
1694:     PetscUseTypeMethod(v, restorelocalvectorread, w);
1695:   } else {
1696:     const PetscScalar *a;

1698:     VecGetArrayRead(w, &a);
1699:     VecRestoreArrayRead(v, &a);
1700:     VecResetArray(w);
1701:   }
1702:   VecLockReadPop(v);
1703:   VecLockReadPop(w);
1704:   PetscObjectStateIncrease((PetscObject)w);
1705:   return 0;
1706: }

1708: /*@
1709:    VecGetLocalVector - Maps the local portion of a vector into a
1710:    vector.

1712:    Collective on v, not collective on w.

1714:    Input parameter:
1715: .  v - The vector for which the local vector is desired.

1717:    Output parameter:
1718: .  w - Upon exit this contains the local vector.

1720:    Level: beginner

1722:    Notes:
1723:    This function is similar to VecGetArray() which maps the local
1724:    portion into a raw pointer.  VecGetLocalVector() is usually about as
1725:    efficient as VecGetArray() but in certain circumstances
1726:    VecGetLocalVector() can be much more efficient than VecGetArray().
1727:    This is because the construction of a contiguous array representing
1728:    the vector data required by VecGetArray() can be an expensive
1729:    operation for certain vector types.  For example, for GPU vectors
1730:    VecGetArray() requires that the data between device and host is
1731:    synchronized.

1733: .seealso: `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1734: @*/
1735: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1736: {
1739:   VecCheckSameLocalSize(v, 1, w, 2);
1740:   if (v->ops->getlocalvector) {
1741:     PetscUseTypeMethod(v, getlocalvector, w);
1742:   } else {
1743:     PetscScalar *a;

1745:     VecGetArray(v, &a);
1746:     VecPlaceArray(w, a);
1747:   }
1748:   PetscObjectStateIncrease((PetscObject)w);
1749:   return 0;
1750: }

1752: /*@
1753:    VecRestoreLocalVector - Unmaps the local portion of a vector
1754:    previously mapped into a vector using VecGetLocalVector().

1756:    Logically collective.

1758:    Input parameter:
1759: +  v - The local portion of this vector was previously mapped into w using VecGetLocalVector().
1760: -  w - The vector into which the local portion of v was mapped.

1762:    Level: beginner

1764: .seealso: `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1765: @*/
1766: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1767: {
1770:   if (v->ops->restorelocalvector) {
1771:     PetscUseTypeMethod(v, restorelocalvector, w);
1772:   } else {
1773:     PetscScalar *a;
1774:     VecGetArray(w, &a);
1775:     VecRestoreArray(v, &a);
1776:     VecResetArray(w);
1777:   }
1778:   PetscObjectStateIncrease((PetscObject)w);
1779:   PetscObjectStateIncrease((PetscObject)v);
1780:   return 0;
1781: }

1783: /*@C
1784:    VecGetArray - Returns a pointer to a contiguous array that contains this
1785:    processor's portion of the vector data. For the standard PETSc
1786:    vectors, VecGetArray() returns a pointer to the local data array and
1787:    does not use any copies. If the underlying vector data is not stored
1788:    in a contiguous array this routine will copy the data to a contiguous
1789:    array and return a pointer to that. You MUST call VecRestoreArray()
1790:    when you no longer need access to the array.

1792:    Logically Collective on Vec

1794:    Input Parameter:
1795: .  x - the vector

1797:    Output Parameter:
1798: .  a - location to put pointer to the array

1800:    Fortran Note:
1801:    This routine is used differently from Fortran 77
1802: .vb
1803:     Vec         x
1804:     PetscScalar x_array(1)
1805:     PetscOffset i_x
1806:     PetscErrorCode ierr
1807:     call VecGetArray(x,x_array,i_x,ierr)

1809:     Access first local entry in vector with
1810:     value = x_array(i_x + 1)

1812:     ...... other code
1813:     call VecRestoreArray(x,x_array,i_x,ierr)
1814: .ve
1815:    For Fortran 90 see `VecGetArrayF90()`

1817:    See the Fortran chapter of the users manual and
1818:    petsc/src/snes/tutorials/ex5f.F for details.

1820:    Level: beginner

1822: .seealso: `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
1823:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
1824: @*/
1825: PetscErrorCode VecGetArray(Vec x, PetscScalar **a)
1826: {
1828:   VecSetErrorIfLocked(x, 1);
1829:   if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
1830:     PetscUseTypeMethod(x, getarray, a);
1831:   } else if (x->petscnative) { /* VECSTANDARD */
1832:     *a = *((PetscScalar **)x->data);
1833:   } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
1834:   return 0;
1835: }

1837: /*@C
1838:    VecRestoreArray - Restores a vector after VecGetArray() has been called.

1840:    Logically Collective on Vec

1842:    Input Parameters:
1843: +  x - the vector
1844: -  a - location of pointer to array obtained from VecGetArray()

1846:    Level: beginner

1848: .seealso: `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
1849:           `VecGetArrayPair()`, `VecRestoreArrayPair()`
1850: @*/
1851: PetscErrorCode VecRestoreArray(Vec x, PetscScalar **a)
1852: {
1855:   if (x->ops->restorearray) {
1856:     PetscUseTypeMethod(x, restorearray, a);
1858:   if (a) *a = NULL;
1859:   PetscObjectStateIncrease((PetscObject)x);
1860:   return 0;
1861: }
1862: /*@C
1863:    VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.

1865:    Not Collective

1867:    Input Parameter:
1868: .  x - the vector

1870:    Output Parameter:
1871: .  a - the array

1873:    Level: beginner

1875:    Notes:
1876:    The array must be returned using a matching call to VecRestoreArrayRead().

1878:    Unlike VecGetArray(), this routine is not collective and preserves cached information like vector norms.

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

1884: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
1885: @*/
1886: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar **a)
1887: {
1890:   if (x->ops->getarrayread) {
1891:     PetscUseTypeMethod(x, getarrayread, a);
1892:   } else if (x->ops->getarray) {
1893:     /* VECNEST, VECCUDA, VECKOKKOS etc */
1894:     PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
1895:   } else if (x->petscnative) {
1896:     /* VECSTANDARD */
1897:     *a = *((PetscScalar **)x->data);
1898:   } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
1899:   return 0;
1900: }

1902: /*@C
1903:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

1905:    Not Collective

1907:    Input Parameters:
1908: +  vec - the vector
1909: -  array - the array

1911:    Level: beginner

1913: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
1914: @*/
1915: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar **a)
1916: {
1919:   if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
1920:     /* nothing */
1921:   } else if (x->ops->restorearrayread) { /* VECNEST */
1922:     PetscUseTypeMethod(x, restorearrayread, a);
1923:   } else { /* No one? */
1924:     PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
1925:   }
1926:   if (a) *a = NULL;
1927:   return 0;
1928: }

1930: /*@C
1931:    VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contains this
1932:    processor's portion of the vector data. The values in this array are NOT valid, the routine calling this
1933:    routine is responsible for putting values into the array; any values it does not set will be invalid

1935:    Logically Collective on Vec

1937:    Input Parameter:
1938: .  x - the vector

1940:    Output Parameter:
1941: .  a - location to put pointer to the array

1943:    Level: intermediate

1945:    This is for vectors associate with GPUs, the vector is not copied up before giving access. If you need correct
1946:    values in the array use VecGetArray()

1948: .seealso: `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
1949:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
1950: @*/
1951: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar **a)
1952: {
1955:   VecSetErrorIfLocked(x, 1);
1956:   if (x->ops->getarraywrite) {
1957:     PetscUseTypeMethod(x, getarraywrite, a);
1958:   } else {
1959:     VecGetArray(x, a);
1960:   }
1961:   return 0;
1962: }

1964: /*@C
1965:    VecRestoreArrayWrite - Restores a vector after VecGetArrayWrite() has been called.

1967:    Logically Collective on Vec

1969:    Input Parameters:
1970: +  x - the vector
1971: -  a - location of pointer to array obtained from VecGetArray()

1973:    Level: beginner

1975: .seealso: `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
1976:           `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
1977: @*/
1978: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar **a)
1979: {
1982:   if (x->ops->restorearraywrite) {
1983:     PetscUseTypeMethod(x, restorearraywrite, a);
1984:   } else if (x->ops->restorearray) {
1985:     PetscUseTypeMethod(x, restorearray, a);
1986:   }
1987:   if (a) *a = NULL;
1988:   PetscObjectStateIncrease((PetscObject)x);
1989:   return 0;
1990: }

1992: /*@C
1993:    VecGetArrays - Returns a pointer to the arrays in a set of vectors
1994:    that were created by a call to VecDuplicateVecs().  You MUST call
1995:    VecRestoreArrays() when you no longer need access to the array.

1997:    Logically Collective on Vec

1999:    Input Parameters:
2000: +  x - the vectors
2001: -  n - the number of vectors

2003:    Output Parameter:
2004: .  a - location to put pointer to the array

2006:    Fortran Note:
2007:    This routine is not supported in Fortran.

2009:    Level: intermediate

2011: .seealso: `VecGetArray()`, `VecRestoreArrays()`
2012: @*/
2013: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2014: {
2015:   PetscInt      i;
2016:   PetscScalar **q;

2022:   PetscMalloc1(n, &q);
2023:   for (i = 0; i < n; ++i) VecGetArray(x[i], &q[i]);
2024:   *a = q;
2025:   return 0;
2026: }

2028: /*@C
2029:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
2030:    has been called.

2032:    Logically Collective on Vec

2034:    Input Parameters:
2035: +  x - the vector
2036: .  n - the number of vectors
2037: -  a - location of pointer to arrays obtained from VecGetArrays()

2039:    Notes:
2040:    For regular PETSc vectors this routine does not involve any copies. For
2041:    any special vectors that do not store local vector data in a contiguous
2042:    array, this routine will copy the data back into the underlying
2043:    vector data structure from the arrays obtained with VecGetArrays().

2045:    Fortran Note:
2046:    This routine is not supported in Fortran.

2048:    Level: intermediate

2050: .seealso: `VecGetArrays()`, `VecRestoreArray()`
2051: @*/
2052: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2053: {
2054:   PetscInt      i;
2055:   PetscScalar **q = *a;


2061:   for (i = 0; i < n; ++i) VecRestoreArray(x[i], &q[i]);
2062:   PetscFree(q);
2063:   return 0;
2064: }

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

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

2074:    Logically Collective on Vec

2076:    Input Parameter:
2077: .  x - the vector

2079:    Output Parameters:
2080: +  a - location to put pointer to the array
2081: -  mtype - memory type of the array

2083:    Level: beginner

2085: .seealso: `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`,
2086:           `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2087: @*/
2088: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2089: {
2094:   VecSetErrorIfLocked(x, 1);
2095:   if (x->ops->getarrayandmemtype) {
2096:     /* VECCUDA, VECKOKKOS etc */
2097:     PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2098:   } else {
2099:     /* VECSTANDARD, VECNEST, VECVIENNACL */
2100:     VecGetArray(x, a);
2101:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2102:   }
2103:   return 0;
2104: }

2106: /*@C
2107:    VecRestoreArrayAndMemType - Restores a vector after VecGetArrayAndMemType() has been called.

2109:    Logically Collective on Vec

2111:    Input Parameters:
2112: +  x - the vector
2113: -  a - location of pointer to array obtained from VecGetArrayAndMemType()

2115:    Level: beginner

2117: .seealso: `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`,
2118:           `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2119: @*/
2120: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar **a)
2121: {
2125:   if (x->ops->restorearrayandmemtype) {
2126:     /* VECCUDA, VECKOKKOS etc */
2127:     PetscUseTypeMethod(x, restorearrayandmemtype, a);
2128:   } else {
2129:     /* VECNEST, VECVIENNACL */
2130:     VecRestoreArray(x, a);
2131:   } /* VECSTANDARD does nothing */
2132:   if (a) *a = NULL;
2133:   PetscObjectStateIncrease((PetscObject)x);
2134:   return 0;
2135: }

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

2140:    Not Collective

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

2145:    Output Parameters:
2146: +  a - the array
2147: -  mtype - memory type of the array

2149:    Level: beginner

2151:    Notes:
2152:    The array must be returned using a matching call to VecRestoreArrayReadAndMemType().

2154: .seealso: `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayAndMemType()`
2155: @*/
2156: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar **a, PetscMemType *mtype)
2157: {
2162:   if (x->ops->getarrayreadandmemtype) {
2163:     /* VECCUDA/VECHIP though they are also petscnative */
2164:     PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2165:   } else if (x->ops->getarrayandmemtype) {
2166:     /* VECKOKKOS */
2167:     PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2168:   } else {
2169:     VecGetArrayRead(x, a);
2170:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2171:   }
2172:   return 0;
2173: }

2175: /*@C
2176:    VecRestoreArrayReadAndMemType - Restore array obtained with VecGetArrayReadAndMemType()

2178:    Not Collective

2180:    Input Parameters:
2181: +  vec - the vector
2182: -  array - the array

2184:    Level: beginner

2186: .seealso: `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2187: @*/
2188: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar **a)
2189: {
2193:   if (x->ops->restorearrayreadandmemtype) {
2194:     /* VECCUDA/VECHIP */
2195:     PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2196:   } else if (!x->petscnative) {
2197:     /* VECNEST */
2198:     VecRestoreArrayRead(x, a);
2199:   }
2200:   if (a) *a = NULL;
2201:   return 0;
2202: }

2204: /*@C
2205:    VecGetArrayWriteAndMemType - Like VecGetArrayWrite(), but if this is a device vector it will aways return
2206:     a device pointer to the device memory that contains this processor's portion of the vector data.

2208:    Not Collective

2210:    Input Parameter:
2211: .  x - the vector

2213:    Output Parameters:
2214: +  a - the array
2215: -  mtype - memory type of the array

2217:    Level: beginner

2219:    Notes:
2220:    The array must be returned using a matching call to VecRestoreArrayWriteAndMemType(), where it will label the device memory as most recent.

2222: .seealso: `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2223: @*/
2224: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2225: {
2228:   VecSetErrorIfLocked(x, 1);
2231:   if (x->ops->getarraywriteandmemtype) {
2232:     /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2233:     PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2234:   } else if (x->ops->getarrayandmemtype) {
2235:     VecGetArrayAndMemType(x, a, mtype);
2236:   } else {
2237:     /* VECNEST, VECVIENNACL */
2238:     VecGetArrayWrite(x, a);
2239:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2240:   }
2241:   return 0;
2242: }

2244: /*@C
2245:    VecRestoreArrayWriteAndMemType - Restore array obtained with VecGetArrayWriteAndMemType()

2247:    Not Collective

2249:    Input Parameters:
2250: +  vec - the vector
2251: -  array - the array

2253:    Level: beginner

2255: .seealso: `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2256: @*/
2257: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar **a)
2258: {
2261:   VecSetErrorIfLocked(x, 1);
2263:   if (x->ops->restorearraywriteandmemtype) {
2264:     /* VECCUDA/VECHIP */
2265:     PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2266:     PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2267:   } else if (x->ops->restorearrayandmemtype) {
2268:     VecRestoreArrayAndMemType(x, a);
2269:   } else {
2270:     VecRestoreArray(x, a);
2271:   }
2272:   if (a) *a = NULL;
2273:   return 0;
2274: }

2276: /*@
2277:    VecPlaceArray - Allows one to replace the array in a vector with an
2278:    array provided by the user. This is useful to avoid copying an array
2279:    into a vector.

2281:    Not Collective

2283:    Input Parameters:
2284: +  vec - the vector
2285: -  array - the array

2287:    Notes:
2288:    You can return to the original array with a call to `VecResetArray()`. `vec` does not take
2289:    ownership of `array` in any way. The user must free `array` themselves but be careful not to
2290:    do so before the vector has either been destroyed, had its original array restored with
2291:    `VecResetArray()` or permanently replaced with `VecReplaceArray()`.

2293:    Level: developer

2295: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`

2297: @*/
2298: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2299: {
2303:   PetscUseTypeMethod(vec, placearray, array);
2304:   PetscObjectStateIncrease((PetscObject)vec);
2305:   return 0;
2306: }

2308: /*@C
2309:    VecReplaceArray - Allows one to replace the array in a vector with an
2310:    array provided by the user. This is useful to avoid copying an array
2311:    into a vector.

2313:    Not Collective

2315:    Input Parameters:
2316: +  vec - the vector
2317: -  array - the array

2319:    Notes:
2320:    This permanently replaces the array and frees the memory associated
2321:    with the old array.

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

2326:    Not supported from Fortran

2328:    Level: developer

2330: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`

2332: @*/
2333: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2334: {
2337:   PetscUseTypeMethod(vec, replacearray, array);
2338:   PetscObjectStateIncrease((PetscObject)vec);
2339:   return 0;
2340: }

2342: /*@C
2343:    VecCUDAGetArray - Provides access to the CUDA buffer inside a vector.

2345:    This function has semantics similar to VecGetArray():  the pointer
2346:    returned by this function points to a consistent view of the vector
2347:    data.  This may involve a copy operation of data from the host to the
2348:    device if the data on the device is out of date.  If the device
2349:    memory hasn't been allocated previously it will be allocated as part
2350:    of this function call.  VecCUDAGetArray() assumes that
2351:    the user will modify the vector data.  This is similar to
2352:    intent(inout) in fortran.

2354:    The CUDA device pointer has to be released by calling
2355:    VecCUDARestoreArray().  Upon restoring the vector data
2356:    the data on the host will be marked as out of date.  A subsequent
2357:    access of the host data will thus incur a data transfer from the
2358:    device to the host.

2360:    Input Parameter:
2361: .  v - the vector

2363:    Output Parameter:
2364: .  a - the CUDA device pointer

2366:    Fortran note:
2367:    This function is not currently available from Fortran.

2369:    Level: intermediate

2371: .seealso: `VecCUDARestoreArray()`, `VecCUDAGetArrayRead()`, `VecCUDAGetArrayWrite()`, `VecGetArray()`, `VecGetArrayRead()`
2372: @*/
2373: PETSC_EXTERN PetscErrorCode VecCUDAGetArray(Vec v, PetscScalar **a)
2374: {
2376: #if defined(PETSC_HAVE_CUDA)
2377:   {
2378:     VecCUDACopyToGPU(v);
2379:     *a = ((Vec_CUDA *)v->spptr)->GPUarray;
2380:   }
2381: #endif
2382:   return 0;
2383: }

2385: /*@C
2386:    VecCUDARestoreArray - Restore a CUDA device pointer previously acquired with VecCUDAGetArray().

2388:    This marks the host data as out of date.  Subsequent access to the
2389:    vector data on the host side with for instance VecGetArray() incurs a
2390:    data transfer.

2392:    Input Parameters:
2393: +  v - the vector
2394: -  a - the CUDA device pointer.  This pointer is invalid after
2395:        VecCUDARestoreArray() returns.

2397:    Fortran note:
2398:    This function is not currently available from Fortran.

2400:    Level: intermediate

2402: .seealso: `VecCUDAGetArray()`, `VecCUDAGetArrayRead()`, `VecCUDAGetArrayWrite()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`
2403: @*/
2404: PETSC_EXTERN PetscErrorCode VecCUDARestoreArray(Vec v, PetscScalar **a)
2405: {
2407: #if defined(PETSC_HAVE_CUDA)
2408:   v->offloadmask = PETSC_OFFLOAD_GPU;
2409: #endif
2410:   PetscObjectStateIncrease((PetscObject)v);
2411:   return 0;
2412: }

2414: /*@C
2415:    VecCUDAGetArrayRead - Provides read access to the CUDA buffer inside a vector.

2417:    This function is analogous to VecGetArrayRead():  The pointer
2418:    returned by this function points to a consistent view of the vector
2419:    data.  This may involve a copy operation of data from the host to the
2420:    device if the data on the device is out of date.  If the device
2421:    memory hasn't been allocated previously it will be allocated as part
2422:    of this function call.  VecCUDAGetArrayRead() assumes that the
2423:    user will not modify the vector data.  This is analgogous to
2424:    intent(in) in Fortran.

2426:    The CUDA device pointer has to be released by calling
2427:    VecCUDARestoreArrayRead().  If the data on the host side was
2428:    previously up to date it will remain so, i.e. data on both the device
2429:    and the host is up to date.  Accessing data on the host side does not
2430:    incur a device to host data transfer.

2432:    Input Parameter:
2433: .  v - the vector

2435:    Output Parameter:
2436: .  a - the CUDA pointer.

2438:    Fortran note:
2439:    This function is not currently available from Fortran.

2441:    Level: intermediate

2443: .seealso: `VecCUDARestoreArrayRead()`, `VecCUDAGetArray()`, `VecCUDAGetArrayWrite()`, `VecGetArray()`, `VecGetArrayRead()`
2444: @*/
2445: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayRead(Vec v, const PetscScalar **a)
2446: {
2447:   VecCUDAGetArray(v, (PetscScalar **)a);
2448:   return 0;
2449: }

2451: /*@C
2452:    VecCUDARestoreArrayRead - Restore a CUDA device pointer previously acquired with VecCUDAGetArrayRead().

2454:    If the data on the host side was previously up to date it will remain
2455:    so, i.e. data on both the device and the host is up to date.
2456:    Accessing data on the host side e.g. with VecGetArray() does not
2457:    incur a device to host data transfer.

2459:    Input Parameters:
2460: +  v - the vector
2461: -  a - the CUDA device pointer.  This pointer is invalid after
2462:        VecCUDARestoreArrayRead() returns.

2464:    Fortran note:
2465:    This function is not currently available from Fortran.

2467:    Level: intermediate

2469: .seealso: `VecCUDAGetArrayRead()`, `VecCUDAGetArrayWrite()`, `VecCUDAGetArray()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`
2470: @*/
2471: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayRead(Vec v, const PetscScalar **a)
2472: {
2474:   *a = NULL;
2475:   return 0;
2476: }

2478: /*@C
2479:    VecCUDAGetArrayWrite - Provides write access to the CUDA buffer inside a vector.

2481:    The data pointed to by the device pointer is uninitialized.  The user
2482:    may not read from this data.  Furthermore, the entire array needs to
2483:    be filled by the user to obtain well-defined behaviour.  The device
2484:    memory will be allocated by this function if it hasn't been allocated
2485:    previously.  This is analogous to intent(out) in Fortran.

2487:    The device pointer needs to be released with
2488:    VecCUDARestoreArrayWrite().  When the pointer is released the
2489:    host data of the vector is marked as out of data.  Subsequent access
2490:    of the host data with e.g. VecGetArray() incurs a device to host data
2491:    transfer.

2493:    Input Parameter:
2494: .  v - the vector

2496:    Output Parameter:
2497: .  a - the CUDA pointer

2499:    Fortran note:
2500:    This function is not currently available from Fortran.

2502:    Level: advanced

2504: .seealso: `VecCUDARestoreArrayWrite()`, `VecCUDAGetArray()`, `VecCUDAGetArrayRead()`, `VecCUDAGetArrayWrite()`, `VecGetArray()`, `VecGetArrayRead()`
2505: @*/
2506: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayWrite(Vec v, PetscScalar **a)
2507: {
2509: #if defined(PETSC_HAVE_CUDA)
2510:   {
2511:     VecCUDAAllocateCheck(v);
2512:     *a = ((Vec_CUDA *)v->spptr)->GPUarray;
2513:   }
2514: #endif
2515:   return 0;
2516: }

2518: /*@C
2519:    VecCUDARestoreArrayWrite - Restore a CUDA device pointer previously acquired with VecCUDAGetArrayWrite().

2521:    Data on the host will be marked as out of date.  Subsequent access of
2522:    the data on the host side e.g. with VecGetArray() will incur a device
2523:    to host data transfer.

2525:    Input Parameters:
2526: +  v - the vector
2527: -  a - the CUDA device pointer.  This pointer is invalid after
2528:        VecCUDARestoreArrayWrite() returns.

2530:    Fortran note:
2531:    This function is not currently available from Fortran.

2533:    Level: intermediate

2535: .seealso: `VecCUDAGetArrayWrite()`, `VecCUDAGetArray()`, `VecCUDAGetArrayRead()`, `VecCUDAGetArrayWrite()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`
2536: @*/
2537: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayWrite(Vec v, PetscScalar **a)
2538: {
2540: #if defined(PETSC_HAVE_CUDA)
2541:   v->offloadmask = PETSC_OFFLOAD_GPU;
2542:   if (a) *a = NULL;
2543: #endif
2544:   PetscObjectStateIncrease((PetscObject)v);
2545:   return 0;
2546: }

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

2553:    Not Collective

2555:    Input Parameters:
2556: +  vec - the vector
2557: -  array - the GPU array

2559:    Notes:
2560:    You can return to the original GPU array with a call to VecCUDAResetArray()
2561:    It is not possible to use VecCUDAPlaceArray() and VecPlaceArray() at the
2562:    same time on the same vector.

2564:    `vec` does not take ownership of `array` in any way. The user must free `array` themselves
2565:    but be careful not to do so before the vector has either been destroyed, had its original
2566:    array restored with `VecCUDAResetArray()` or permanently replaced with
2567:    `VecCUDAReplaceArray()`.

2569:    Level: developer

2571: .seealso: `VecPlaceArray()`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`, `VecCUDAResetArray()`, `VecCUDAReplaceArray()`

2573: @*/
2574: PetscErrorCode VecCUDAPlaceArray(Vec vin, const PetscScalar a[])
2575: {
2577: #if defined(PETSC_HAVE_CUDA)
2578:   VecCUDACopyToGPU(vin);
2580:   ((Vec_Seq *)vin->data)->unplacedarray = (PetscScalar *)((Vec_CUDA *)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2581:   ((Vec_CUDA *)vin->spptr)->GPUarray    = (PetscScalar *)a;
2582:   vin->offloadmask                      = PETSC_OFFLOAD_GPU;
2583: #endif
2584:   PetscObjectStateIncrease((PetscObject)vin);
2585:   return 0;
2586: }

2588: /*@C
2589:    VecCUDAReplaceArray - Allows one to replace the GPU array in a vector
2590:    with a GPU array provided by the user. This is useful to avoid copying
2591:    a GPU array into a vector.

2593:    Not Collective

2595:    Input Parameters:
2596: +  vec - the vector
2597: -  array - the GPU array

2599:    Notes:
2600:    This permanently replaces the GPU array and frees the memory associated
2601:    with the old GPU array.

2603:    The memory passed in CANNOT be freed by the user. It will be freed
2604:    when the vector is destroyed.

2606:    Not supported from Fortran

2608:    Level: developer

2610: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`, `VecCUDAResetArray()`, `VecCUDAPlaceArray()`, `VecReplaceArray()`

2612: @*/
2613: PetscErrorCode VecCUDAReplaceArray(Vec vin, const PetscScalar a[])
2614: {
2615: #if defined(PETSC_HAVE_CUDA)
2616: #endif

2619: #if defined(PETSC_HAVE_CUDA)
2620:   if (((Vec_CUDA *)vin->spptr)->GPUarray_allocated) cudaFree(((Vec_CUDA *)vin->spptr)->GPUarray_allocated);
2621:   ((Vec_CUDA *)vin->spptr)->GPUarray_allocated = ((Vec_CUDA *)vin->spptr)->GPUarray = (PetscScalar *)a;
2622:   vin->offloadmask                                                                  = PETSC_OFFLOAD_GPU;
2623: #endif
2624:   PetscObjectStateIncrease((PetscObject)vin);
2625:   return 0;
2626: }

2628: /*@C
2629:    VecCUDAResetArray - Resets a vector to use its default memory. Call this
2630:    after the use of VecCUDAPlaceArray().

2632:    Not Collective

2634:    Input Parameters:
2635: .  vec - the vector

2637:    Level: developer

2639: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecPlaceArray()`, `VecResetArray()`, `VecCUDAPlaceArray()`, `VecCUDAReplaceArray()`

2641: @*/
2642: PetscErrorCode VecCUDAResetArray(Vec vin)
2643: {
2645: #if defined(PETSC_HAVE_CUDA)
2646:   VecCUDACopyToGPU(vin);
2647:   ((Vec_CUDA *)vin->spptr)->GPUarray    = (PetscScalar *)((Vec_Seq *)vin->data)->unplacedarray;
2648:   ((Vec_Seq *)vin->data)->unplacedarray = 0;
2649:   vin->offloadmask                      = PETSC_OFFLOAD_GPU;
2650: #endif
2651:   PetscObjectStateIncrease((PetscObject)vin);
2652:   return 0;
2653: }

2655: /*@C
2656:    VecHIPGetArray - Provides access to the HIP buffer inside a vector.

2658:    This function has semantics similar to VecGetArray():  the pointer
2659:    returned by this function points to a consistent view of the vector
2660:    data.  This may involve a copy operation of data from the host to the
2661:    device if the data on the device is out of date.  If the device
2662:    memory hasn't been allocated previously it will be allocated as part
2663:    of this function call.  VecHIPGetArray() assumes that
2664:    the user will modify the vector data.  This is similar to
2665:    intent(inout) in fortran.

2667:    The HIP device pointer has to be released by calling
2668:    VecHIPRestoreArray().  Upon restoring the vector data
2669:    the data on the host will be marked as out of date.  A subsequent
2670:    access of the host data will thus incur a data transfer from the
2671:    device to the host.

2673:    Input Parameter:
2674: .  v - the vector

2676:    Output Parameter:
2677: .  a - the HIP device pointer

2679:    Fortran note:
2680:    This function is not currently available from Fortran.

2682:    Level: intermediate

2684: .seealso: `VecHIPRestoreArray()`, `VecHIPGetArrayRead()`, `VecHIPGetArrayWrite()`, `VecGetArray()`, `VecGetArrayRead()`
2685: @*/
2686: PETSC_EXTERN PetscErrorCode VecHIPGetArray(Vec v, PetscScalar **a)
2687: {
2689: #if defined(PETSC_HAVE_HIP)
2690:   *a = 0;
2691:   VecHIPCopyToGPU(v);
2692:   *a = ((Vec_HIP *)v->spptr)->GPUarray;
2693: #endif
2694:   return 0;
2695: }

2697: /*@C
2698:    VecHIPRestoreArray - Restore a HIP device pointer previously acquired with VecHIPGetArray().

2700:    This marks the host data as out of date.  Subsequent access to the
2701:    vector data on the host side with for instance VecGetArray() incurs a
2702:    data transfer.

2704:    Input Parameters:
2705: +  v - the vector
2706: -  a - the HIP device pointer.  This pointer is invalid after
2707:        VecHIPRestoreArray() returns.

2709:    Fortran note:
2710:    This function is not currently available from Fortran.

2712:    Level: intermediate

2714: .seealso: `VecHIPGetArray()`, `VecHIPGetArrayRead()`, `VecHIPGetArrayWrite()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`
2715: @*/
2716: PETSC_EXTERN PetscErrorCode VecHIPRestoreArray(Vec v, PetscScalar **a)
2717: {
2719: #if defined(PETSC_HAVE_HIP)
2720:   v->offloadmask = PETSC_OFFLOAD_GPU;
2721: #endif

2723:   PetscObjectStateIncrease((PetscObject)v);
2724:   return 0;
2725: }

2727: /*@C
2728:    VecHIPGetArrayRead - Provides read access to the HIP buffer inside a vector.

2730:    This function is analogous to VecGetArrayRead():  The pointer
2731:    returned by this function points to a consistent view of the vector
2732:    data.  This may involve a copy operation of data from the host to the
2733:    device if the data on the device is out of date.  If the device
2734:    memory hasn't been allocated previously it will be allocated as part
2735:    of this function call.  VecHIPGetArrayRead() assumes that the
2736:    user will not modify the vector data.  This is analgogous to
2737:    intent(in) in Fortran.

2739:    The HIP device pointer has to be released by calling
2740:    VecHIPRestoreArrayRead().  If the data on the host side was
2741:    previously up to date it will remain so, i.e. data on both the device
2742:    and the host is up to date.  Accessing data on the host side does not
2743:    incur a device to host data transfer.

2745:    Input Parameter:
2746: .  v - the vector

2748:    Output Parameter:
2749: .  a - the HIP pointer.

2751:    Fortran note:
2752:    This function is not currently available from Fortran.

2754:    Level: intermediate

2756: .seealso: `VecHIPRestoreArrayRead()`, `VecHIPGetArray()`, `VecHIPGetArrayWrite()`, `VecGetArray()`, `VecGetArrayRead()`
2757: @*/
2758: PETSC_EXTERN PetscErrorCode VecHIPGetArrayRead(Vec v, const PetscScalar **a)
2759: {
2761: #if defined(PETSC_HAVE_HIP)
2762:   *a = 0;
2763:   VecHIPCopyToGPU(v);
2764:   *a = ((Vec_HIP *)v->spptr)->GPUarray;
2765: #endif
2766:   return 0;
2767: }

2769: /*@C
2770:    VecHIPRestoreArrayRead - Restore a HIP device pointer previously acquired with VecHIPGetArrayRead().

2772:    If the data on the host side was previously up to date it will remain
2773:    so, i.e. data on both the device and the host is up to date.
2774:    Accessing data on the host side e.g. with VecGetArray() does not
2775:    incur a device to host data transfer.

2777:    Input Parameters:
2778: +  v - the vector
2779: -  a - the HIP device pointer.  This pointer is invalid after
2780:        VecHIPRestoreArrayRead() returns.

2782:    Fortran note:
2783:    This function is not currently available from Fortran.

2785:    Level: intermediate

2787: .seealso: `VecHIPGetArrayRead()`, `VecHIPGetArrayWrite()`, `VecHIPGetArray()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`
2788: @*/
2789: PETSC_EXTERN PetscErrorCode VecHIPRestoreArrayRead(Vec v, const PetscScalar **a)
2790: {
2792:   *a = NULL;
2793:   return 0;
2794: }

2796: /*@C
2797:    VecHIPGetArrayWrite - Provides write access to the HIP buffer inside a vector.

2799:    The data pointed to by the device pointer is uninitialized.  The user
2800:    may not read from this data.  Furthermore, the entire array needs to
2801:    be filled by the user to obtain well-defined behaviour.  The device
2802:    memory will be allocated by this function if it hasn't been allocated
2803:    previously.  This is analogous to intent(out) in Fortran.

2805:    The device pointer needs to be released with
2806:    VecHIPRestoreArrayWrite().  When the pointer is released the
2807:    host data of the vector is marked as out of data.  Subsequent access
2808:    of the host data with e.g. VecGetArray() incurs a device to host data
2809:    transfer.

2811:    Input Parameter:
2812: .  v - the vector

2814:    Output Parameter:
2815: .  a - the HIP pointer

2817:    Fortran note:
2818:    This function is not currently available from Fortran.

2820:    Level: advanced

2822: .seealso: `VecHIPRestoreArrayWrite()`, `VecHIPGetArray()`, `VecHIPGetArrayRead()`, `VecHIPGetArrayWrite()`, `VecGetArray()`, `VecGetArrayRead()`
2823: @*/
2824: PETSC_EXTERN PetscErrorCode VecHIPGetArrayWrite(Vec v, PetscScalar **a)
2825: {
2827: #if defined(PETSC_HAVE_HIP)
2828:   *a = 0;
2829:   VecHIPAllocateCheck(v);
2830:   *a = ((Vec_HIP *)v->spptr)->GPUarray;
2831: #endif
2832:   return 0;
2833: }

2835: /*@C
2836:    VecHIPRestoreArrayWrite - Restore a HIP device pointer previously acquired with VecHIPGetArrayWrite().

2838:    Data on the host will be marked as out of date.  Subsequent access of
2839:    the data on the host side e.g. with VecGetArray() will incur a device
2840:    to host data transfer.

2842:    Input Parameters:
2843: +  v - the vector
2844: -  a - the HIP device pointer.  This pointer is invalid after
2845:        VecHIPRestoreArrayWrite() returns.

2847:    Fortran note:
2848:    This function is not currently available from Fortran.

2850:    Level: intermediate

2852: .seealso: `VecHIPGetArrayWrite()`, `VecHIPGetArray()`, `VecHIPGetArrayRead()`, `VecHIPGetArrayWrite()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`
2853: @*/
2854: PETSC_EXTERN PetscErrorCode VecHIPRestoreArrayWrite(Vec v, PetscScalar **a)
2855: {
2857: #if defined(PETSC_HAVE_HIP)
2858:   v->offloadmask = PETSC_OFFLOAD_GPU;
2859: #endif

2861:   PetscObjectStateIncrease((PetscObject)v);
2862:   return 0;
2863: }

2865: /*@C
2866:    VecHIPPlaceArray - Allows one to replace the GPU array in a vector with a
2867:    GPU array provided by the user. This is useful to avoid copying an
2868:    array into a vector.

2870:    Not Collective

2872:    Input Parameters:
2873: +  vec - the vector
2874: -  array - the GPU array

2876:    Notes:
2877:    You can return to the original GPU array with a call to VecHIPResetArray()
2878:    It is not possible to use VecHIPPlaceArray() and VecPlaceArray() at the
2879:    same time on the same vector.

2881:    `vec` does not take ownership of `array` in any way. The user must free `array` themselves
2882:    but be careful not to do so before the vector has either been destroyed, had its original
2883:    array restored with `VecHIPResetArray()` or permanently replaced with
2884:    `VecHIPReplaceArray()`.

2886:    Level: developer

2888: .seealso: `VecPlaceArray()`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`, `VecHIPResetArray()`, `VecHIPReplaceArray()`

2890: @*/
2891: PetscErrorCode VecHIPPlaceArray(Vec vin, const PetscScalar a[])
2892: {
2894: #if defined(PETSC_HAVE_HIP)
2895:   VecHIPCopyToGPU(vin);
2897:   ((Vec_Seq *)vin->data)->unplacedarray = (PetscScalar *)((Vec_HIP *)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2898:   ((Vec_HIP *)vin->spptr)->GPUarray     = (PetscScalar *)a;
2899:   vin->offloadmask                      = PETSC_OFFLOAD_GPU;
2900: #endif
2901:   PetscObjectStateIncrease((PetscObject)vin);
2902:   return 0;
2903: }

2905: /*@C
2906:    VecHIPReplaceArray - Allows one to replace the GPU array in a vector
2907:    with a GPU array provided by the user. This is useful to avoid copying
2908:    a GPU array into a vector.

2910:    Not Collective

2912:    Input Parameters:
2913: +  vec - the vector
2914: -  array - the GPU array

2916:    Notes:
2917:    This permanently replaces the GPU array and frees the memory associated
2918:    with the old GPU array.

2920:    The memory passed in CANNOT be freed by the user. It will be freed
2921:    when the vector is destroyed.

2923:    Not supported from Fortran

2925:    Level: developer

2927: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`, `VecHIPResetArray()`, `VecHIPPlaceArray()`, `VecReplaceArray()`

2929: @*/
2930: PetscErrorCode VecHIPReplaceArray(Vec vin, const PetscScalar a[])
2931: {
2933: #if defined(PETSC_HAVE_HIP)
2934:   hipFree(((Vec_HIP *)vin->spptr)->GPUarray);
2935:   ((Vec_HIP *)vin->spptr)->GPUarray = (PetscScalar *)a;
2936:   vin->offloadmask                  = PETSC_OFFLOAD_GPU;
2937: #endif
2938:   PetscObjectStateIncrease((PetscObject)vin);
2939:   return 0;
2940: }

2942: /*@C
2943:    VecHIPResetArray - Resets a vector to use its default memory. Call this
2944:    after the use of VecHIPPlaceArray().

2946:    Not Collective

2948:    Input Parameters:
2949: .  vec - the vector

2951:    Level: developer

2953: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecPlaceArray()`, `VecResetArray()`, `VecHIPPlaceArray()`, `VecHIPReplaceArray()`

2955: @*/
2956: PetscErrorCode VecHIPResetArray(Vec vin)
2957: {
2959: #if defined(PETSC_HAVE_HIP)
2960:   VecHIPCopyToGPU(vin);
2961:   ((Vec_HIP *)vin->spptr)->GPUarray     = (PetscScalar *)((Vec_Seq *)vin->data)->unplacedarray;
2962:   ((Vec_Seq *)vin->data)->unplacedarray = 0;
2963:   vin->offloadmask                      = PETSC_OFFLOAD_GPU;
2964: #endif
2965:   PetscObjectStateIncrease((PetscObject)vin);
2966:   return 0;
2967: }

2969: /*MC
2970:     VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
2971:     and makes them accessible via a Fortran90 pointer.

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

2976:     Collective on Vec

2978:     Input Parameters:
2979: +   x - a vector to mimic
2980: -   n - the number of vectors to obtain

2982:     Output Parameters:
2983: +   y - Fortran90 pointer to the array of vectors
2984: -   ierr - error code

2986:     Example of Usage:
2987: .vb
2988: #include <petsc/finclude/petscvec.h>
2989:     use petscvec

2991:     Vec x
2992:     Vec, pointer :: y(:)
2993:     ....
2994:     call VecDuplicateVecsF90(x,2,y,ierr)
2995:     call VecSet(y(2),alpha,ierr)
2996:     call VecSet(y(2),alpha,ierr)
2997:     ....
2998:     call VecDestroyVecsF90(2,y,ierr)
2999: .ve

3001:     Notes:
3002:     Not yet supported for all F90 compilers

3004:     Use VecDestroyVecsF90() to free the space.

3006:     Level: beginner

3008: .seealso: `VecDestroyVecsF90()`, `VecDuplicateVecs()`

3010: M*/

3012: /*MC
3013:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
3014:     VecGetArrayF90().

3016:     Synopsis:
3017:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

3019:     Logically Collective on Vec

3021:     Input Parameters:
3022: +   x - vector
3023: -   xx_v - the Fortran90 pointer to the array

3025:     Output Parameter:
3026: .   ierr - error code

3028:     Example of Usage:
3029: .vb
3030: #include <petsc/finclude/petscvec.h>
3031:     use petscvec

3033:     PetscScalar, pointer :: xx_v(:)
3034:     ....
3035:     call VecGetArrayF90(x,xx_v,ierr)
3036:     xx_v(3) = a
3037:     call VecRestoreArrayF90(x,xx_v,ierr)
3038: .ve

3040:     Level: beginner

3042: .seealso: `VecGetArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrayReadF90()`

3044: M*/

3046: /*MC
3047:     VecDestroyVecsF90 - Frees a block of vectors obtained with VecDuplicateVecsF90().

3049:     Synopsis:
3050:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

3052:     Collective on Vec

3054:     Input Parameters:
3055: +   n - the number of vectors previously obtained
3056: -   x - pointer to array of vector pointers

3058:     Output Parameter:
3059: .   ierr - error code

3061:     Notes:
3062:     Not yet supported for all F90 compilers

3064:     Level: beginner

3066: .seealso: `VecDestroyVecs()`, `VecDuplicateVecsF90()`

3068: M*/

3070: /*MC
3071:     VecGetArrayF90 - Accesses a vector array from Fortran90. For default PETSc
3072:     vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
3073:     this routine is implementation dependent. You MUST call VecRestoreArrayF90()
3074:     when you no longer need access to the array.

3076:     Synopsis:
3077:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

3079:     Logically Collective on Vec

3081:     Input Parameter:
3082: .   x - vector

3084:     Output Parameters:
3085: +   xx_v - the Fortran90 pointer to the array
3086: -   ierr - error code

3088:     Example of Usage:
3089: .vb
3090: #include <petsc/finclude/petscvec.h>
3091:     use petscvec

3093:     PetscScalar, pointer :: xx_v(:)
3094:     ....
3095:     call VecGetArrayF90(x,xx_v,ierr)
3096:     xx_v(3) = a
3097:     call VecRestoreArrayF90(x,xx_v,ierr)
3098: .ve

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

3102:     Level: beginner

3104: .seealso: `VecRestoreArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayReadF90()`

3106: M*/

3108: /*MC
3109:     VecGetArrayReadF90 - Accesses a read only array from Fortran90. For default PETSc
3110:     vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
3111:     this routine is implementation dependent. You MUST call VecRestoreArrayReadF90()
3112:     when you no longer need access to the array.

3114:     Synopsis:
3115:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

3117:     Logically Collective on Vec

3119:     Input Parameter:
3120: .   x - vector

3122:     Output Parameters:
3123: +   xx_v - the Fortran90 pointer to the array
3124: -   ierr - error code

3126:     Example of Usage:
3127: .vb
3128: #include <petsc/finclude/petscvec.h>
3129:     use petscvec

3131:     PetscScalar, pointer :: xx_v(:)
3132:     ....
3133:     call VecGetArrayReadF90(x,xx_v,ierr)
3134:     a = xx_v(3)
3135:     call VecRestoreArrayReadF90(x,xx_v,ierr)
3136: .ve

3138:     If you intend to write entries into the array you must use VecGetArrayF90().

3140:     Level: beginner

3142: .seealso: `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecGetArrayF90()`

3144: M*/

3146: /*MC
3147:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
3148:     VecGetArrayReadF90().

3150:     Synopsis:
3151:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

3153:     Logically Collective on Vec

3155:     Input Parameters:
3156: +   x - vector
3157: -   xx_v - the Fortran90 pointer to the array

3159:     Output Parameter:
3160: .   ierr - error code

3162:     Example of Usage:
3163: .vb
3164: #include <petsc/finclude/petscvec.h>
3165:     use petscvec

3167:     PetscScalar, pointer :: xx_v(:)
3168:     ....
3169:     call VecGetArrayReadF90(x,xx_v,ierr)
3170:     a = xx_v(3)
3171:     call VecRestoreArrayReadF90(x,xx_v,ierr)
3172: .ve

3174:     Level: beginner

3176: .seealso: `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecRestoreArrayF90()`

3178: M*/

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

3185:    Logically Collective

3187:    Input Parameters:
3188: +  x - the vector
3189: .  m - first dimension of two dimensional array
3190: .  n - second dimension of two dimensional array
3191: .  mstart - first index you will use in first coordinate direction (often 0)
3192: -  nstart - first index in the second coordinate direction (often 0)

3194:    Output Parameter:
3195: .  a - location to put pointer to the array

3197:    Level: developer

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

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

3207: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3208:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3209:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3210: @*/
3211: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3212: {
3213:   PetscInt     i, N;
3214:   PetscScalar *aa;

3219:   VecGetLocalSize(x, &N);
3221:   VecGetArray(x, &aa);

3223:   PetscMalloc1(m, a);
3224:   for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
3225:   *a -= mstart;
3226:   return 0;
3227: }

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

3234:    Logically Collective

3236:    Input Parameters:
3237: +  x - the vector
3238: .  m - first dimension of two dimensional array
3239: .  n - second dimension of two dimensional array
3240: .  mstart - first index you will use in first coordinate direction (often 0)
3241: -  nstart - first index in the second coordinate direction (often 0)

3243:    Output Parameter:
3244: .  a - location to put pointer to the array

3246:    Level: developer

3248:   Notes:
3249:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
3250:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3251:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3252:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

3256: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3257:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3258:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3259: @*/
3260: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3261: {
3262:   PetscInt     i, N;
3263:   PetscScalar *aa;

3268:   VecGetLocalSize(x, &N);
3270:   VecGetArrayWrite(x, &aa);

3272:   PetscMalloc1(m, a);
3273:   for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
3274:   *a -= mstart;
3275:   return 0;
3276: }

3278: /*@C
3279:    VecRestoreArray2d - Restores a vector after VecGetArray2d() has been called.

3281:    Logically Collective

3283:    Input Parameters:
3284: +  x - the vector
3285: .  m - first dimension of two dimensional array
3286: .  n - second dimension of the two dimensional array
3287: .  mstart - first index you will use in first coordinate direction (often 0)
3288: .  nstart - first index in the second coordinate direction (often 0)
3289: -  a - location of pointer to array obtained from VecGetArray2d()

3291:    Level: developer

3293:    Notes:
3294:    For regular PETSc vectors this routine does not involve any copies. For
3295:    any special vectors that do not store local vector data in a contiguous
3296:    array, this routine will copy the data back into the underlying
3297:    vector data structure from the array obtained with VecGetArray().

3299:    This routine actually zeros out the a pointer.

3301: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3302:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3303:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3304: @*/
3305: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3306: {
3307:   void *dummy;

3312:   dummy = (void *)(*a + mstart);
3313:   PetscFree(dummy);
3314:   VecRestoreArray(x, NULL);
3315:   return 0;
3316: }

3318: /*@C
3319:    VecRestoreArray2dWrite - Restores a vector after VecGetArray2dWrite() has been called.

3321:    Logically Collective

3323:    Input Parameters:
3324: +  x - the vector
3325: .  m - first dimension of two dimensional array
3326: .  n - second dimension of the two dimensional array
3327: .  mstart - first index you will use in first coordinate direction (often 0)
3328: .  nstart - first index in the second coordinate direction (often 0)
3329: -  a - location of pointer to array obtained from VecGetArray2d()

3331:    Level: developer

3333:    Notes:
3334:    For regular PETSc vectors this routine does not involve any copies. For
3335:    any special vectors that do not store local vector data in a contiguous
3336:    array, this routine will copy the data back into the underlying
3337:    vector data structure from the array obtained with VecGetArray().

3339:    This routine actually zeros out the a pointer.

3341: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3342:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3343:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3344: @*/
3345: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3346: {
3347:   void *dummy;

3352:   dummy = (void *)(*a + mstart);
3353:   PetscFree(dummy);
3354:   VecRestoreArrayWrite(x, NULL);
3355:   return 0;
3356: }

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

3363:    Logically Collective

3365:    Input Parameters:
3366: +  x - the vector
3367: .  m - first dimension of two dimensional array
3368: -  mstart - first index you will use in first coordinate direction (often 0)

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

3373:    Level: developer

3375:   Notes:
3376:    For a vector obtained from DMCreateLocalVector() mstart are likely
3377:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3378:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3382: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3383:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3384:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3385: @*/
3386: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3387: {
3388:   PetscInt N;

3393:   VecGetLocalSize(x, &N);
3395:   VecGetArray(x, a);
3396:   *a -= mstart;
3397:   return 0;
3398: }

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

3405:    Logically Collective

3407:    Input Parameters:
3408: +  x - the vector
3409: .  m - first dimension of two dimensional array
3410: -  mstart - first index you will use in first coordinate direction (often 0)

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

3415:    Level: developer

3417:   Notes:
3418:    For a vector obtained from DMCreateLocalVector() mstart are likely
3419:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3420:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3424: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3425:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3426:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3427: @*/
3428: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3429: {
3430:   PetscInt N;

3435:   VecGetLocalSize(x, &N);
3437:   VecGetArrayWrite(x, a);
3438:   *a -= mstart;
3439:   return 0;
3440: }

3442: /*@C
3443:    VecRestoreArray1d - Restores a vector after VecGetArray1d() has been called.

3445:    Logically Collective

3447:    Input Parameters:
3448: +  x - the vector
3449: .  m - first dimension of two dimensional array
3450: .  mstart - first index you will use in first coordinate direction (often 0)
3451: -  a - location of pointer to array obtained from VecGetArray21()

3453:    Level: developer

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

3461:    This routine actually zeros out the a pointer.

3463: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3464:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3465:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3466: @*/
3467: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3468: {
3471:   VecRestoreArray(x, NULL);
3472:   return 0;
3473: }

3475: /*@C
3476:    VecRestoreArray1dWrite - Restores a vector after VecGetArray1dWrite() has been called.

3478:    Logically Collective

3480:    Input Parameters:
3481: +  x - the vector
3482: .  m - first dimension of two dimensional array
3483: .  mstart - first index you will use in first coordinate direction (often 0)
3484: -  a - location of pointer to array obtained from VecGetArray21()

3486:    Level: developer

3488:    Notes:
3489:    For regular PETSc vectors this routine does not involve any copies. For
3490:    any special vectors that do not store local vector data in a contiguous
3491:    array, this routine will copy the data back into the underlying
3492:    vector data structure from the array obtained with VecGetArray1d().

3494:    This routine actually zeros out the a pointer.

3496: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3497:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3498:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3499: @*/
3500: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3501: {
3504:   VecRestoreArrayWrite(x, NULL);
3505:   return 0;
3506: }

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

3513:    Logically Collective

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

3524:    Output Parameter:
3525: .  a - location to put pointer to the array

3527:    Level: developer

3529:   Notes:
3530:    For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3531:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3532:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3533:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().

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

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

3549:   VecGetLocalSize(x, &N);
3551:   VecGetArray(x, &aa);

3553:   PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a);
3554:   b = (PetscScalar **)((*a) + m);
3555:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3556:   for (i = 0; i < m; i++)
3557:     for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3558:   *a -= mstart;
3559:   return 0;
3560: }

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

3567:    Logically Collective

3569:    Input Parameters:
3570: +  x - the vector
3571: .  m - first dimension of three dimensional array
3572: .  n - second dimension of three dimensional array
3573: .  p - third dimension of three dimensional array
3574: .  mstart - first index you will use in first coordinate direction (often 0)
3575: .  nstart - first index in the second coordinate direction (often 0)
3576: -  pstart - first index in the third coordinate direction (often 0)

3578:    Output Parameter:
3579: .  a - location to put pointer to the array

3581:    Level: developer

3583:   Notes:
3584:    For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3585:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3586:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3587:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().

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

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

3603:   VecGetLocalSize(x, &N);
3605:   VecGetArrayWrite(x, &aa);

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

3613:   *a -= mstart;
3614:   return 0;
3615: }

3617: /*@C
3618:    VecRestoreArray3d - Restores a vector after VecGetArray3d() has been called.

3620:    Logically Collective

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

3632:    Level: developer

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

3640:    This routine actually zeros out the a pointer.

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

3653:   dummy = (void *)(*a + mstart);
3654:   PetscFree(dummy);
3655:   VecRestoreArray(x, NULL);
3656:   return 0;
3657: }

3659: /*@C
3660:    VecRestoreArray3dWrite - Restores a vector after VecGetArray3dWrite() has been called.

3662:    Logically Collective

3664:    Input Parameters:
3665: +  x - the vector
3666: .  m - first dimension of three dimensional array
3667: .  n - second dimension of the three dimensional array
3668: .  p - third dimension of the three dimensional array
3669: .  mstart - first index you will use in first coordinate direction (often 0)
3670: .  nstart - first index in the second coordinate direction (often 0)
3671: .  pstart - first index in the third coordinate direction (often 0)
3672: -  a - location of pointer to array obtained from VecGetArray3d()

3674:    Level: developer

3676:    Notes:
3677:    For regular PETSc vectors this routine does not involve any copies. For
3678:    any special vectors that do not store local vector data in a contiguous
3679:    array, this routine will copy the data back into the underlying
3680:    vector data structure from the array obtained with VecGetArray().

3682:    This routine actually zeros out the a pointer.

3684: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3685:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3686:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3687: @*/
3688: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3689: {
3690:   void *dummy;

3695:   dummy = (void *)(*a + mstart);
3696:   PetscFree(dummy);
3697:   VecRestoreArrayWrite(x, NULL);
3698:   return 0;
3699: }

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

3706:    Logically Collective

3708:    Input Parameters:
3709: +  x - the vector
3710: .  m - first dimension of four dimensional array
3711: .  n - second dimension of four dimensional array
3712: .  p - third dimension of four dimensional array
3713: .  q - fourth dimension of four dimensional array
3714: .  mstart - first index you will use in first coordinate direction (often 0)
3715: .  nstart - first index in the second coordinate direction (often 0)
3716: .  pstart - first index in the third coordinate direction (often 0)
3717: -  qstart - first index in the fourth coordinate direction (often 0)

3719:    Output Parameter:
3720: .  a - location to put pointer to the array

3722:    Level: beginner

3724:   Notes:
3725:    For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3726:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3727:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3728:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().

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

3732: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3733:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3734:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3735: @*/
3736: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3737: {
3738:   PetscInt     i, N, j, k;
3739:   PetscScalar *aa, ***b, **c;

3744:   VecGetLocalSize(x, &N);
3746:   VecGetArray(x, &aa);

3748:   PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a);
3749:   b = (PetscScalar ***)((*a) + m);
3750:   c = (PetscScalar **)(b + m * n);
3751:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3752:   for (i = 0; i < m; i++)
3753:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3754:   for (i = 0; i < m; i++)
3755:     for (j = 0; j < n; j++)
3756:       for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3757:   *a -= mstart;
3758:   return 0;
3759: }

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

3766:    Logically Collective

3768:    Input Parameters:
3769: +  x - the vector
3770: .  m - first dimension of four dimensional array
3771: .  n - second dimension of four dimensional array
3772: .  p - third dimension of four dimensional array
3773: .  q - fourth dimension of four dimensional array
3774: .  mstart - first index you will use in first coordinate direction (often 0)
3775: .  nstart - first index in the second coordinate direction (often 0)
3776: .  pstart - first index in the third coordinate direction (often 0)
3777: -  qstart - first index in the fourth coordinate direction (often 0)

3779:    Output Parameter:
3780: .  a - location to put pointer to the array

3782:    Level: beginner

3784:   Notes:
3785:    For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3786:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3787:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3788:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().

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

3792: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3793:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3794:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3795: @*/
3796: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3797: {
3798:   PetscInt     i, N, j, k;
3799:   PetscScalar *aa, ***b, **c;

3804:   VecGetLocalSize(x, &N);
3806:   VecGetArrayWrite(x, &aa);

3808:   PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a);
3809:   b = (PetscScalar ***)((*a) + m);
3810:   c = (PetscScalar **)(b + m * n);
3811:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3812:   for (i = 0; i < m; i++)
3813:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3814:   for (i = 0; i < m; i++)
3815:     for (j = 0; j < n; j++)
3816:       for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3817:   *a -= mstart;
3818:   return 0;
3819: }

3821: /*@C
3822:    VecRestoreArray4d - Restores a vector after VecGetArray3d() has been called.

3824:    Logically Collective

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

3838:    Level: beginner

3840:    Notes:
3841:    For regular PETSc vectors this routine does not involve any copies. For
3842:    any special vectors that do not store local vector data in a contiguous
3843:    array, this routine will copy the data back into the underlying
3844:    vector data structure from the array obtained with VecGetArray().

3846:    This routine actually zeros out the a pointer.

3848: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3849:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3850:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3851: @*/
3852: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3853: {
3854:   void *dummy;

3859:   dummy = (void *)(*a + mstart);
3860:   PetscFree(dummy);
3861:   VecRestoreArray(x, NULL);
3862:   return 0;
3863: }

3865: /*@C
3866:    VecRestoreArray4dWrite - Restores a vector after VecGetArray3dWrite() has been called.

3868:    Logically Collective

3870:    Input Parameters:
3871: +  x - the vector
3872: .  m - first dimension of four dimensional array
3873: .  n - second dimension of the four dimensional array
3874: .  p - third dimension of the four dimensional array
3875: .  q - fourth dimension of the four dimensional array
3876: .  mstart - first index you will use in first coordinate direction (often 0)
3877: .  nstart - first index in the second coordinate direction (often 0)
3878: .  pstart - first index in the third coordinate direction (often 0)
3879: .  qstart - first index in the fourth coordinate direction (often 0)
3880: -  a - location of pointer to array obtained from VecGetArray4d()

3882:    Level: beginner

3884:    Notes:
3885:    For regular PETSc vectors this routine does not involve any copies. For
3886:    any special vectors that do not store local vector data in a contiguous
3887:    array, this routine will copy the data back into the underlying
3888:    vector data structure from the array obtained with VecGetArray().

3890:    This routine actually zeros out the a pointer.

3892: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3893:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3894:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3895: @*/
3896: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3897: {
3898:   void *dummy;

3903:   dummy = (void *)(*a + mstart);
3904:   PetscFree(dummy);
3905:   VecRestoreArrayWrite(x, NULL);
3906:   return 0;
3907: }

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

3914:    Logically Collective

3916:    Input Parameters:
3917: +  x - the vector
3918: .  m - first dimension of two dimensional array
3919: .  n - second dimension of two dimensional array
3920: .  mstart - first index you will use in first coordinate direction (often 0)
3921: -  nstart - first index in the second coordinate direction (often 0)

3923:    Output Parameter:
3924: .  a - location to put pointer to the array

3926:    Level: developer

3928:   Notes:
3929:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
3930:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3931:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3932:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

3936: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3937:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3938:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3939: @*/
3940: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3941: {
3942:   PetscInt           i, N;
3943:   const PetscScalar *aa;

3948:   VecGetLocalSize(x, &N);
3950:   VecGetArrayRead(x, &aa);

3952:   PetscMalloc1(m, a);
3953:   for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3954:   *a -= mstart;
3955:   return 0;
3956: }

3958: /*@C
3959:    VecRestoreArray2dRead - Restores a vector after VecGetArray2dRead() has been called.

3961:    Logically Collective

3963:    Input Parameters:
3964: +  x - the vector
3965: .  m - first dimension of two dimensional array
3966: .  n - second dimension of the two dimensional array
3967: .  mstart - first index you will use in first coordinate direction (often 0)
3968: .  nstart - first index in the second coordinate direction (often 0)
3969: -  a - location of pointer to array obtained from VecGetArray2d()

3971:    Level: developer

3973:    Notes:
3974:    For regular PETSc vectors this routine does not involve any copies. For
3975:    any special vectors that do not store local vector data in a contiguous
3976:    array, this routine will copy the data back into the underlying
3977:    vector data structure from the array obtained with VecGetArray().

3979:    This routine actually zeros out the a pointer.

3981: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3982:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3983:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3984: @*/
3985: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3986: {
3987:   void *dummy;

3992:   dummy = (void *)(*a + mstart);
3993:   PetscFree(dummy);
3994:   VecRestoreArrayRead(x, NULL);
3995:   return 0;
3996: }

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

4003:    Logically Collective

4005:    Input Parameters:
4006: +  x - the vector
4007: .  m - first dimension of two dimensional array
4008: -  mstart - first index you will use in first coordinate direction (often 0)

4010:    Output Parameter:
4011: .  a - location to put pointer to the array

4013:    Level: developer

4015:   Notes:
4016:    For a vector obtained from DMCreateLocalVector() mstart are likely
4017:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
4018:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

4022: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
4023:           `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
4024:           `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
4025: @*/
4026: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
4027: {
4028:   PetscInt N;

4033:   VecGetLocalSize(x, &N);
4035:   VecGetArrayRead(x, (const PetscScalar **)a);
4036:   *a -= mstart;
4037:   return 0;
4038: }

4040: /*@C
4041:    VecRestoreArray1dRead - Restores a vector after VecGetArray1dRead() has been called.

4043:    Logically Collective

4045:    Input Parameters:
4046: +  x - the vector
4047: .  m - first dimension of two dimensional array
4048: .  mstart - first index you will use in first coordinate direction (often 0)
4049: -  a - location of pointer to array obtained from VecGetArray21()

4051:    Level: developer

4053:    Notes:
4054:    For regular PETSc vectors this routine does not involve any copies. For
4055:    any special vectors that do not store local vector data in a contiguous
4056:    array, this routine will copy the data back into the underlying
4057:    vector data structure from the array obtained with VecGetArray1dRead().

4059:    This routine actually zeros out the a pointer.

4061: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
4062:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
4063:           `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
4064: @*/
4065: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
4066: {
4069:   VecRestoreArrayRead(x, NULL);
4070:   return 0;
4071: }

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

4078:    Logically Collective

4080:    Input Parameters:
4081: +  x - the vector
4082: .  m - first dimension of three dimensional array
4083: .  n - second dimension of three dimensional array
4084: .  p - third dimension of three dimensional array
4085: .  mstart - first index you will use in first coordinate direction (often 0)
4086: .  nstart - first index in the second coordinate direction (often 0)
4087: -  pstart - first index in the third coordinate direction (often 0)

4089:    Output Parameter:
4090: .  a - location to put pointer to the array

4092:    Level: developer

4094:   Notes:
4095:    For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
4096:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
4097:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
4098:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3dRead().

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

4102: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
4103:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
4104:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
4105: @*/
4106: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
4107: {
4108:   PetscInt           i, N, j;
4109:   const PetscScalar *aa;
4110:   PetscScalar      **b;

4115:   VecGetLocalSize(x, &N);
4117:   VecGetArrayRead(x, &aa);

4119:   PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a);
4120:   b = (PetscScalar **)((*a) + m);
4121:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
4122:   for (i = 0; i < m; i++)
4123:     for (j = 0; j < n; j++) b[i * n + j] = (PetscScalar *)aa + i * n * p + j * p - pstart;
4124:   *a -= mstart;
4125:   return 0;
4126: }

4128: /*@C
4129:    VecRestoreArray3dRead - Restores a vector after VecGetArray3dRead() has been called.

4131:    Logically Collective

4133:    Input Parameters:
4134: +  x - the vector
4135: .  m - first dimension of three dimensional array
4136: .  n - second dimension of the three dimensional array
4137: .  p - third dimension of the three dimensional array
4138: .  mstart - first index you will use in first coordinate direction (often 0)
4139: .  nstart - first index in the second coordinate direction (often 0)
4140: .  pstart - first index in the third coordinate direction (often 0)
4141: -  a - location of pointer to array obtained from VecGetArray3dRead()

4143:    Level: developer

4145:    Notes:
4146:    For regular PETSc vectors this routine does not involve any copies. For
4147:    any special vectors that do not store local vector data in a contiguous
4148:    array, this routine will copy the data back into the underlying
4149:    vector data structure from the array obtained with VecGetArray().

4151:    This routine actually zeros out the a pointer.

4153: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
4154:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
4155:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
4156: @*/
4157: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
4158: {
4159:   void *dummy;

4164:   dummy = (void *)(*a + mstart);
4165:   PetscFree(dummy);
4166:   VecRestoreArrayRead(x, NULL);
4167:   return 0;
4168: }

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

4175:    Logically Collective

4177:    Input Parameters:
4178: +  x - the vector
4179: .  m - first dimension of four dimensional array
4180: .  n - second dimension of four dimensional array
4181: .  p - third dimension of four dimensional array
4182: .  q - fourth dimension of four dimensional array
4183: .  mstart - first index you will use in first coordinate direction (often 0)
4184: .  nstart - first index in the second coordinate direction (often 0)
4185: .  pstart - first index in the third coordinate direction (often 0)
4186: -  qstart - first index in the fourth coordinate direction (often 0)

4188:    Output Parameter:
4189: .  a - location to put pointer to the array

4191:    Level: beginner

4193:   Notes:
4194:    For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
4195:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
4196:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
4197:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().

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

4201: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
4202:           `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
4203:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
4204: @*/
4205: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
4206: {
4207:   PetscInt           i, N, j, k;
4208:   const PetscScalar *aa;
4209:   PetscScalar     ***b, **c;

4214:   VecGetLocalSize(x, &N);
4216:   VecGetArrayRead(x, &aa);

4218:   PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a);
4219:   b = (PetscScalar ***)((*a) + m);
4220:   c = (PetscScalar **)(b + m * n);
4221:   for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
4222:   for (i = 0; i < m; i++)
4223:     for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
4224:   for (i = 0; i < m; i++)
4225:     for (j = 0; j < n; j++)
4226:       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;
4227:   *a -= mstart;
4228:   return 0;
4229: }

4231: /*@C
4232:    VecRestoreArray4dRead - Restores a vector after VecGetArray3d() has been called.

4234:    Logically Collective

4236:    Input Parameters:
4237: +  x - the vector
4238: .  m - first dimension of four dimensional array
4239: .  n - second dimension of the four dimensional array
4240: .  p - third dimension of the four dimensional array
4241: .  q - fourth dimension of the four dimensional array
4242: .  mstart - first index you will use in first coordinate direction (often 0)
4243: .  nstart - first index in the second coordinate direction (often 0)
4244: .  pstart - first index in the third coordinate direction (often 0)
4245: .  qstart - first index in the fourth coordinate direction (often 0)
4246: -  a - location of pointer to array obtained from VecGetArray4dRead()

4248:    Level: beginner

4250:    Notes:
4251:    For regular PETSc vectors this routine does not involve any copies. For
4252:    any special vectors that do not store local vector data in a contiguous
4253:    array, this routine will copy the data back into the underlying
4254:    vector data structure from the array obtained with VecGetArray().

4256:    This routine actually zeros out the a pointer.

4258: .seealso: `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
4259:           `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
4260:           `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
4261: @*/
4262: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
4263: {
4264:   void *dummy;

4269:   dummy = (void *)(*a + mstart);
4270:   PetscFree(dummy);
4271:   VecRestoreArrayRead(x, NULL);
4272:   return 0;
4273: }

4275: #if defined(PETSC_USE_DEBUG)

4277: /*@
4278:    VecLockGet  - Gets the current lock status of a vector

4280:    Logically Collective on Vec

4282:    Input Parameter:
4283: .  x - the vector

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

4289:    Level: beginner

4291: .seealso: `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
4292: @*/
4293: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
4294: {
4296:   *state = x->lock;
4297:   return 0;
4298: }

4300: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
4301: {
4306:   {
4307:     const int index = x->lockstack.currentsize - 1;

4310:     *file = x->lockstack.file[index];
4311:     *func = x->lockstack.function[index];
4312:     *line = x->lockstack.line[index];
4313:   }
4314:   return 0;
4315: }

4317: /*@
4318:    VecLockReadPush  - Pushes a read-only lock on a vector to prevent it from writing

4320:    Logically Collective on Vec

4322:    Input Parameter:
4323: .  x - the vector

4325:    Notes:
4326:     If this is set then calls to VecGetArray() or VecSetValues() or any other routines that change the vectors values will fail.

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

4331:    Level: beginner

4333: .seealso: `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
4334: @*/
4335: PetscErrorCode VecLockReadPush(Vec x)
4336: {
4337:   const char *file, *func;
4338:   int         index, line;

4342:   if ((index = petscstack.currentsize - 2) == -1) {
4343:     // vec was locked "outside" of petsc, either in user-land or main. the error message will
4344:     // now show this function as the culprit, but it will include the stacktrace
4345:     file = "unknown user-file";
4346:     func = "unknown_user_function";
4347:     line = 0;
4348:   } else {
4350:     file = petscstack.file[index];
4351:     func = petscstack.function[index];
4352:     line = petscstack.line[index];
4353:   }
4354:   PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
4355:   return 0;
4356: }

4358: /*@
4359:    VecLockReadPop  - Pops a read-only lock from a vector

4361:    Logically Collective on Vec

4363:    Input Parameter:
4364: .  x - the vector

4366:    Level: beginner

4368: .seealso: `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
4369: @*/
4370: PetscErrorCode VecLockReadPop(Vec x)
4371: {
4374:   {
4375:     const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];

4377:     PetscStackPop_Private(x->lockstack, previous);
4378:   }
4379:   return 0;
4380: }

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

4385:    Logically Collective on Vec

4387:    Input Parameters:
4388: +  x   - the vector
4389: -  flg - PETSC_TRUE to lock the vector for exclusive read/write access; PETSC_FALSE to unlock it.

4391:    Notes:
4392:     The function is usefull in split-phase computations, which usually have a begin phase and an end phase.
4393:     One can call VecLockWriteSet(x,PETSC_TRUE) in the begin phase to lock a vector for exclusive
4394:     access, and call VecLockWriteSet(x,PETSC_FALSE) in the end phase to unlock the vector from exclusive
4395:     access. In this way, one is ensured no other operations can access the vector in between. The code may like

4397:        VecGetArray(x,&xdata); // begin phase
4398:        VecLockWriteSet(v,PETSC_TRUE);

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

4402:        VecRestoreArray(x,&vdata); // end phase
4403:        VecLockWriteSet(v,PETSC_FALSE);

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

4408:    Level: beginner

4410: .seealso: `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
4411: @*/
4412: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
4413: {
4415:   if (flg) {
4418:     x->lock = -1;
4419:   } else {
4421:     x->lock = 0;
4422:   }
4423:   return 0;
4424: }

4426: /*@
4427:    VecLockPush  - Pushes a read-only lock on a vector to prevent it from writing

4429:    Level: deprecated

4431: .seealso: `VecLockReadPush()`
4432: @*/
4433: PetscErrorCode VecLockPush(Vec x)
4434: {
4435:   VecLockReadPush(x);
4436:   return 0;
4437: }

4439: /*@
4440:    VecLockPop  - Pops a read-only lock from a vector

4442:    Level: deprecated

4444: .seealso: `VecLockReadPop()`
4445: @*/
4446: PetscErrorCode VecLockPop(Vec x)
4447: {
4448:   VecLockReadPop(x);
4449:   return 0;
4450: }

4452: #endif