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: PETSC_EXTERN PetscErrorCode VecValidValues(Vec vec,PetscInt argnum,PetscBool begin)
 19: {
 20: #if defined(PETSC_USE_DEBUG)
 21:   PetscInt          n,i;
 22:   const PetscScalar *x;

 24: #if defined(PETSC_HAVE_DEVICE)
 25:   if ((vec->petscnative || vec->ops->getarray) && (vec->offloadmask & PETSC_OFFLOAD_CPU)) {
 26: #else
 27:   if (vec->petscnative || vec->ops->getarray) {
 28: #endif
 29:     VecGetLocalSize(vec,&n);
 30:     VecGetArrayRead(vec,&x);
 31:     for (i=0; i<n; i++) {
 32:       if (begin) {
 34:       } else {
 36:       }
 37:     }
 38:     VecRestoreArrayRead(vec,&x);
 39:   }
 40: #else
 41: #endif
 42:   return 0;
 43: }

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

 48:    Logically Collective on Vec

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

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

 56:    Level: advanced

 58:    Notes:
 59:     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: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs()
 63: @*/
 64: PetscErrorCode  VecMaxPointwiseDivide(Vec x,Vec y,PetscReal *max)
 65: {
 72:   VecCheckSameSize(x,1,y,2);
 73:   (*x->ops->maxpointwisedivide)(x,y,max);
 74:   return 0;
 75: }

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

 80:    Collective on Vec

 82:    Input Parameters:
 83: .  x, y - the vectors

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

 88:    Performance Issues:
 89: $    per-processor memory bandwidth
 90: $    interprocessor latency
 91: $    work load imbalance that causes certain processes to arrive much earlier than others

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

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

104:    Level: intermediate

106: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDotRealPart()
107: @*/
108: PetscErrorCode  VecDot(Vec x,Vec y,PetscScalar *val)
109: {
116:   VecCheckSameSize(x,1,y,2);

118:   PetscLogEventBegin(VEC_Dot,x,y,0,0);
119:   (*x->ops->dot)(x,y,val);
120:   PetscLogEventEnd(VEC_Dot,x,y,0,0);
121:   return 0;
122: }

124: /*@
125:    VecDotRealPart - Computes the real part of the vector dot product.

127:    Collective on Vec

129:    Input Parameters:
130: .  x, y - the vectors

132:    Output Parameter:
133: .  val - the real part of the dot product;

135:    Performance Issues:
136: $    per-processor memory bandwidth
137: $    interprocessor latency
138: $    work load imbalance that causes certain processes to arrive much earlier than others

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

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

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

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

150:    Level: intermediate

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

158:   VecDot(x,y,&fdot);
159:   *val = PetscRealPart(fdot);
160:   return 0;
161: }

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

166:    Collective on Vec

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

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

175:    Values of NormType:
176: +     NORM_1 - sum_i |x_i|
177: .     NORM_2 - sqrt(sum_i |x_i|^2)
178: .     NORM_INFINITY - max_i |x_i|
179: -     NORM_1_AND_2 - computes efficiently both  NORM_1 and NORM_2 and stores them each in an output array

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

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

192:    Level: intermediate

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

200: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNormAvailable(),
201:           VecNormBegin(), VecNormEnd(), NormType()

203: @*/
204: PetscErrorCode  VecNorm(Vec x,NormType type,PetscReal *val)
205: {
206:   PetscBool      flg;


212:   /*
213:    * Cached data?
214:    */
215:   if (type!=NORM_1_AND_2) {
216:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,flg);
217:     if (flg) return 0;
218:   }
219:   PetscLogEventBegin(VEC_Norm,x,0,0,0);
220:   (*x->ops->norm)(x,type,val);
221:   PetscLogEventEnd(VEC_Norm,x,0,0,0);
222:   if (type!=NORM_1_AND_2) {
223:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[type],*val);
224:   }
225:   return 0;
226: }

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

231:    Not Collective

233:    Input Parameters:
234: +  x - the vector
235: -  type - one of NORM_1, NORM_2, NORM_INFINITY.  Also available
236:           NORM_1_AND_2, which computes both norms and stores them
237:           in a two element array.

239:    Output Parameters:
240: +  available - PETSC_TRUE if the val returned is valid
241: -  val - the norm

243:    Notes:
244: $     NORM_1 denotes sum_i |x_i|
245: $     NORM_2 denotes sqrt(sum_i (x_i)^2)
246: $     NORM_INFINITY denotes max_i |x_i|

248:    Level: intermediate

250:    Performance Issues:
251: $    per-processor memory bandwidth
252: $    interprocessor latency
253: $    work load imbalance that causes certain processes to arrive much earlier than others

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

260: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNorm()
261:           VecNormBegin(), VecNormEnd()

263: @*/
264: PetscErrorCode  VecNormAvailable(Vec x,NormType type,PetscBool  *available,PetscReal *val)
265: {

270:   *available = PETSC_FALSE;
271:   if (type!=NORM_1_AND_2) {
272:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,*available);
273:   }
274:   return 0;
275: }

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

280:    Collective on Vec

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

285:    Output Parameter:
286: .  val - the vector norm before normalization

288:    Level: intermediate

290: @*/
291: PetscErrorCode  VecNormalize(Vec x,PetscReal *val)
292: {
293:   PetscReal      norm;

297:   PetscLogEventBegin(VEC_Normalize,x,0,0,0);
298:   VecNorm(x,NORM_2,&norm);
299:   if (norm == 0.0) {
300:     PetscInfo(x,"Vector of zero norm can not be normalized; Returning only the zero norm\n");
301:   } else if (norm != 1.0) {
302:     PetscScalar tmp = 1.0/norm;
303:     VecScale(x,tmp);
304:   }
305:   if (val) *val = norm;
306:   PetscLogEventEnd(VEC_Normalize,x,0,0,0);
307:   return 0;
308: }

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

313:    Collective on Vec

315:    Input Parameter:
316: .  x - the vector

318:    Output Parameters:
319: +  p - the location of val (pass NULL if you don't want this)
320: -  val - the maximum component

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

325:    Returns the smallest index with the maximum value
326:    Level: intermediate

328: .seealso: VecNorm(), VecMin()
329: @*/
330: PetscErrorCode  VecMax(Vec x,PetscInt *p,PetscReal *val)
331: {
335:   PetscLogEventBegin(VEC_Max,x,0,0,0);
336:   (*x->ops->max)(x,p,val);
337:   PetscLogEventEnd(VEC_Max,x,0,0,0);
338:   return 0;
339: }

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

344:    Collective on Vec

346:    Input Parameter:
347: .  x - the vector

349:    Output Parameters:
350: +  p - the location of val (pass NULL if you don't want this location)
351: -  val - the minimum component

353:    Level: intermediate

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

358:    This returns the smallest index with the minumum value

360: .seealso: VecMax()
361: @*/
362: PetscErrorCode  VecMin(Vec x,PetscInt *p,PetscReal *val)
363: {
367:   PetscLogEventBegin(VEC_Min,x,0,0,0);
368:   (*x->ops->min)(x,p,val);
369:   PetscLogEventEnd(VEC_Min,x,0,0,0);
370:   return 0;
371: }

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

377:    Collective on Vec

379:    Input Parameters:
380: .  x, y - the vectors

382:    Output Parameter:
383: .  val - the dot product

385:    Notes for Users of Complex Numbers:
386:    For complex vectors, VecTDot() computes the indefinite form
387: $     val = (x,y) = y^T x,
388:    where y^T denotes the transpose of y.

390:    Use VecDot() for the inner product
391: $     val = (x,y) = y^H x,
392:    where y^H denotes the conjugate transpose of y.

394:    Level: intermediate

396: .seealso: VecDot(), VecMTDot()
397: @*/
398: PetscErrorCode  VecTDot(Vec x,Vec y,PetscScalar *val)
399: {
406:   VecCheckSameSize(x,1,y,2);

408:   PetscLogEventBegin(VEC_TDot,x,y,0,0);
409:   (*x->ops->tdot)(x,y,val);
410:   PetscLogEventEnd(VEC_TDot,x,y,0,0);
411:   return 0;
412: }

414: /*@
415:    VecScale - Scales a vector.

417:    Not collective on Vec

419:    Input Parameters:
420: +  x - the vector
421: -  alpha - the scalar

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

427:    Level: intermediate

429: @*/
430: PetscErrorCode  VecScale(Vec x, PetscScalar alpha)
431: {
432:   PetscReal      norms[4] = {0.0,0.0,0.0, 0.0};
433:   PetscBool      flgs[4];
434:   PetscInt       i;

439:   PetscLogEventBegin(VEC_Scale,x,0,0,0);
440:   if (alpha != (PetscScalar)1.0) {
441:     VecSetErrorIfLocked(x,1);
442:     /* get current stashed norms */
443:     for (i=0; i<4; i++) {
444:       PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
445:     }
446:     (*x->ops->scale)(x,alpha);
447:     PetscObjectStateIncrease((PetscObject)x);
448:     /* put the scaled stashed norms back into the Vec */
449:     for (i=0; i<4; i++) {
450:       if (flgs[i]) {
451:         PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],PetscAbsScalar(alpha)*norms[i]);
452:       }
453:     }
454:   }
455:   PetscLogEventEnd(VEC_Scale,x,0,0,0);
456:   return 0;
457: }

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

462:    Logically Collective on Vec

464:    Input Parameters:
465: +  x  - the vector
466: -  alpha - the scalar

468:    Output Parameter:
469: .  x  - the vector

471:    Note:
472:    For a vector of dimension n, VecSet() computes
473: $     x[i] = alpha, for i=1,...,n,
474:    so that all vector entries then equal the identical
475:    scalar value, alpha.  Use the more general routine
476:    VecSetValues() to set different vector entries.

478:    You CANNOT call this after you have called VecSetValues() but before you call
479:    VecAssemblyBegin/End().

481:    Level: beginner

483: .seealso VecSetValues(), VecSetValuesBlocked(), VecSetRandom()

485: @*/
486: PetscErrorCode  VecSet(Vec x,PetscScalar alpha)
487: {
488:   PetscReal val;

494:   VecSetErrorIfLocked(x,1);

496:   PetscLogEventBegin(VEC_Set,x,0,0,0);
497:   (*x->ops->set)(x,alpha);
498:   PetscLogEventEnd(VEC_Set,x,0,0,0);
499:   PetscObjectStateIncrease((PetscObject)x);

501:   /*  norms can be simply set (if |alpha|*N not too large) */
502:   val = PetscAbsScalar(alpha);
503:   if (x->map->N == 0) {
504:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],0.0l);
505:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],0.0);
506:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],0.0);
507:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],0.0);
508:   } else if (val > PETSC_MAX_REAL/x->map->N) {
509:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
510:   } else {
511:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],x->map->N * val);
512:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
513:     val  = PetscSqrtReal((PetscReal)x->map->N) * val;
514:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],val);
515:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],val);
516:   }
517:   return 0;
518: }

520: /*@
521:    VecAXPY - Computes y = alpha x + y.

523:    Logically Collective on Vec

525:    Input Parameters:
526: +  alpha - the scalar
527: -  x, y  - the vectors

529:    Output Parameter:
530: .  y - output vector

532:    Level: intermediate

534:    Notes:
535:     x and y MUST be different vectors
536:     This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine

538: $    VecAXPY(y,alpha,x)                   y = alpha x           +      y
539: $    VecAYPX(y,beta,x)                    y =       x           + beta y
540: $    VecAXPBY(y,alpha,beta,x)             y = alpha x           + beta y
541: $    VecWAXPY(w,alpha,x,y)                w = alpha x           +      y
542: $    VecAXPBYPCZ(w,alpha,beta,gamma,x,y)  z = alpha x           + beta y + gamma z
543: $    VecMAXPY(y,nv,alpha[],x[])           y = sum alpha[i] x[i] +      y

545: .seealso:  VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPBYPCZ(), VecAXPBY()
546: @*/
547: PetscErrorCode  VecAXPY(Vec y,PetscScalar alpha,Vec x)
548: {
554:   VecCheckSameSize(x,3,y,1);
557:   if (alpha == (PetscScalar)0.0) return 0;
558:   VecSetErrorIfLocked(y,1);

560:   VecLockReadPush(x);
561:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
562:   (*y->ops->axpy)(y,alpha,x);
563:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
564:   VecLockReadPop(x);
565:   PetscObjectStateIncrease((PetscObject)y);
566:   return 0;
567: }

569: /*@
570:    VecAXPBY - Computes y = alpha x + beta y.

572:    Logically Collective on Vec

574:    Input Parameters:
575: +  alpha,beta - the scalars
576: -  x, y  - the vectors

578:    Output Parameter:
579: .  y - output vector

581:    Level: intermediate

583:    Notes:
584:     x and y MUST be different vectors
585:     The implementation is optimized for alpha and/or beta values of 0.0 and 1.0

587: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ()
588: @*/
589: PetscErrorCode  VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
590: {
596:   VecCheckSameSize(y,1,x,4);
600:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) return 0;
601:   VecSetErrorIfLocked(y,1);
602:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
603:   (*y->ops->axpby)(y,alpha,beta,x);
604:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
605:   PetscObjectStateIncrease((PetscObject)y);
606:   return 0;
607: }

609: /*@
610:    VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z

612:    Logically Collective on Vec

614:    Input Parameters:
615: +  alpha,beta, gamma - the scalars
616: -  x, y, z  - the vectors

618:    Output Parameter:
619: .  z - output vector

621:    Level: intermediate

623:    Notes:
624:     x, y and z must be different vectors
625:     The implementation is optimized for alpha of 1.0 and gamma of 1.0 or 0.0

627: .seealso:  VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBY()
628: @*/
629: PetscErrorCode  VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
630: {
639:   VecCheckSameSize(x,1,y,5);
640:   VecCheckSameSize(x,1,z,6);
646:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) return 0;
647:   VecSetErrorIfLocked(z,1);

649:   PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
650:   (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
651:   PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
652:   PetscObjectStateIncrease((PetscObject)z);
653:   return 0;
654: }

656: /*@
657:    VecAYPX - Computes y = x + beta y.

659:    Logically Collective on Vec

661:    Input Parameters:
662: +  beta - the scalar
663: -  x, y  - the vectors

665:    Output Parameter:
666: .  y - output vector

668:    Level: intermediate

670:    Notes:
671:     x and y MUST be different vectors
672:     The implementation is optimized for beta of -1.0, 0.0, and 1.0

674: .seealso:  VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
675: @*/
676: PetscErrorCode  VecAYPX(Vec y,PetscScalar beta,Vec x)
677: {
683:   VecCheckSameSize(x,1,y,3);
686:   VecSetErrorIfLocked(y,1);

688:   PetscLogEventBegin(VEC_AYPX,x,y,0,0);
689:   (*y->ops->aypx)(y,beta,x);
690:   PetscLogEventEnd(VEC_AYPX,x,y,0,0);
691:   PetscObjectStateIncrease((PetscObject)y);
692:   return 0;
693: }

695: /*@
696:    VecWAXPY - Computes w = alpha x + y.

698:    Logically Collective on Vec

700:    Input Parameters:
701: +  alpha - the scalar
702: -  x, y  - the vectors

704:    Output Parameter:
705: .  w - the result

707:    Level: intermediate

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

713: .seealso: VecAXPY(), VecAYPX(), VecAXPBY(), VecMAXPY(), VecAXPBYPCZ()
714: @*/
715: PetscErrorCode  VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
716: {
725:   VecCheckSameSize(x,3,y,4);
726:   VecCheckSameSize(x,3,w,1);
730:   VecSetErrorIfLocked(w,1);

732:   PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
733:   (*w->ops->waxpy)(w,alpha,x,y);
734:   PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
735:   PetscObjectStateIncrease((PetscObject)w);
736:   return 0;
737: }

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

742:    Not Collective

744:    Input Parameters:
745: +  x - vector to insert in
746: .  ni - number of elements to add
747: .  ix - indices where to add
748: .  y - array of values
749: -  iora - either INSERT_VALUES or ADD_VALUES, where
750:    ADD_VALUES adds values to any existing entries, and
751:    INSERT_VALUES replaces existing entries with new values

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

756:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
757:    options cannot be mixed without intervening calls to the assembly
758:    routines.

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

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

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

771:    Level: beginner

773: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
774:            VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
775: @*/
776: PetscErrorCode  VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
777: {
780:   if (!ni) return 0;

785:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
786:   (*x->ops->setvalues)(x,ni,ix,y,iora);
787:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
788:   PetscObjectStateIncrease((PetscObject)x);
789:   return 0;
790: }

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

796:     Not Collective

798:    Input Parameters:
799: +  x - vector to get values from
800: .  ni - number of elements to get
801: -  ix - indices where to get them from (in global 1d numbering)

803:    Output Parameter:
804: .   y - array of values

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

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

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

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

815:    If you call VecSetOption(x, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE),
816:    negative indices may be passed in ix. These rows are
817:    simply ignored.

819:    Level: beginner

821: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues()
822: @*/
823: PetscErrorCode  VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
824: {
826:   if (!ni) return 0;
830:   (*x->ops->getvalues)(x,ni,ix,y);
831:   return 0;
832: }

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

837:    Not Collective

839:    Input Parameters:
840: +  x - vector to insert in
841: .  ni - number of blocks to add
842: .  ix - indices where to add in block count, rather than element count
843: .  y - array of values
844: -  iora - either INSERT_VALUES or ADD_VALUES, where
845:    ADD_VALUES adds values to any existing entries, and
846:    INSERT_VALUES replaces existing entries with new values

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

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

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

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

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

866:    Level: intermediate

868: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
869:            VecSetValues()
870: @*/
871: PetscErrorCode  VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
872: {
875:   if (!ni) return 0;

880:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
881:   (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
882:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
883:   PetscObjectStateIncrease((PetscObject)x);
884:   return 0;
885: }

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

891:    Not Collective

893:    Input Parameters:
894: +  x - vector to insert in
895: .  ni - number of elements to add
896: .  ix - indices where to add
897: .  y - array of values
898: -  iora - either INSERT_VALUES or ADD_VALUES, where
899:    ADD_VALUES adds values to any existing entries, and
900:    INSERT_VALUES replaces existing entries with new values

902:    Level: intermediate

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

907:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
908:    options cannot be mixed without intervening calls to the assembly
909:    routines.

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

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

916: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetLocalToGlobalMapping(),
917:            VecSetValuesBlockedLocal()
918: @*/
919: PetscErrorCode  VecSetValuesLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
920: {
921:   PetscInt       lixp[128],*lix = lixp;

925:   if (!ni) return 0;

930:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
931:   if (!x->ops->setvalueslocal) {
932:     if (x->map->mapping) {
933:       if (ni > 128) {
934:         PetscMalloc1(ni,&lix);
935:       }
936:       ISLocalToGlobalMappingApply(x->map->mapping,ni,(PetscInt*)ix,lix);
937:       (*x->ops->setvalues)(x,ni,lix,y,iora);
938:       if (ni > 128) {
939:         PetscFree(lix);
940:       }
941:     } else {
942:       (*x->ops->setvalues)(x,ni,ix,y,iora);
943:     }
944:   } else {
945:     (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
946:   }
947:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
948:   PetscObjectStateIncrease((PetscObject)x);
949:   return 0;
950: }

952: /*@
953:    VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
954:    using a local ordering of the nodes.

956:    Not Collective

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

967:    Level: intermediate

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

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

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

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

982: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
983:            VecSetLocalToGlobalMapping()
984: @*/
985: PetscErrorCode  VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
986: {
987:   PetscInt       lixp[128],*lix = lixp;

991:   if (!ni) return 0;
995:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
996:   if (x->map->mapping) {
997:     if (ni > 128) {
998:       PetscMalloc1(ni,&lix);
999:     }
1000:     ISLocalToGlobalMappingApplyBlock(x->map->mapping,ni,(PetscInt*)ix,lix);
1001:     (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1002:     if (ni > 128) {
1003:       PetscFree(lix);
1004:     }
1005:   } else {
1006:     (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
1007:   }
1008:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1009:   PetscObjectStateIncrease((PetscObject)x);
1010:   return 0;
1011: }

1013: /*@
1014:    VecMTDot - Computes indefinite vector multiple dot products.
1015:    That is, it does NOT use the complex conjugate.

1017:    Collective on Vec

1019:    Input Parameters:
1020: +  x - one vector
1021: .  nv - number of vectors
1022: -  y - array of vectors.  Note that vectors are pointers

1024:    Output Parameter:
1025: .  val - array of the dot products

1027:    Notes for Users of Complex Numbers:
1028:    For complex vectors, VecMTDot() computes the indefinite form
1029: $      val = (x,y) = y^T x,
1030:    where y^T denotes the transpose of y.

1032:    Use VecMDot() for the inner product
1033: $      val = (x,y) = y^H x,
1034:    where y^H denotes the conjugate transpose of y.

1036:    Level: intermediate

1038: .seealso: VecMDot(), VecTDot()
1039: @*/
1040: PetscErrorCode  VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1041: {
1044:   if (!nv) return 0;
1051:   VecCheckSameSize(x,1,*y,3);

1053:   PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1054:   (*x->ops->mtdot)(x,nv,y,val);
1055:   PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1056:   return 0;
1057: }

1059: /*@
1060:    VecMDot - Computes vector multiple dot products.

1062:    Collective on Vec

1064:    Input Parameters:
1065: +  x - one vector
1066: .  nv - number of vectors
1067: -  y - array of vectors.

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

1072:    Notes for Users of Complex Numbers:
1073:    For complex vectors, VecMDot() computes
1074: $     val = (x,y) = y^H x,
1075:    where y^H denotes the conjugate transpose of y.

1077:    Use VecMTDot() for the indefinite form
1078: $     val = (x,y) = y^T x,
1079:    where y^T denotes the transpose of y.

1081:    Level: intermediate

1083: .seealso: VecMTDot(), VecDot()
1084: @*/
1085: PetscErrorCode  VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1086: {
1089:   if (!nv) return 0;
1097:   VecCheckSameSize(x,1,*y,3);

1099:   PetscLogEventBegin(VEC_MDot,x,*y,0,0);
1100:   (*x->ops->mdot)(x,nv,y,val);
1101:   PetscLogEventEnd(VEC_MDot,x,*y,0,0);
1102:   return 0;
1103: }

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

1108:    Logically Collective on Vec

1110:    Input Parameters:
1111: +  nv - number of scalars and x-vectors
1112: .  alpha - array of scalars
1113: .  y - one vector
1114: -  x - array of vectors

1116:    Level: intermediate

1118:    Notes:
1119:     y cannot be any of the x vectors

1121: .seealso:  VecAYPX(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
1122: @*/
1123: PetscErrorCode  VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1124: {
1125:   PetscInt       i;
1126:   PetscBool      nonzero;

1130:   if (!nv) return 0;
1138:   VecCheckSameSize(y,1,*x,4);
1140:   for (i=0, nonzero = PETSC_FALSE; i<nv && !nonzero; i++) nonzero = (PetscBool)(nonzero || alpha[i] != (PetscScalar)0.0);
1141:   if (!nonzero) return 0;
1142:   VecSetErrorIfLocked(y,1);
1143:   PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1144:   (*y->ops->maxpy)(y,nv,alpha,x);
1145:   PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1146:   PetscObjectStateIncrease((PetscObject)y);
1147:   return 0;
1148: }

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

1155:    Collective on X

1157:    Input Parameters:
1158: +  nx   - number of vectors to be concatenated
1159: -  X    - array containing the vectors to be concatenated in the order of concatenation

1161:    Output Parameters:
1162: +  Y    - concatenated vector
1163: -  x_is - array of index sets corresponding to the concatenated components of Y (NULL if not needed)

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

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

1175:    Level: advanced

1177: .seealso: VECNEST, VECSCATTER, VecScatterCreate()
1178: @*/
1179: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1180: {
1181:   MPI_Comm       comm;
1182:   VecType        vec_type;
1183:   Vec            Ytmp, Xtmp;
1184:   IS             *is_tmp;
1185:   PetscInt       i, shift=0, Xnl, Xng, Xbegin;


1192:   if ((*X)->ops->concatenate) {
1193:     /* use the dedicated concatenation function if available */
1194:     (*(*X)->ops->concatenate)(nx,X,Y,x_is);
1195:   } else {
1196:     /* loop over vectors and start creating IS */
1197:     comm = PetscObjectComm((PetscObject)(*X));
1198:     VecGetType(*X, &vec_type);
1199:     PetscMalloc1(nx, &is_tmp);
1200:     for (i=0; i<nx; i++) {
1201:       VecGetSize(X[i], &Xng);
1202:       VecGetLocalSize(X[i], &Xnl);
1203:       VecGetOwnershipRange(X[i], &Xbegin, NULL);
1204:       ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]);
1205:       shift += Xng;
1206:     }
1207:     /* create the concatenated vector */
1208:     VecCreate(comm, &Ytmp);
1209:     VecSetType(Ytmp, vec_type);
1210:     VecSetSizes(Ytmp, PETSC_DECIDE, shift);
1211:     VecSetUp(Ytmp);
1212:     /* copy data from X array to Y and return */
1213:     for (i=0; i<nx; i++) {
1214:       VecGetSubVector(Ytmp, is_tmp[i], &Xtmp);
1215:       VecCopy(X[i], Xtmp);
1216:       VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp);
1217:     }
1218:     *Y = Ytmp;
1219:     if (x_is) {
1220:       *x_is = is_tmp;
1221:     } else {
1222:       for (i=0; i<nx; i++) {
1223:         ISDestroy(&is_tmp[i]);
1224:       }
1225:       PetscFree(is_tmp);
1226:     }
1227:   }
1228:   return 0;
1229: }

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

1234:     Input Parameters:
1235: +   X - the original vector
1236: -   is - the index set of the subvector

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

1243: */
1244: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X,IS is,PetscBool *contig,PetscInt *start,PetscInt *blocksize)
1245: {
1246:   PetscInt         gstart,gend,lstart;
1247:   PetscBool        red[2] = {PETSC_TRUE/*contiguous*/,PETSC_TRUE/*validVBS*/};
1248:   PetscInt         n,N,ibs,vbs,bs = -1;

1250:   ISGetLocalSize(is,&n);
1251:   ISGetSize(is,&N);
1252:   ISGetBlockSize(is,&ibs);
1253:   VecGetBlockSize(X,&vbs);
1254:   VecGetOwnershipRange(X,&gstart,&gend);
1255:   ISContiguousLocal(is,gstart,gend,&lstart,&red[0]);
1256:   /* block size is given by IS if ibs > 1; otherwise, check the vector */
1257:   if (ibs > 1) {
1258:     MPIU_Allreduce(MPI_IN_PLACE,red,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1259:     bs   = ibs;
1260:   } else {
1261:     if (n%vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1262:     MPIU_Allreduce(MPI_IN_PLACE,red,2,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1263:     if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1264:   }

1266:   *contig     = red[0];
1267:   *start      = lstart;
1268:   *blocksize  = bs;
1269:   return 0;
1270: }

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

1274:     Input Parameters:
1275: +   X - the original vector
1276: .   is - the index set of the subvector
1277: -   bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()

1279:     Output Parameters:
1280: .   Z  - the subvector, which will compose the VecScatter context on output
1281: */
1282: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X,IS is,PetscInt bs,Vec *Z)
1283: {
1284:   PetscInt       n,N;
1285:   VecScatter     vscat;
1286:   Vec            Y;

1288:   ISGetLocalSize(is,&n);
1289:   ISGetSize(is,&N);
1290:   VecCreate(PetscObjectComm((PetscObject)is),&Y);
1291:   VecSetSizes(Y,n,N);
1292:   VecSetBlockSize(Y,bs);
1293:   VecSetType(Y,((PetscObject)X)->type_name);
1294:   VecScatterCreate(X,is,Y,NULL,&vscat);
1295:   VecScatterBegin(vscat,X,Y,INSERT_VALUES,SCATTER_FORWARD);
1296:   VecScatterEnd(vscat,X,Y,INSERT_VALUES,SCATTER_FORWARD);
1297:   PetscObjectCompose((PetscObject)Y,"VecGetSubVector_Scatter",(PetscObject)vscat);
1298:   VecScatterDestroy(&vscat);
1299:   *Z   = Y;
1300:   return 0;
1301: }

1303: /*@
1304:    VecGetSubVector - Gets a vector representing part of another vector

1306:    Collective on X and IS

1308:    Input Parameters:
1309: + X - vector from which to extract a subvector
1310: - is - index set representing portion of X to extract

1312:    Output Parameter:
1313: . Y - subvector corresponding to is

1315:    Level: advanced

1317:    Notes:
1318:    The subvector Y should be returned with VecRestoreSubVector().
1319:    X and is must be defined on the same communicator

1321:    This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1322:    modifying the subvector.  Other non-overlapping subvectors can still be obtained from X using this function.
1323:    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.

1325: .seealso: MatCreateSubMatrix()
1326: @*/
1327: PetscErrorCode  VecGetSubVector(Vec X,IS is,Vec *Y)
1328: {
1329:   Vec              Z;

1335:   if (X->ops->getsubvector) {
1336:     (*X->ops->getsubvector)(X,is,&Z);
1337:   } else { /* Default implementation currently does no caching */
1338:     PetscBool   contig;
1339:     PetscInt    n,N,start,bs;

1341:     ISGetLocalSize(is,&n);
1342:     ISGetSize(is,&N);
1343:     VecGetSubVectorContiguityAndBS_Private(X,is,&contig,&start,&bs);
1344:     if (contig) { /* We can do a no-copy implementation */
1345:       const PetscScalar *x;
1346:       PetscInt          state = 0;
1347:       PetscBool         isstd,iscuda,iship;

1349:       PetscObjectTypeCompareAny((PetscObject)X,&isstd,VECSEQ,VECMPI,VECSTANDARD,"");
1350:       PetscObjectTypeCompareAny((PetscObject)X,&iscuda,VECSEQCUDA,VECMPICUDA,"");
1351:       PetscObjectTypeCompareAny((PetscObject)X,&iship,VECSEQHIP,VECMPIHIP,"");
1352:       if (iscuda) {
1353: #if defined(PETSC_HAVE_CUDA)
1354:         const PetscScalar *x_d;
1355:         PetscMPIInt       size;
1356:         PetscOffloadMask  flg;

1358:         VecCUDAGetArrays_Private(X,&x,&x_d,&flg);
1361:         if (x) x += start;
1362:         if (x_d) x_d += start;
1363:         MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1364:         if (size == 1) {
1365:           VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X),bs,n,x,x_d,&Z);
1366:         } else {
1367:           VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X),bs,n,N,x,x_d,&Z);
1368:         }
1369:         Z->offloadmask = flg;
1370: #endif
1371:       } else if (iship) {
1372: #if defined(PETSC_HAVE_HIP)
1373:         const PetscScalar *x_d;
1374:         PetscMPIInt       size;
1375:         PetscOffloadMask  flg;

1377:         VecHIPGetArrays_Private(X,&x,&x_d,&flg);
1380:         if (x) x += start;
1381:         if (x_d) x_d += start;
1382:         MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1383:         if (size == 1) {
1384:           VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X),bs,n,x,x_d,&Z);
1385:         } else {
1386:           VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X),bs,n,N,x,x_d,&Z);
1387:         }
1388:         Z->offloadmask = flg;
1389: #endif
1390:       } else if (isstd) {
1391:         PetscMPIInt size;

1393:         MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1394:         VecGetArrayRead(X,&x);
1395:         if (x) x += start;
1396:         if (size == 1) {
1397:           VecCreateSeqWithArray(PetscObjectComm((PetscObject)X),bs,n,x,&Z);
1398:         } else {
1399:           VecCreateMPIWithArray(PetscObjectComm((PetscObject)X),bs,n,N,x,&Z);
1400:         }
1401:         VecRestoreArrayRead(X,&x);
1402:       } else { /* default implementation: use place array */
1403:         VecGetArrayRead(X,&x);
1404:         VecCreate(PetscObjectComm((PetscObject)X),&Z);
1405:         VecSetType(Z,((PetscObject)X)->type_name);
1406:         VecSetSizes(Z,n,N);
1407:         VecSetBlockSize(Z,bs);
1408:         VecPlaceArray(Z,x ? x+start : NULL);
1409:         VecRestoreArrayRead(X,&x);
1410:       }

1412:       /* this is relevant only in debug mode */
1413:       VecLockGet(X,&state);
1414:       if (state) {
1415:         VecLockReadPush(Z);
1416:       }
1417:       Z->ops->placearray = NULL;
1418:       Z->ops->replacearray = NULL;
1419:     } else { /* Have to create a scatter and do a copy */
1420:       VecGetSubVectorThroughVecScatter_Private(X,is,bs,&Z);
1421:     }
1422:   }
1423:   /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1424:   if (VecGetSubVectorSavedStateId < 0) PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);
1425:   PetscObjectComposedDataSetInt((PetscObject)Z,VecGetSubVectorSavedStateId,1);
1426:   *Y   = Z;
1427:   return 0;
1428: }

1430: /*@
1431:    VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()

1433:    Collective on IS

1435:    Input Parameters:
1436: + X - vector from which subvector was obtained
1437: . is - index set representing the subset of X
1438: - Y - subvector being restored

1440:    Level: advanced

1442: .seealso: VecGetSubVector()
1443: @*/
1444: PetscErrorCode  VecRestoreSubVector(Vec X,IS is,Vec *Y)
1445: {
1446:   PETSC_UNUSED PetscObjectState dummystate = 0;
1447:   PetscBool                     unchanged;


1455:   if (X->ops->restoresubvector) {
1456:     (*X->ops->restoresubvector)(X,is,Y);
1457:   } else {
1458:     PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,dummystate,unchanged);
1459:     if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1460:       VecScatter scatter;
1461:       PetscInt   state;

1463:       VecLockGet(X,&state);

1466:       PetscObjectQuery((PetscObject)*Y,"VecGetSubVector_Scatter",(PetscObject*)&scatter);
1467:       if (scatter) {
1468:         VecScatterBegin(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1469:         VecScatterEnd(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1470:       } else {
1471:         PetscBool         iscuda,iship;
1472:         PetscObjectTypeCompareAny((PetscObject)X,&iscuda,VECSEQCUDA,VECMPICUDA,"");
1473:         PetscObjectTypeCompareAny((PetscObject)X,&iship,VECSEQHIP,VECMPIHIP,"");

1475:         if (iscuda) {
1476: #if defined(PETSC_HAVE_CUDA)
1477:           PetscOffloadMask ymask = (*Y)->offloadmask;

1479:           /* The offloadmask of X dictates where to move memory
1480:               If X GPU data is valid, then move Y data on GPU if needed
1481:               Otherwise, move back to the CPU */
1482:           switch (X->offloadmask) {
1483:           case PETSC_OFFLOAD_BOTH:
1484:             if (ymask == PETSC_OFFLOAD_CPU) {
1485:               VecCUDAResetArray(*Y);
1486:             } else if (ymask == PETSC_OFFLOAD_GPU) {
1487:               X->offloadmask = PETSC_OFFLOAD_GPU;
1488:             }
1489:             break;
1490:           case PETSC_OFFLOAD_GPU:
1491:             if (ymask == PETSC_OFFLOAD_CPU) {
1492:               VecCUDAResetArray(*Y);
1493:             }
1494:             break;
1495:           case PETSC_OFFLOAD_CPU:
1496:             if (ymask == PETSC_OFFLOAD_GPU) {
1497:               VecResetArray(*Y);
1498:             }
1499:             break;
1500:           case PETSC_OFFLOAD_UNALLOCATED:
1501:           case PETSC_OFFLOAD_KOKKOS:
1502:             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"This should not happen");
1503:           }
1504: #endif
1505:         } else if (iship) {
1506: #if defined(PETSC_HAVE_HIP)
1507:           PetscOffloadMask ymask = (*Y)->offloadmask;

1509:           /* The offloadmask of X dictates where to move memory
1510:               If X GPU data is valid, then move Y data on GPU if needed
1511:               Otherwise, move back to the CPU */
1512:           switch (X->offloadmask) {
1513:           case PETSC_OFFLOAD_BOTH:
1514:             if (ymask == PETSC_OFFLOAD_CPU) {
1515:               VecHIPResetArray(*Y);
1516:             } else if (ymask == PETSC_OFFLOAD_GPU) {
1517:               X->offloadmask = PETSC_OFFLOAD_GPU;
1518:             }
1519:             break;
1520:           case PETSC_OFFLOAD_GPU:
1521:             if (ymask == PETSC_OFFLOAD_CPU) {
1522:               VecHIPResetArray(*Y);
1523:             }
1524:             break;
1525:           case PETSC_OFFLOAD_CPU:
1526:             if (ymask == PETSC_OFFLOAD_GPU) {
1527:               VecResetArray(*Y);
1528:             }
1529:             break;
1530:           case PETSC_OFFLOAD_UNALLOCATED:
1531:           case PETSC_OFFLOAD_KOKKOS:
1532:             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"This should not happen");
1533:           }
1534: #endif
1535:         } else {
1536:           /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1537:           VecResetArray(*Y);
1538:         }
1539:         PetscObjectStateIncrease((PetscObject)X);
1540:       }
1541:     }
1542:   }
1543:   VecDestroy(Y);
1544:   return 0;
1545: }

1547: /*@
1548:    VecGetLocalVectorRead - Maps the local portion of a vector into a
1549:    vector.  You must call VecRestoreLocalVectorRead() when the local
1550:    vector is no longer needed.

1552:    Not collective.

1554:    Input parameter:
1555: .  v - The vector for which the local vector is desired.

1557:    Output parameter:
1558: .  w - Upon exit this contains the local vector.

1560:    Level: beginner

1562:    Notes:
1563:    This function is similar to VecGetArrayRead() which maps the local
1564:    portion into a raw pointer.  VecGetLocalVectorRead() is usually
1565:    almost as efficient as VecGetArrayRead() but in certain circumstances
1566:    VecGetLocalVectorRead() can be much more efficient than
1567:    VecGetArrayRead().  This is because the construction of a contiguous
1568:    array representing the vector data required by VecGetArrayRead() can
1569:    be an expensive operation for certain vector types.  For example, for
1570:    GPU vectors VecGetArrayRead() requires that the data between device
1571:    and host is synchronized.

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

1576: .seealso: VecRestoreLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1577: @*/
1578: PetscErrorCode VecGetLocalVectorRead(Vec v,Vec w)
1579: {
1580:   PetscScalar    *a;

1584:   VecCheckSameLocalSize(v,1,w,2);
1585:   if (v->ops->getlocalvectorread) {
1586:     (*v->ops->getlocalvectorread)(v,w);
1587:   } else {
1588:     VecGetArrayRead(v,(const PetscScalar**)&a);
1589:     VecPlaceArray(w,a);
1590:   }
1591:   PetscObjectStateIncrease((PetscObject)w);
1592:   VecLockReadPush(v);
1593:   VecLockReadPush(w);
1594:   return 0;
1595: }

1597: /*@
1598:    VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1599:    previously mapped into a vector using VecGetLocalVectorRead().

1601:    Not collective.

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

1607:    Level: beginner

1609: .seealso: VecGetLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1610: @*/
1611: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec w)
1612: {
1613:   PetscScalar    *a;

1617:   if (v->ops->restorelocalvectorread) {
1618:     (*v->ops->restorelocalvectorread)(v,w);
1619:   } else {
1620:     VecGetArrayRead(w,(const PetscScalar**)&a);
1621:     VecRestoreArrayRead(v,(const PetscScalar**)&a);
1622:     VecResetArray(w);
1623:   }
1624:   VecLockReadPop(v);
1625:   VecLockReadPop(w);
1626:   PetscObjectStateIncrease((PetscObject)w);
1627:   return 0;
1628: }

1630: /*@
1631:    VecGetLocalVector - Maps the local portion of a vector into a
1632:    vector.

1634:    Collective on v, not collective on w.

1636:    Input parameter:
1637: .  v - The vector for which the local vector is desired.

1639:    Output parameter:
1640: .  w - Upon exit this contains the local vector.

1642:    Level: beginner

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

1655: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1656: @*/
1657: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1658: {
1659:   PetscScalar    *a;

1663:   VecCheckSameLocalSize(v,1,w,2);
1664:   if (v->ops->getlocalvector) {
1665:     (*v->ops->getlocalvector)(v,w);
1666:   } else {
1667:     VecGetArray(v,&a);
1668:     VecPlaceArray(w,a);
1669:   }
1670:   PetscObjectStateIncrease((PetscObject)w);
1671:   return 0;
1672: }

1674: /*@
1675:    VecRestoreLocalVector - Unmaps the local portion of a vector
1676:    previously mapped into a vector using VecGetLocalVector().

1678:    Logically collective.

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

1684:    Level: beginner

1686: .seealso: VecGetLocalVector(), VecGetLocalVectorRead(), VecRestoreLocalVectorRead(), LocalVectorRead(), VecGetArrayRead(), VecGetArray()
1687: @*/
1688: PetscErrorCode VecRestoreLocalVector(Vec v,Vec w)
1689: {
1690:   PetscScalar    *a;

1694:   if (v->ops->restorelocalvector) {
1695:     (*v->ops->restorelocalvector)(v,w);
1696:   } else {
1697:     VecGetArray(w,&a);
1698:     VecRestoreArray(v,&a);
1699:     VecResetArray(w);
1700:   }
1701:   PetscObjectStateIncrease((PetscObject)w);
1702:   PetscObjectStateIncrease((PetscObject)v);
1703:   return 0;
1704: }

1706: /*@C
1707:    VecGetArray - Returns a pointer to a contiguous array that contains this
1708:    processor's portion of the vector data. For the standard PETSc
1709:    vectors, VecGetArray() returns a pointer to the local data array and
1710:    does not use any copies. If the underlying vector data is not stored
1711:    in a contiguous array this routine will copy the data to a contiguous
1712:    array and return a pointer to that. You MUST call VecRestoreArray()
1713:    when you no longer need access to the array.

1715:    Logically Collective on Vec

1717:    Input Parameter:
1718: .  x - the vector

1720:    Output Parameter:
1721: .  a - location to put pointer to the array

1723:    Fortran Note:
1724:    This routine is used differently from Fortran 77
1725: $    Vec         x
1726: $    PetscScalar x_array(1)
1727: $    PetscOffset i_x
1728: $    PetscErrorCode ierr
1729: $       call VecGetArray(x,x_array,i_x,ierr)
1730: $
1731: $   Access first local entry in vector with
1732: $      value = x_array(i_x + 1)
1733: $
1734: $      ...... other code
1735: $       call VecRestoreArray(x,x_array,i_x,ierr)
1736:    For Fortran 90 see VecGetArrayF90()

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

1741:    Level: beginner

1743: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1744:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1745: @*/
1746: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1747: {
1749:   VecSetErrorIfLocked(x,1);
1750:   if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
1751:     (*x->ops->getarray)(x,a);
1752:   } else if (x->petscnative) { /* VECSTANDARD */
1753:     *a = *((PetscScalar**)x->data);
1754:   } else SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array for vector type \"%s\"",((PetscObject)x)->type_name);
1755:   return 0;
1756: }

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

1761:    Logically Collective on Vec

1763:    Input Parameters:
1764: +  x - the vector
1765: -  a - location of pointer to array obtained from VecGetArray()

1767:    Level: beginner

1769: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1770:           VecGetArrayPair(), VecRestoreArrayPair()
1771: @*/
1772: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1773: {
1775:   if (x->ops->restorearray) { /* VECNEST, VECCUDA etc */
1776:     (*x->ops->restorearray)(x,a);
1777:   } else if (x->petscnative) { /* VECSTANDARD */
1778:     /* nothing */
1779:   } else SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot restore array for vector type \"%s\"",((PetscObject)x)->type_name);
1780:   if (a) *a = NULL;
1781:   PetscObjectStateIncrease((PetscObject)x);
1782:   return 0;
1783: }
1784: /*@C
1785:    VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.

1787:    Not Collective

1789:    Input Parameter:
1790: .  x - the vector

1792:    Output Parameter:
1793: .  a - the array

1795:    Level: beginner

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

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

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

1806: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1807: @*/
1808: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1809: {
1811:   if (x->ops->getarray) { /* VECNEST, VECCUDA, VECKOKKOS etc */
1812:     (*x->ops->getarray)(x,(PetscScalar**)a);
1813:   } else if (x->petscnative) { /* VECSTANDARD */
1814:     *a = *((PetscScalar**)x->data);
1815:   } else SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array read for vector type \"%s\"",((PetscObject)x)->type_name);
1816:   return 0;
1817: }

1819: /*@C
1820:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

1822:    Not Collective

1824:    Input Parameters:
1825: +  vec - the vector
1826: -  array - the array

1828:    Level: beginner

1830: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1831: @*/
1832: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1833: {
1835:   if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
1836:     /* nothing */
1837:   } else if (x->ops->restorearrayread) { /* VECNEST */
1838:     (*x->ops->restorearrayread)(x,a);
1839:   } else { /* No one? */
1840:     (*x->ops->restorearray)(x,(PetscScalar**)a);
1841:   }
1842:   if (a) *a = NULL;
1843:   return 0;
1844: }

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

1851:    Logically Collective on Vec

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

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

1859:    Level: intermediate

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

1864: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1865:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArray(), VecRestoreArrayWrite()
1866: @*/
1867: PetscErrorCode VecGetArrayWrite(Vec x,PetscScalar **a)
1868: {
1870:   VecSetErrorIfLocked(x,1);
1871:   if (x->ops->getarraywrite) {
1872:     (*x->ops->getarraywrite)(x,a);
1873:   } else {
1874:     VecGetArray(x,a);
1875:   }
1876:   return 0;
1877: }

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

1882:    Logically Collective on Vec

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

1888:    Level: beginner

1890: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1891:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite()
1892: @*/
1893: PetscErrorCode VecRestoreArrayWrite(Vec x,PetscScalar **a)
1894: {
1896:   if (x->ops->restorearraywrite) {
1897:     (*x->ops->restorearraywrite)(x,a);
1898:   } else if (x->ops->restorearray) {
1899:     (*x->ops->restorearray)(x,a);
1900:   }
1901:   if (a) *a = NULL;
1902:   PetscObjectStateIncrease((PetscObject)x);
1903:   return 0;
1904: }

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

1911:    Logically Collective on Vec

1913:    Input Parameters:
1914: +  x - the vectors
1915: -  n - the number of vectors

1917:    Output Parameter:
1918: .  a - location to put pointer to the array

1920:    Fortran Note:
1921:    This routine is not supported in Fortran.

1923:    Level: intermediate

1925: .seealso: VecGetArray(), VecRestoreArrays()
1926: @*/
1927: PetscErrorCode  VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1928: {
1929:   PetscInt       i;
1930:   PetscScalar    **q;

1936:   PetscMalloc1(n,&q);
1937:   for (i=0; i<n; ++i) {
1938:     VecGetArray(x[i],&q[i]);
1939:   }
1940:   *a = q;
1941:   return 0;
1942: }

1944: /*@C
1945:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1946:    has been called.

1948:    Logically Collective on Vec

1950:    Input Parameters:
1951: +  x - the vector
1952: .  n - the number of vectors
1953: -  a - location of pointer to arrays obtained from VecGetArrays()

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

1961:    Fortran Note:
1962:    This routine is not supported in Fortran.

1964:    Level: intermediate

1966: .seealso: VecGetArrays(), VecRestoreArray()
1967: @*/
1968: PetscErrorCode  VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1969: {
1970:   PetscInt       i;
1971:   PetscScalar    **q = *a;


1977:   for (i=0; i<n; ++i) {
1978:     VecRestoreArray(x[i],&q[i]);
1979:   }
1980:   PetscFree(q);
1981:   return 0;
1982: }

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

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

1992:    Logically Collective on Vec

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

1997:    Output Parameters:
1998: +  a - location to put pointer to the array
1999: -  mtype - memory type of the array

2001:    Level: beginner

2003: .seealso: VecRestoreArrayAndMemType(), VecGetArrayReadAndMemType(), VecGetArrayWriteAndMemType(), VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(),
2004:           VecPlaceArray(), VecGetArray2d(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
2005: @*/
2006: PetscErrorCode VecGetArrayAndMemType(Vec x,PetscScalar **a,PetscMemType *mtype)
2007: {
2008:   PetscMemType   omtype;

2012:   VecSetErrorIfLocked(x,1);
2013:   if (x->ops->getarrayandmemtype) { /* VECCUDA, VECKOKKOS etc */
2014:     (*x->ops->getarrayandmemtype)(x,a,&omtype);
2015:   } else { /* VECSTANDARD, VECNEST, VECVIENNACL */
2016:     VecGetArray(x,a);
2017:     omtype = PETSC_MEMTYPE_HOST;
2018:   }
2019:   if (mtype) *mtype = omtype;
2020:   return 0;
2021: }

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

2026:    Logically Collective on Vec

2028:    Input Parameters:
2029: +  x - the vector
2030: -  a - location of pointer to array obtained from VecGetArrayAndMemType()

2032:    Level: beginner

2034: .seealso: VecGetArrayAndMemType(), VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(),
2035:           VecPlaceArray(), VecRestoreArray2d(), VecGetArrayPair(), VecRestoreArrayPair()
2036: @*/
2037: PetscErrorCode VecRestoreArrayAndMemType(Vec x,PetscScalar **a)
2038: {
2041:   if (x->ops->restorearrayandmemtype) { /* VECCUDA, VECKOKKOS etc */
2042:     (*x->ops->restorearrayandmemtype)(x,a);
2043:   } else if (x->ops->restorearray) { /* VECNEST, VECVIENNACL */
2044:     (*x->ops->restorearray)(x,a);
2045:   } /* VECSTANDARD does nothing */
2046:   if (a) *a = NULL;
2047:   PetscObjectStateIncrease((PetscObject)x);
2048:   return 0;
2049: }

2051: /*@C
2052:    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().

2054:    Not Collective

2056:    Input Parameter:
2057: .  x - the vector

2059:    Output Parameters:
2060: +  a - the array
2061: -  mtype - memory type of the array

2063:    Level: beginner

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

2068: .seealso: VecRestoreArrayReadAndMemType(), VecGetArrayAndMemType(), VecGetArrayWriteAndMemType(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayAndMemType()
2069: @*/
2070: PetscErrorCode VecGetArrayReadAndMemType(Vec x,const PetscScalar **a,PetscMemType *mtype)
2071: {
2072:   PetscMemType   omtype;

2076:  #if defined(PETSC_USE_DEBUG)
2078:  #endif

2080:   if (x->ops->getarrayandmemtype) { /* VECCUDA, VECKOKKOS etc, though they are also petscnative */
2081:     (*x->ops->getarrayandmemtype)(x,(PetscScalar**)a,&omtype);
2082:   } else if (x->ops->getarray) { /* VECNEST, VECVIENNACL */
2083:     (*x->ops->getarray)(x,(PetscScalar**)a);
2084:     omtype = PETSC_MEMTYPE_HOST;
2085:   } else if (x->petscnative) { /* VECSTANDARD */
2086:     *a = *((PetscScalar**)x->data);
2087:     omtype = PETSC_MEMTYPE_HOST;
2088:   } else SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array read in place for vector type \"%s\"",((PetscObject)x)->type_name);
2089:   if (mtype) *mtype = omtype;
2090:   return 0;
2091: }

2093: /*@C
2094:    VecRestoreArrayReadAndMemType - Restore array obtained with VecGetArrayReadAndMemType()

2096:    Not Collective

2098:    Input Parameters:
2099: +  vec - the vector
2100: -  array - the array

2102:    Level: beginner

2104: .seealso: VecGetArrayReadAndMemType(), VecRestoreArrayAndMemType(), VecRestoreArrayWriteAndMemType(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
2105: @*/
2106: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x,const PetscScalar **a)
2107: {
2110:   if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS, VECVIENNACL etc */
2111:     /* nothing */
2112:   } else if (x->ops->restorearrayread) { /* VECNEST */
2113:     (*x->ops->restorearrayread)(x,a);
2114:   } else SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot restore array read in place for vector type \"%s\"",((PetscObject)x)->type_name);
2115:   if (a) *a = NULL;
2116:   return 0;
2117: }

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

2123:    Not Collective

2125:    Input Parameter:
2126: .  x - the vector

2128:    Output Parameters:
2129: +  a - the array
2130: -  mtype - memory type of the array

2132:    Level: beginner

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

2137: .seealso: VecRestoreArrayWriteAndMemType(), VecGetArrayReadAndMemType(), VecGetArrayAndMemType(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair(),
2138: @*/
2139: PetscErrorCode VecGetArrayWriteAndMemType(Vec x,PetscScalar **a,PetscMemType *mtype)
2140: {
2141:   PetscMemType   omtype;

2145:   if (x->ops->getarraywriteandmemtype) { /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2146:     (*x->ops->getarrayandmemtype)(x,a,&omtype);
2147:   } else if (x->ops->getarraywrite) { /* VECNEST, VECVIENNACL */
2148:     (*x->ops->getarraywrite)(x,a);
2149:     omtype = PETSC_MEMTYPE_HOST;
2150:   } else if (x->petscnative) { /* VECSTANDARD */
2151:     *a = *((PetscScalar**)x->data);
2152:     omtype = PETSC_MEMTYPE_HOST;
2153:   } else SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array read in place for vector type \"%s\"",((PetscObject)x)->type_name);
2154:   if (mtype) *mtype = omtype;
2155:   return 0;
2156: }

2158: /*@C
2159:    VecRestoreArrayWriteAndMemType - Restore array obtained with VecGetArrayWriteAndMemType()

2161:    Not Collective

2163:    Input Parameters:
2164: +  vec - the vector
2165: -  array - the array

2167:    Level: beginner

2169: .seealso: VecGetArrayWriteAndMemType(), VecRestoreArrayAndMemType(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
2170: @*/
2171: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x,PetscScalar **a)
2172: {
2173:   VecRestoreArrayAndMemType(x,a);
2174:   return 0;
2175: }

2177: /*@
2178:    VecPlaceArray - Allows one to replace the array in a vector with an
2179:    array provided by the user. This is useful to avoid copying an array
2180:    into a vector.

2182:    Not Collective

2184:    Input Parameters:
2185: +  vec - the vector
2186: -  array - the array

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

2194:    Level: developer

2196: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecResetArray()

2198: @*/
2199: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
2200: {
2204:   if (vec->ops->placearray) {
2205:     (*vec->ops->placearray)(vec,array);
2206:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
2207:   PetscObjectStateIncrease((PetscObject)vec);
2208:   return 0;
2209: }

2211: /*@C
2212:    VecReplaceArray - Allows one to replace the array in a vector with an
2213:    array provided by the user. This is useful to avoid copying an array
2214:    into a vector.

2216:    Not Collective

2218:    Input Parameters:
2219: +  vec - the vector
2220: -  array - the array

2222:    Notes:
2223:    This permanently replaces the array and frees the memory associated
2224:    with the old array.

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

2229:    Not supported from Fortran

2231:    Level: developer

2233: .seealso: VecGetArray(), VecRestoreArray(), VecPlaceArray(), VecResetArray()

2235: @*/
2236: PetscErrorCode  VecReplaceArray(Vec vec,const PetscScalar array[])
2237: {
2240:   if (vec->ops->replacearray) {
2241:     (*vec->ops->replacearray)(vec,array);
2242:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
2243:   PetscObjectStateIncrease((PetscObject)vec);
2244:   return 0;
2245: }

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

2250:    This function has semantics similar to VecGetArray():  the pointer
2251:    returned by this function points to a consistent view of the vector
2252:    data.  This may involve a copy operation of data from the host to the
2253:    device if the data on the device is out of date.  If the device
2254:    memory hasn't been allocated previously it will be allocated as part
2255:    of this function call.  VecCUDAGetArray() assumes that
2256:    the user will modify the vector data.  This is similar to
2257:    intent(inout) in fortran.

2259:    The CUDA device pointer has to be released by calling
2260:    VecCUDARestoreArray().  Upon restoring the vector data
2261:    the data on the host will be marked as out of date.  A subsequent
2262:    access of the host data will thus incur a data transfer from the
2263:    device to the host.

2265:    Input Parameter:
2266: .  v - the vector

2268:    Output Parameter:
2269: .  a - the CUDA device pointer

2271:    Fortran note:
2272:    This function is not currently available from Fortran.

2274:    Level: intermediate

2276: .seealso: VecCUDARestoreArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2277: @*/
2278: PETSC_EXTERN PetscErrorCode VecCUDAGetArray(Vec v, PetscScalar **a)
2279: {
2281:  #if defined(PETSC_HAVE_CUDA)
2282:   {
2283:     VecCUDACopyToGPU(v);
2284:     *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2285:   }
2286:  #endif
2287:   return 0;
2288: }

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

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

2297:    Input Parameters:
2298: +  v - the vector
2299: -  a - the CUDA device pointer.  This pointer is invalid after
2300:        VecCUDARestoreArray() returns.

2302:    Fortran note:
2303:    This function is not currently available from Fortran.

2305:    Level: intermediate

2307: .seealso: VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2308: @*/
2309: PETSC_EXTERN PetscErrorCode VecCUDARestoreArray(Vec v, PetscScalar **a)
2310: {
2312: #if defined(PETSC_HAVE_CUDA)
2313:   v->offloadmask = PETSC_OFFLOAD_GPU;
2314: #endif
2315:   PetscObjectStateIncrease((PetscObject)v);
2316:   return 0;
2317: }

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

2322:    This function is analogous to VecGetArrayRead():  The pointer
2323:    returned by this function points to a consistent view of the vector
2324:    data.  This may involve a copy operation of data from the host to the
2325:    device if the data on the device is out of date.  If the device
2326:    memory hasn't been allocated previously it will be allocated as part
2327:    of this function call.  VecCUDAGetArrayRead() assumes that the
2328:    user will not modify the vector data.  This is analgogous to
2329:    intent(in) in Fortran.

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

2337:    Input Parameter:
2338: .  v - the vector

2340:    Output Parameter:
2341: .  a - the CUDA pointer.

2343:    Fortran note:
2344:    This function is not currently available from Fortran.

2346:    Level: intermediate

2348: .seealso: VecCUDARestoreArrayRead(), VecCUDAGetArray(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2349: @*/
2350: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayRead(Vec v,const PetscScalar** a)
2351: {
2352:    VecCUDAGetArray(v,(PetscScalar**)a);
2353:    return 0;
2354: }

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

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

2364:    Input Parameters:
2365: +  v - the vector
2366: -  a - the CUDA device pointer.  This pointer is invalid after
2367:        VecCUDARestoreArrayRead() returns.

2369:    Fortran note:
2370:    This function is not currently available from Fortran.

2372:    Level: intermediate

2374: .seealso: VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2375: @*/
2376: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayRead(Vec v, const PetscScalar **a)
2377: {
2379:   *a = NULL;
2380:   return 0;
2381: }

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

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

2392:    The device pointer needs to be released with
2393:    VecCUDARestoreArrayWrite().  When the pointer is released the
2394:    host data of the vector is marked as out of data.  Subsequent access
2395:    of the host data with e.g. VecGetArray() incurs a device to host data
2396:    transfer.

2398:    Input Parameter:
2399: .  v - the vector

2401:    Output Parameter:
2402: .  a - the CUDA pointer

2404:    Fortran note:
2405:    This function is not currently available from Fortran.

2407:    Level: advanced

2409: .seealso: VecCUDARestoreArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2410: @*/
2411: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayWrite(Vec v, PetscScalar **a)
2412: {
2414:  #if defined(PETSC_HAVE_CUDA)
2415:   {
2416:     VecCUDAAllocateCheck(v);
2417:     *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2418:   }
2419:  #endif
2420:   return 0;
2421: }

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

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

2430:    Input Parameters:
2431: +  v - the vector
2432: -  a - the CUDA device pointer.  This pointer is invalid after
2433:        VecCUDARestoreArrayWrite() returns.

2435:    Fortran note:
2436:    This function is not currently available from Fortran.

2438:    Level: intermediate

2440: .seealso: VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2441: @*/
2442: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayWrite(Vec v, PetscScalar **a)
2443: {
2445:  #if defined(PETSC_HAVE_CUDA)
2446:   v->offloadmask = PETSC_OFFLOAD_GPU;
2447:   if (a) *a = NULL;
2448:  #endif
2449:   PetscObjectStateIncrease((PetscObject)v);
2450:   return 0;
2451: }

2453: /*@C
2454:    VecCUDAPlaceArray - Allows one to replace the GPU array in a vector with a
2455:    GPU array provided by the user. This is useful to avoid copying an
2456:    array into a vector.

2458:    Not Collective

2460:    Input Parameters:
2461: +  vec - the vector
2462: -  array - the GPU array

2464:    Notes:
2465:    You can return to the original GPU array with a call to VecCUDAResetArray()
2466:    It is not possible to use VecCUDAPlaceArray() and VecPlaceArray() at the
2467:    same time on the same vector.

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

2474:    Level: developer

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

2478: @*/
2479: PetscErrorCode VecCUDAPlaceArray(Vec vin,const PetscScalar a[])
2480: {
2482: #if defined(PETSC_HAVE_CUDA)
2483:   VecCUDACopyToGPU(vin);
2485:   ((Vec_Seq*)vin->data)->unplacedarray  = (PetscScalar *) ((Vec_CUDA*)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2486:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2487:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2488: #endif
2489:   PetscObjectStateIncrease((PetscObject)vin);
2490:   return 0;
2491: }

2493: /*@C
2494:    VecCUDAReplaceArray - Allows one to replace the GPU array in a vector
2495:    with a GPU array provided by the user. This is useful to avoid copying
2496:    a GPU array into a vector.

2498:    Not Collective

2500:    Input Parameters:
2501: +  vec - the vector
2502: -  array - the GPU array

2504:    Notes:
2505:    This permanently replaces the GPU array and frees the memory associated
2506:    with the old GPU array.

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

2511:    Not supported from Fortran

2513:    Level: developer

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

2517: @*/
2518: PetscErrorCode VecCUDAReplaceArray(Vec vin,const PetscScalar a[])
2519: {
2520: #if defined(PETSC_HAVE_CUDA)
2521: #endif

2524: #if defined(PETSC_HAVE_CUDA)
2525:   if (((Vec_CUDA*)vin->spptr)->GPUarray_allocated) {
2526:     cudaFree(((Vec_CUDA*)vin->spptr)->GPUarray_allocated);
2527:   }
2528:   ((Vec_CUDA*)vin->spptr)->GPUarray_allocated = ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2529:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2530: #endif
2531:   PetscObjectStateIncrease((PetscObject)vin);
2532:   return 0;
2533: }

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

2539:    Not Collective

2541:    Input Parameters:
2542: .  vec - the vector

2544:    Level: developer

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

2548: @*/
2549: PetscErrorCode VecCUDAResetArray(Vec vin)
2550: {
2552: #if defined(PETSC_HAVE_CUDA)
2553:   VecCUDACopyToGPU(vin);
2554:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar *) ((Vec_Seq*)vin->data)->unplacedarray;
2555:   ((Vec_Seq*)vin->data)->unplacedarray = 0;
2556:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2557: #endif
2558:   PetscObjectStateIncrease((PetscObject)vin);
2559:   return 0;
2560: }

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

2565:    This function has semantics similar to VecGetArray():  the pointer
2566:    returned by this function points to a consistent view of the vector
2567:    data.  This may involve a copy operation of data from the host to the
2568:    device if the data on the device is out of date.  If the device
2569:    memory hasn't been allocated previously it will be allocated as part
2570:    of this function call.  VecHIPGetArray() assumes that
2571:    the user will modify the vector data.  This is similar to
2572:    intent(inout) in fortran.

2574:    The HIP device pointer has to be released by calling
2575:    VecHIPRestoreArray().  Upon restoring the vector data
2576:    the data on the host will be marked as out of date.  A subsequent
2577:    access of the host data will thus incur a data transfer from the
2578:    device to the host.

2580:    Input Parameter:
2581: .  v - the vector

2583:    Output Parameter:
2584: .  a - the HIP device pointer

2586:    Fortran note:
2587:    This function is not currently available from Fortran.

2589:    Level: intermediate

2591: .seealso: VecHIPRestoreArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2592: @*/
2593: PETSC_EXTERN PetscErrorCode VecHIPGetArray(Vec v, PetscScalar **a)
2594: {
2596: #if defined(PETSC_HAVE_HIP)
2597:   *a   = 0;
2598:   VecHIPCopyToGPU(v);
2599:   *a   = ((Vec_HIP*)v->spptr)->GPUarray;
2600: #endif
2601:   return 0;
2602: }

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

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

2611:    Input Parameters:
2612: +  v - the vector
2613: -  a - the HIP device pointer.  This pointer is invalid after
2614:        VecHIPRestoreArray() returns.

2616:    Fortran note:
2617:    This function is not currently available from Fortran.

2619:    Level: intermediate

2621: .seealso: VecHIPGetArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2622: @*/
2623: PETSC_EXTERN PetscErrorCode VecHIPRestoreArray(Vec v, PetscScalar **a)
2624: {
2626: #if defined(PETSC_HAVE_HIP)
2627:   v->offloadmask = PETSC_OFFLOAD_GPU;
2628: #endif

2630:   PetscObjectStateIncrease((PetscObject)v);
2631:   return 0;
2632: }

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

2637:    This function is analogous to VecGetArrayRead():  The pointer
2638:    returned by this function points to a consistent view of the vector
2639:    data.  This may involve a copy operation of data from the host to the
2640:    device if the data on the device is out of date.  If the device
2641:    memory hasn't been allocated previously it will be allocated as part
2642:    of this function call.  VecHIPGetArrayRead() assumes that the
2643:    user will not modify the vector data.  This is analgogous to
2644:    intent(in) in Fortran.

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

2652:    Input Parameter:
2653: .  v - the vector

2655:    Output Parameter:
2656: .  a - the HIP pointer.

2658:    Fortran note:
2659:    This function is not currently available from Fortran.

2661:    Level: intermediate

2663: .seealso: VecHIPRestoreArrayRead(), VecHIPGetArray(), VecHIPGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2664: @*/
2665: PETSC_EXTERN PetscErrorCode VecHIPGetArrayRead(Vec v, const PetscScalar **a)
2666: {
2668: #if defined(PETSC_HAVE_HIP)
2669:   *a   = 0;
2670:   VecHIPCopyToGPU(v);
2671:   *a   = ((Vec_HIP*)v->spptr)->GPUarray;
2672: #endif
2673:   return 0;
2674: }

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

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

2684:    Input Parameters:
2685: +  v - the vector
2686: -  a - the HIP device pointer.  This pointer is invalid after
2687:        VecHIPRestoreArrayRead() returns.

2689:    Fortran note:
2690:    This function is not currently available from Fortran.

2692:    Level: intermediate

2694: .seealso: VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecHIPGetArray(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2695: @*/
2696: PETSC_EXTERN PetscErrorCode VecHIPRestoreArrayRead(Vec v, const PetscScalar **a)
2697: {
2699:   *a = NULL;
2700:   return 0;
2701: }

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

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

2712:    The device pointer needs to be released with
2713:    VecHIPRestoreArrayWrite().  When the pointer is released the
2714:    host data of the vector is marked as out of data.  Subsequent access
2715:    of the host data with e.g. VecGetArray() incurs a device to host data
2716:    transfer.

2718:    Input Parameter:
2719: .  v - the vector

2721:    Output Parameter:
2722: .  a - the HIP pointer

2724:    Fortran note:
2725:    This function is not currently available from Fortran.

2727:    Level: advanced

2729: .seealso: VecHIPRestoreArrayWrite(), VecHIPGetArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2730: @*/
2731: PETSC_EXTERN PetscErrorCode VecHIPGetArrayWrite(Vec v, PetscScalar **a)
2732: {
2734: #if defined(PETSC_HAVE_HIP)
2735:   *a   = 0;
2736:   VecHIPAllocateCheck(v);
2737:   *a   = ((Vec_HIP*)v->spptr)->GPUarray;
2738: #endif
2739:   return 0;
2740: }

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

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

2749:    Input Parameters:
2750: +  v - the vector
2751: -  a - the HIP device pointer.  This pointer is invalid after
2752:        VecHIPRestoreArrayWrite() returns.

2754:    Fortran note:
2755:    This function is not currently available from Fortran.

2757:    Level: intermediate

2759: .seealso: VecHIPGetArrayWrite(), VecHIPGetArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2760: @*/
2761: PETSC_EXTERN PetscErrorCode VecHIPRestoreArrayWrite(Vec v, PetscScalar **a)
2762: {
2764: #if defined(PETSC_HAVE_HIP)
2765:   v->offloadmask = PETSC_OFFLOAD_GPU;
2766: #endif

2768:   PetscObjectStateIncrease((PetscObject)v);
2769:   return 0;
2770: }

2772: /*@C
2773:    VecHIPPlaceArray - Allows one to replace the GPU array in a vector with a
2774:    GPU array provided by the user. This is useful to avoid copying an
2775:    array into a vector.

2777:    Not Collective

2779:    Input Parameters:
2780: +  vec - the vector
2781: -  array - the GPU array

2783:    Notes:
2784:    You can return to the original GPU array with a call to VecHIPResetArray()
2785:    It is not possible to use VecHIPPlaceArray() and VecPlaceArray() at the
2786:    same time on the same vector.

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

2793:    Level: developer

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

2797: @*/
2798: PetscErrorCode VecHIPPlaceArray(Vec vin,const PetscScalar a[])
2799: {
2801: #if defined(PETSC_HAVE_HIP)
2802:   VecHIPCopyToGPU(vin);
2804:   ((Vec_Seq*)vin->data)->unplacedarray  = (PetscScalar *) ((Vec_HIP*)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2805:   ((Vec_HIP*)vin->spptr)->GPUarray = (PetscScalar*)a;
2806:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2807: #endif
2808:   PetscObjectStateIncrease((PetscObject)vin);
2809:   return 0;
2810: }

2812: /*@C
2813:    VecHIPReplaceArray - Allows one to replace the GPU array in a vector
2814:    with a GPU array provided by the user. This is useful to avoid copying
2815:    a GPU array into a vector.

2817:    Not Collective

2819:    Input Parameters:
2820: +  vec - the vector
2821: -  array - the GPU array

2823:    Notes:
2824:    This permanently replaces the GPU array and frees the memory associated
2825:    with the old GPU array.

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

2830:    Not supported from Fortran

2832:    Level: developer

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

2836: @*/
2837: PetscErrorCode VecHIPReplaceArray(Vec vin,const PetscScalar a[])
2838: {
2840: #if defined(PETSC_HAVE_HIP)
2841:   hipFree(((Vec_HIP*)vin->spptr)->GPUarray);
2842:   ((Vec_HIP*)vin->spptr)->GPUarray = (PetscScalar*)a;
2843:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2844: #endif
2845:   PetscObjectStateIncrease((PetscObject)vin);
2846:   return 0;
2847: }

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

2853:    Not Collective

2855:    Input Parameters:
2856: .  vec - the vector

2858:    Level: developer

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

2862: @*/
2863: PetscErrorCode VecHIPResetArray(Vec vin)
2864: {
2866: #if defined(PETSC_HAVE_HIP)
2867:   VecHIPCopyToGPU(vin);
2868:   ((Vec_HIP*)vin->spptr)->GPUarray = (PetscScalar *) ((Vec_Seq*)vin->data)->unplacedarray;
2869:   ((Vec_Seq*)vin->data)->unplacedarray = 0;
2870:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2871: #endif
2872:   PetscObjectStateIncrease((PetscObject)vin);
2873:   return 0;
2874: }

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

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

2883:     Collective on Vec

2885:     Input Parameters:
2886: +   x - a vector to mimic
2887: -   n - the number of vectors to obtain

2889:     Output Parameters:
2890: +   y - Fortran90 pointer to the array of vectors
2891: -   ierr - error code

2893:     Example of Usage:
2894: .vb
2895: #include <petsc/finclude/petscvec.h>
2896:     use petscvec

2898:     Vec x
2899:     Vec, pointer :: y(:)
2900:     ....
2901:     call VecDuplicateVecsF90(x,2,y,ierr)
2902:     call VecSet(y(2),alpha,ierr)
2903:     call VecSet(y(2),alpha,ierr)
2904:     ....
2905:     call VecDestroyVecsF90(2,y,ierr)
2906: .ve

2908:     Notes:
2909:     Not yet supported for all F90 compilers

2911:     Use VecDestroyVecsF90() to free the space.

2913:     Level: beginner

2915: .seealso:  VecDestroyVecsF90(), VecDuplicateVecs()

2917: M*/

2919: /*MC
2920:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2921:     VecGetArrayF90().

2923:     Synopsis:
2924:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2926:     Logically Collective on Vec

2928:     Input Parameters:
2929: +   x - vector
2930: -   xx_v - the Fortran90 pointer to the array

2932:     Output Parameter:
2933: .   ierr - error code

2935:     Example of Usage:
2936: .vb
2937: #include <petsc/finclude/petscvec.h>
2938:     use petscvec

2940:     PetscScalar, pointer :: xx_v(:)
2941:     ....
2942:     call VecGetArrayF90(x,xx_v,ierr)
2943:     xx_v(3) = a
2944:     call VecRestoreArrayF90(x,xx_v,ierr)
2945: .ve

2947:     Level: beginner

2949: .seealso:  VecGetArrayF90(), VecGetArray(), VecRestoreArray(), VecRestoreArrayReadF90()

2951: M*/

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

2956:     Synopsis:
2957:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

2959:     Collective on Vec

2961:     Input Parameters:
2962: +   n - the number of vectors previously obtained
2963: -   x - pointer to array of vector pointers

2965:     Output Parameter:
2966: .   ierr - error code

2968:     Notes:
2969:     Not yet supported for all F90 compilers

2971:     Level: beginner

2973: .seealso:  VecDestroyVecs(), VecDuplicateVecsF90()

2975: M*/

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

2983:     Synopsis:
2984:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2986:     Logically Collective on Vec

2988:     Input Parameter:
2989: .   x - vector

2991:     Output Parameters:
2992: +   xx_v - the Fortran90 pointer to the array
2993: -   ierr - error code

2995:     Example of Usage:
2996: .vb
2997: #include <petsc/finclude/petscvec.h>
2998:     use petscvec

3000:     PetscScalar, pointer :: xx_v(:)
3001:     ....
3002:     call VecGetArrayF90(x,xx_v,ierr)
3003:     xx_v(3) = a
3004:     call VecRestoreArrayF90(x,xx_v,ierr)
3005: .ve

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

3009:     Level: beginner

3011: .seealso:  VecRestoreArrayF90(), VecGetArray(), VecRestoreArray(), VecGetArrayReadF90()

3013: M*/

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

3021:     Synopsis:
3022:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

3024:     Logically Collective on Vec

3026:     Input Parameter:
3027: .   x - vector

3029:     Output Parameters:
3030: +   xx_v - the Fortran90 pointer to the array
3031: -   ierr - error code

3033:     Example of Usage:
3034: .vb
3035: #include <petsc/finclude/petscvec.h>
3036:     use petscvec

3038:     PetscScalar, pointer :: xx_v(:)
3039:     ....
3040:     call VecGetArrayReadF90(x,xx_v,ierr)
3041:     a = xx_v(3)
3042:     call VecRestoreArrayReadF90(x,xx_v,ierr)
3043: .ve

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

3047:     Level: beginner

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

3051: M*/

3053: /*MC
3054:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
3055:     VecGetArrayReadF90().

3057:     Synopsis:
3058:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

3060:     Logically Collective on Vec

3062:     Input Parameters:
3063: +   x - vector
3064: -   xx_v - the Fortran90 pointer to the array

3066:     Output Parameter:
3067: .   ierr - error code

3069:     Example of Usage:
3070: .vb
3071: #include <petsc/finclude/petscvec.h>
3072:     use petscvec

3074:     PetscScalar, pointer :: xx_v(:)
3075:     ....
3076:     call VecGetArrayReadF90(x,xx_v,ierr)
3077:     a = xx_v(3)
3078:     call VecRestoreArrayReadF90(x,xx_v,ierr)
3079: .ve

3081:     Level: beginner

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

3085: M*/

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

3092:    Logically Collective

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

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

3104:    Level: developer

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

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

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

3126:   VecGetLocalSize(x,&N);
3128:   VecGetArray(x,&aa);

3130:   PetscMalloc1(m,a);
3131:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
3132:   *a -= mstart;
3133:   return 0;
3134: }

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

3141:    Logically Collective

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

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

3153:    Level: developer

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

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

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

3175:   VecGetLocalSize(x,&N);
3177:   VecGetArrayWrite(x,&aa);

3179:   PetscMalloc1(m,a);
3180:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
3181:   *a -= mstart;
3182:   return 0;
3183: }

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

3188:    Logically Collective

3190:    Input Parameters:
3191: +  x - the vector
3192: .  m - first dimension of two dimensional array
3193: .  n - second dimension of the two dimensional array
3194: .  mstart - first index you will use in first coordinate direction (often 0)
3195: .  nstart - first index in the second coordinate direction (often 0)
3196: -  a - location of pointer to array obtained from VecGetArray2d()

3198:    Level: developer

3200:    Notes:
3201:    For regular PETSc vectors this routine does not involve any copies. For
3202:    any special vectors that do not store local vector data in a contiguous
3203:    array, this routine will copy the data back into the underlying
3204:    vector data structure from the array obtained with VecGetArray().

3206:    This routine actually zeros out the a pointer.

3208: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3209:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3210:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3211: @*/
3212: PetscErrorCode  VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3213: {
3214:   void           *dummy;

3219:   dummy = (void*)(*a + mstart);
3220:   PetscFree(dummy);
3221:   VecRestoreArray(x,NULL);
3222:   return 0;
3223: }

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

3228:    Logically Collective

3230:    Input Parameters:
3231: +  x - the vector
3232: .  m - first dimension of two dimensional array
3233: .  n - second dimension of the two dimensional array
3234: .  mstart - first index you will use in first coordinate direction (often 0)
3235: .  nstart - first index in the second coordinate direction (often 0)
3236: -  a - location of pointer to array obtained from VecGetArray2d()

3238:    Level: developer

3240:    Notes:
3241:    For regular PETSc vectors this routine does not involve any copies. For
3242:    any special vectors that do not store local vector data in a contiguous
3243:    array, this routine will copy the data back into the underlying
3244:    vector data structure from the array obtained with VecGetArray().

3246:    This routine actually zeros out the a pointer.

3248: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3249:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3250:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3251: @*/
3252: PetscErrorCode  VecRestoreArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3253: {
3254:   void           *dummy;

3259:   dummy = (void*)(*a + mstart);
3260:   PetscFree(dummy);
3261:   VecRestoreArrayWrite(x,NULL);
3262:   return 0;
3263: }

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

3270:    Logically Collective

3272:    Input Parameters:
3273: +  x - the vector
3274: .  m - first dimension of two dimensional array
3275: -  mstart - first index you will use in first coordinate direction (often 0)

3277:    Output Parameter:
3278: .  a - location to put pointer to the array

3280:    Level: developer

3282:   Notes:
3283:    For a vector obtained from DMCreateLocalVector() mstart are likely
3284:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3285:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3289: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3290:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3291:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3292: @*/
3293: PetscErrorCode  VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3294: {
3295:   PetscInt       N;

3300:   VecGetLocalSize(x,&N);
3302:   VecGetArray(x,a);
3303:   *a  -= mstart;
3304:   return 0;
3305: }

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

3312:    Logically Collective

3314:    Input Parameters:
3315: +  x - the vector
3316: .  m - first dimension of two dimensional array
3317: -  mstart - first index you will use in first coordinate direction (often 0)

3319:    Output Parameter:
3320: .  a - location to put pointer to the array

3322:    Level: developer

3324:   Notes:
3325:    For a vector obtained from DMCreateLocalVector() mstart are likely
3326:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3327:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3331: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3332:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3333:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3334: @*/
3335: PetscErrorCode  VecGetArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3336: {
3337:   PetscInt       N;

3342:   VecGetLocalSize(x,&N);
3344:   VecGetArrayWrite(x,a);
3345:   *a  -= mstart;
3346:   return 0;
3347: }

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

3352:    Logically Collective

3354:    Input Parameters:
3355: +  x - the vector
3356: .  m - first dimension of two dimensional array
3357: .  mstart - first index you will use in first coordinate direction (often 0)
3358: -  a - location of pointer to array obtained from VecGetArray21()

3360:    Level: developer

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

3368:    This routine actually zeros out the a pointer.

3370: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3371:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3372:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3373: @*/
3374: PetscErrorCode  VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3375: {
3378:   VecRestoreArray(x,NULL);
3379:   return 0;
3380: }

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

3385:    Logically Collective

3387:    Input Parameters:
3388: +  x - the vector
3389: .  m - first dimension of two dimensional array
3390: .  mstart - first index you will use in first coordinate direction (often 0)
3391: -  a - location of pointer to array obtained from VecGetArray21()

3393:    Level: developer

3395:    Notes:
3396:    For regular PETSc vectors this routine does not involve any copies. For
3397:    any special vectors that do not store local vector data in a contiguous
3398:    array, this routine will copy the data back into the underlying
3399:    vector data structure from the array obtained with VecGetArray1d().

3401:    This routine actually zeros out the a pointer.

3403: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3404:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3405:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3406: @*/
3407: PetscErrorCode  VecRestoreArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3408: {
3411:   VecRestoreArrayWrite(x,NULL);
3412:   return 0;
3413: }

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

3420:    Logically Collective

3422:    Input Parameters:
3423: +  x - the vector
3424: .  m - first dimension of three dimensional array
3425: .  n - second dimension of three dimensional array
3426: .  p - third dimension of three dimensional array
3427: .  mstart - first index you will use in first coordinate direction (often 0)
3428: .  nstart - first index in the second coordinate direction (often 0)
3429: -  pstart - first index in the third coordinate direction (often 0)

3431:    Output Parameter:
3432: .  a - location to put pointer to the array

3434:    Level: developer

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

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

3444: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3445:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3446:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3447: @*/
3448: PetscErrorCode  VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3449: {
3450:   PetscInt       i,N,j;
3451:   PetscScalar    *aa,**b;

3456:   VecGetLocalSize(x,&N);
3458:   VecGetArray(x,&aa);

3460:   PetscMalloc(m*sizeof(PetscScalar**)+m*n*sizeof(PetscScalar*),a);
3461:   b    = (PetscScalar**)((*a) + m);
3462:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3463:   for (i=0; i<m; i++)
3464:     for (j=0; j<n; j++)
3465:       b[i*n+j] = aa + i*n*p + j*p - pstart;
3466:   *a -= mstart;
3467:   return 0;
3468: }

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

3475:    Logically Collective

3477:    Input Parameters:
3478: +  x - the vector
3479: .  m - first dimension of three dimensional array
3480: .  n - second dimension of three dimensional array
3481: .  p - third dimension of three dimensional array
3482: .  mstart - first index you will use in first coordinate direction (often 0)
3483: .  nstart - first index in the second coordinate direction (often 0)
3484: -  pstart - first index in the third coordinate direction (often 0)

3486:    Output Parameter:
3487: .  a - location to put pointer to the array

3489:    Level: developer

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

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

3499: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3500:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3501:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3502: @*/
3503: PetscErrorCode  VecGetArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3504: {
3505:   PetscInt       i,N,j;
3506:   PetscScalar    *aa,**b;

3511:   VecGetLocalSize(x,&N);
3513:   VecGetArrayWrite(x,&aa);

3515:   PetscMalloc(m*sizeof(PetscScalar**)+m*n*sizeof(PetscScalar*),a);
3516:   b    = (PetscScalar**)((*a) + m);
3517:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3518:   for (i=0; i<m; i++)
3519:     for (j=0; j<n; j++)
3520:       b[i*n+j] = aa + i*n*p + j*p - pstart;

3522:   *a -= mstart;
3523:   return 0;
3524: }

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

3529:    Logically Collective

3531:    Input Parameters:
3532: +  x - the vector
3533: .  m - first dimension of three dimensional array
3534: .  n - second dimension of the three dimensional array
3535: .  p - third dimension of the three dimensional array
3536: .  mstart - first index you will use in first coordinate direction (often 0)
3537: .  nstart - first index in the second coordinate direction (often 0)
3538: .  pstart - first index in the third coordinate direction (often 0)
3539: -  a - location of pointer to array obtained from VecGetArray3d()

3541:    Level: developer

3543:    Notes:
3544:    For regular PETSc vectors this routine does not involve any copies. For
3545:    any special vectors that do not store local vector data in a contiguous
3546:    array, this routine will copy the data back into the underlying
3547:    vector data structure from the array obtained with VecGetArray().

3549:    This routine actually zeros out the a pointer.

3551: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3552:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3553:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3554: @*/
3555: PetscErrorCode  VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3556: {
3557:   void           *dummy;

3562:   dummy = (void*)(*a + mstart);
3563:   PetscFree(dummy);
3564:   VecRestoreArray(x,NULL);
3565:   return 0;
3566: }

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

3571:    Logically Collective

3573:    Input Parameters:
3574: +  x - the vector
3575: .  m - first dimension of three dimensional array
3576: .  n - second dimension of the three dimensional array
3577: .  p - third dimension of the three dimensional array
3578: .  mstart - first index you will use in first coordinate direction (often 0)
3579: .  nstart - first index in the second coordinate direction (often 0)
3580: .  pstart - first index in the third coordinate direction (often 0)
3581: -  a - location of pointer to array obtained from VecGetArray3d()

3583:    Level: developer

3585:    Notes:
3586:    For regular PETSc vectors this routine does not involve any copies. For
3587:    any special vectors that do not store local vector data in a contiguous
3588:    array, this routine will copy the data back into the underlying
3589:    vector data structure from the array obtained with VecGetArray().

3591:    This routine actually zeros out the a pointer.

3593: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3594:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3595:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3596: @*/
3597: PetscErrorCode  VecRestoreArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3598: {
3599:   void           *dummy;

3604:   dummy = (void*)(*a + mstart);
3605:   PetscFree(dummy);
3606:   VecRestoreArrayWrite(x,NULL);
3607:   return 0;
3608: }

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

3615:    Logically Collective

3617:    Input Parameters:
3618: +  x - the vector
3619: .  m - first dimension of four dimensional array
3620: .  n - second dimension of four dimensional array
3621: .  p - third dimension of four dimensional array
3622: .  q - fourth dimension of four dimensional array
3623: .  mstart - first index you will use in first coordinate direction (often 0)
3624: .  nstart - first index in the second coordinate direction (often 0)
3625: .  pstart - first index in the third coordinate direction (often 0)
3626: -  qstart - first index in the fourth coordinate direction (often 0)

3628:    Output Parameter:
3629: .  a - location to put pointer to the array

3631:    Level: beginner

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

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

3641: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3642:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3643:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3644: @*/
3645: PetscErrorCode  VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3646: {
3647:   PetscInt       i,N,j,k;
3648:   PetscScalar    *aa,***b,**c;

3653:   VecGetLocalSize(x,&N);
3655:   VecGetArray(x,&aa);

3657:   PetscMalloc(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p*sizeof(PetscScalar*),a);
3658:   b    = (PetscScalar***)((*a) + m);
3659:   c    = (PetscScalar**)(b + m*n);
3660:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3661:   for (i=0; i<m; i++)
3662:     for (j=0; j<n; j++)
3663:       b[i*n+j] = c + i*n*p + j*p - pstart;
3664:   for (i=0; i<m; i++)
3665:     for (j=0; j<n; j++)
3666:       for (k=0; k<p; k++)
3667:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3668:   *a -= mstart;
3669:   return 0;
3670: }

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

3677:    Logically Collective

3679:    Input Parameters:
3680: +  x - the vector
3681: .  m - first dimension of four dimensional array
3682: .  n - second dimension of four dimensional array
3683: .  p - third dimension of four dimensional array
3684: .  q - fourth dimension of four dimensional array
3685: .  mstart - first index you will use in first coordinate direction (often 0)
3686: .  nstart - first index in the second coordinate direction (often 0)
3687: .  pstart - first index in the third coordinate direction (often 0)
3688: -  qstart - first index in the fourth coordinate direction (often 0)

3690:    Output Parameter:
3691: .  a - location to put pointer to the array

3693:    Level: beginner

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

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

3703: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3704:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3705:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3706: @*/
3707: PetscErrorCode  VecGetArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3708: {
3709:   PetscInt       i,N,j,k;
3710:   PetscScalar    *aa,***b,**c;

3715:   VecGetLocalSize(x,&N);
3717:   VecGetArrayWrite(x,&aa);

3719:   PetscMalloc(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p*sizeof(PetscScalar*),a);
3720:   b    = (PetscScalar***)((*a) + m);
3721:   c    = (PetscScalar**)(b + m*n);
3722:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3723:   for (i=0; i<m; i++)
3724:     for (j=0; j<n; j++)
3725:       b[i*n+j] = c + i*n*p + j*p - pstart;
3726:   for (i=0; i<m; i++)
3727:     for (j=0; j<n; j++)
3728:       for (k=0; k<p; k++)
3729:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3730:   *a -= mstart;
3731:   return 0;
3732: }

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

3737:    Logically Collective

3739:    Input Parameters:
3740: +  x - the vector
3741: .  m - first dimension of four dimensional array
3742: .  n - second dimension of the four dimensional array
3743: .  p - third dimension of the four dimensional array
3744: .  q - fourth dimension of the four dimensional array
3745: .  mstart - first index you will use in first coordinate direction (often 0)
3746: .  nstart - first index in the second coordinate direction (often 0)
3747: .  pstart - first index in the third coordinate direction (often 0)
3748: .  qstart - first index in the fourth coordinate direction (often 0)
3749: -  a - location of pointer to array obtained from VecGetArray4d()

3751:    Level: beginner

3753:    Notes:
3754:    For regular PETSc vectors this routine does not involve any copies. For
3755:    any special vectors that do not store local vector data in a contiguous
3756:    array, this routine will copy the data back into the underlying
3757:    vector data structure from the array obtained with VecGetArray().

3759:    This routine actually zeros out the a pointer.

3761: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3762:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3763:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3764: @*/
3765: PetscErrorCode  VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3766: {
3767:   void           *dummy;

3772:   dummy = (void*)(*a + mstart);
3773:   PetscFree(dummy);
3774:   VecRestoreArray(x,NULL);
3775:   return 0;
3776: }

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

3781:    Logically Collective

3783:    Input Parameters:
3784: +  x - the vector
3785: .  m - first dimension of four dimensional array
3786: .  n - second dimension of the four dimensional array
3787: .  p - third dimension of the four dimensional array
3788: .  q - fourth dimension of the four dimensional array
3789: .  mstart - first index you will use in first coordinate direction (often 0)
3790: .  nstart - first index in the second coordinate direction (often 0)
3791: .  pstart - first index in the third coordinate direction (often 0)
3792: .  qstart - first index in the fourth coordinate direction (often 0)
3793: -  a - location of pointer to array obtained from VecGetArray4d()

3795:    Level: beginner

3797:    Notes:
3798:    For regular PETSc vectors this routine does not involve any copies. For
3799:    any special vectors that do not store local vector data in a contiguous
3800:    array, this routine will copy the data back into the underlying
3801:    vector data structure from the array obtained with VecGetArray().

3803:    This routine actually zeros out the a pointer.

3805: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3806:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3807:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3808: @*/
3809: PetscErrorCode  VecRestoreArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3810: {
3811:   void           *dummy;

3816:   dummy = (void*)(*a + mstart);
3817:   PetscFree(dummy);
3818:   VecRestoreArrayWrite(x,NULL);
3819:   return 0;
3820: }

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

3827:    Logically Collective

3829:    Input Parameters:
3830: +  x - the vector
3831: .  m - first dimension of two dimensional array
3832: .  n - second dimension of two dimensional array
3833: .  mstart - first index you will use in first coordinate direction (often 0)
3834: -  nstart - first index in the second coordinate direction (often 0)

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

3839:    Level: developer

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

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

3849: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3850:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3851:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3852: @*/
3853: PetscErrorCode  VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3854: {
3855:   PetscInt          i,N;
3856:   const PetscScalar *aa;

3861:   VecGetLocalSize(x,&N);
3863:   VecGetArrayRead(x,&aa);

3865:   PetscMalloc1(m,a);
3866:   for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
3867:   *a -= mstart;
3868:   return 0;
3869: }

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

3874:    Logically Collective

3876:    Input Parameters:
3877: +  x - the vector
3878: .  m - first dimension of two dimensional array
3879: .  n - second dimension of the two dimensional array
3880: .  mstart - first index you will use in first coordinate direction (often 0)
3881: .  nstart - first index in the second coordinate direction (often 0)
3882: -  a - location of pointer to array obtained from VecGetArray2d()

3884:    Level: developer

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

3892:    This routine actually zeros out the a pointer.

3894: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3895:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3896:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3897: @*/
3898: PetscErrorCode  VecRestoreArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3899: {
3900:   void           *dummy;

3905:   dummy = (void*)(*a + mstart);
3906:   PetscFree(dummy);
3907:   VecRestoreArrayRead(x,NULL);
3908:   return 0;
3909: }

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

3916:    Logically Collective

3918:    Input Parameters:
3919: +  x - the vector
3920: .  m - first dimension of two dimensional array
3921: -  mstart - first index you will use in first 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 are likely
3930:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3931:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3935: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3936:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3937:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3938: @*/
3939: PetscErrorCode  VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3940: {
3941:   PetscInt       N;

3946:   VecGetLocalSize(x,&N);
3948:   VecGetArrayRead(x,(const PetscScalar**)a);
3949:   *a  -= mstart;
3950:   return 0;
3951: }

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

3956:    Logically Collective

3958:    Input Parameters:
3959: +  x - the vector
3960: .  m - first dimension of two dimensional array
3961: .  mstart - first index you will use in first coordinate direction (often 0)
3962: -  a - location of pointer to array obtained from VecGetArray21()

3964:    Level: developer

3966:    Notes:
3967:    For regular PETSc vectors this routine does not involve any copies. For
3968:    any special vectors that do not store local vector data in a contiguous
3969:    array, this routine will copy the data back into the underlying
3970:    vector data structure from the array obtained with VecGetArray1dRead().

3972:    This routine actually zeros out the a pointer.

3974: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3975:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3976:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3977: @*/
3978: PetscErrorCode  VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3979: {
3982:   VecRestoreArrayRead(x,NULL);
3983:   return 0;
3984: }

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

3991:    Logically Collective

3993:    Input Parameters:
3994: +  x - the vector
3995: .  m - first dimension of three dimensional array
3996: .  n - second dimension of three dimensional array
3997: .  p - third dimension of three dimensional array
3998: .  mstart - first index you will use in first coordinate direction (often 0)
3999: .  nstart - first index in the second coordinate direction (often 0)
4000: -  pstart - first index in the third coordinate direction (often 0)

4002:    Output Parameter:
4003: .  a - location to put pointer to the array

4005:    Level: developer

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

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

4015: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
4016:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
4017:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
4018: @*/
4019: PetscErrorCode  VecGetArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
4020: {
4021:   PetscInt          i,N,j;
4022:   const PetscScalar *aa;
4023:   PetscScalar       **b;

4028:   VecGetLocalSize(x,&N);
4030:   VecGetArrayRead(x,&aa);

4032:   PetscMalloc(m*sizeof(PetscScalar**)+m*n*sizeof(PetscScalar*),a);
4033:   b    = (PetscScalar**)((*a) + m);
4034:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
4035:   for (i=0; i<m; i++)
4036:     for (j=0; j<n; j++)
4037:       b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;
4038:   *a -= mstart;
4039:   return 0;
4040: }

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

4045:    Logically Collective

4047:    Input Parameters:
4048: +  x - the vector
4049: .  m - first dimension of three dimensional array
4050: .  n - second dimension of the three dimensional array
4051: .  p - third dimension of the three dimensional array
4052: .  mstart - first index you will use in first coordinate direction (often 0)
4053: .  nstart - first index in the second coordinate direction (often 0)
4054: .  pstart - first index in the third coordinate direction (often 0)
4055: -  a - location of pointer to array obtained from VecGetArray3dRead()

4057:    Level: developer

4059:    Notes:
4060:    For regular PETSc vectors this routine does not involve any copies. For
4061:    any special vectors that do not store local vector data in a contiguous
4062:    array, this routine will copy the data back into the underlying
4063:    vector data structure from the array obtained with VecGetArray().

4065:    This routine actually zeros out the a pointer.

4067: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
4068:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
4069:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
4070: @*/
4071: PetscErrorCode  VecRestoreArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
4072: {
4073:   void           *dummy;

4078:   dummy = (void*)(*a + mstart);
4079:   PetscFree(dummy);
4080:   VecRestoreArrayRead(x,NULL);
4081:   return 0;
4082: }

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

4089:    Logically Collective

4091:    Input Parameters:
4092: +  x - the vector
4093: .  m - first dimension of four dimensional array
4094: .  n - second dimension of four dimensional array
4095: .  p - third dimension of four dimensional array
4096: .  q - fourth dimension of four dimensional array
4097: .  mstart - first index you will use in first coordinate direction (often 0)
4098: .  nstart - first index in the second coordinate direction (often 0)
4099: .  pstart - first index in the third coordinate direction (often 0)
4100: -  qstart - first index in the fourth coordinate direction (often 0)

4102:    Output Parameter:
4103: .  a - location to put pointer to the array

4105:    Level: beginner

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

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

4115: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
4116:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
4117:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
4118: @*/
4119: PetscErrorCode  VecGetArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
4120: {
4121:   PetscInt          i,N,j,k;
4122:   const PetscScalar *aa;
4123:   PetscScalar       ***b,**c;

4128:   VecGetLocalSize(x,&N);
4130:   VecGetArrayRead(x,&aa);

4132:   PetscMalloc(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p*sizeof(PetscScalar*),a);
4133:   b    = (PetscScalar***)((*a) + m);
4134:   c    = (PetscScalar**)(b + m*n);
4135:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
4136:   for (i=0; i<m; i++)
4137:     for (j=0; j<n; j++)
4138:       b[i*n+j] = c + i*n*p + j*p - pstart;
4139:   for (i=0; i<m; i++)
4140:     for (j=0; j<n; j++)
4141:       for (k=0; k<p; k++)
4142:         c[i*n*p+j*p+k] = (PetscScalar*) aa + i*n*p*q + j*p*q + k*q - qstart;
4143:   *a -= mstart;
4144:   return 0;
4145: }

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

4150:    Logically Collective

4152:    Input Parameters:
4153: +  x - the vector
4154: .  m - first dimension of four dimensional array
4155: .  n - second dimension of the four dimensional array
4156: .  p - third dimension of the four dimensional array
4157: .  q - fourth dimension of the four dimensional array
4158: .  mstart - first index you will use in first coordinate direction (often 0)
4159: .  nstart - first index in the second coordinate direction (often 0)
4160: .  pstart - first index in the third coordinate direction (often 0)
4161: .  qstart - first index in the fourth coordinate direction (often 0)
4162: -  a - location of pointer to array obtained from VecGetArray4dRead()

4164:    Level: beginner

4166:    Notes:
4167:    For regular PETSc vectors this routine does not involve any copies. For
4168:    any special vectors that do not store local vector data in a contiguous
4169:    array, this routine will copy the data back into the underlying
4170:    vector data structure from the array obtained with VecGetArray().

4172:    This routine actually zeros out the a pointer.

4174: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
4175:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
4176:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
4177: @*/
4178: PetscErrorCode  VecRestoreArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
4179: {
4180:   void           *dummy;

4185:   dummy = (void*)(*a + mstart);
4186:   PetscFree(dummy);
4187:   VecRestoreArrayRead(x,NULL);
4188:   return 0;
4189: }

4191: #if defined(PETSC_USE_DEBUG)

4193: /*@
4194:    VecLockGet  - Gets the current lock status of a vector

4196:    Logically Collective on Vec

4198:    Input Parameter:
4199: .  x - the vector

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

4205:    Level: beginner

4207: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop()
4208: @*/
4209: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
4210: {
4212:   *state = x->lock;
4213:   return 0;
4214: }

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

4219:    Logically Collective on Vec

4221:    Input Parameter:
4222: .  x - the vector

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

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

4230:    Level: beginner

4232: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPop(), VecLockGet()
4233: @*/
4234: PetscErrorCode VecLockReadPush(Vec x)
4235: {
4238:   x->lock++;
4239:   return 0;
4240: }

4242: /*@
4243:    VecLockReadPop  - Pops a read-only lock from a vector

4245:    Logically Collective on Vec

4247:    Input Parameter:
4248: .  x - the vector

4250:    Level: beginner

4252: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockGet()
4253: @*/
4254: PetscErrorCode VecLockReadPop(Vec x)
4255: {
4257:   x->lock--;
4259:   return 0;
4260: }

4262: /*@C
4263:    VecLockWriteSet_Private  - Lock or unlock a vector for exclusive read/write access

4265:    Logically Collective on Vec

4267:    Input Parameters:
4268: +  x   - the vector
4269: -  flg - PETSC_TRUE to lock the vector for writing; PETSC_FALSE to unlock it.

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

4277:        VecGetArray(x,&xdata); // begin phase
4278:        VecLockWriteSet_Private(v,PETSC_TRUE);

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

4282:        VecRestoreArray(x,&vdata); // end phase
4283:        VecLockWriteSet_Private(v,PETSC_FALSE);

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

4288:    Level: beginner

4290: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop(), VecLockGet()
4291: @*/
4292: PetscErrorCode VecLockWriteSet_Private(Vec x,PetscBool flg)
4293: {
4295:   if (flg) {
4298:     else x->lock = -1;
4299:   } else {
4301:     x->lock = 0;
4302:   }
4303:   return 0;
4304: }

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

4309:    Level: deprecated

4311: .seealso: VecLockReadPush()
4312: @*/
4313: PetscErrorCode VecLockPush(Vec x)
4314: {
4315:   VecLockReadPush(x);
4316:   return 0;
4317: }

4319: /*@
4320:    VecLockPop  - Pops a read-only lock from a vector

4322:    Level: deprecated

4324: .seealso: VecLockReadPop()
4325: @*/
4326: PetscErrorCode VecLockPop(Vec x)
4327: {
4328:   VecLockReadPop(x);
4329:   return 0;
4330: }

4332: #endif