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: static PetscInt VecGetSubVectorSavedStateId = -1;

 18: PETSC_EXTERN PetscErrorCode VecValidValues(Vec vec,PetscInt argnum,PetscBool begin)
 19: {
 20: #if defined(PETSC_USE_DEBUG)
 21:   PetscErrorCode    ierr;
 22:   PetscInt          n,i;
 23:   const PetscScalar *x;

 26: #if defined(PETSC_HAVE_DEVICE)
 27:   if ((vec->petscnative || vec->ops->getarray) && (vec->offloadmask & PETSC_OFFLOAD_CPU)) {
 28: #else
 29:   if (vec->petscnative || vec->ops->getarray) {
 30: #endif
 31:     VecGetLocalSize(vec,&n);
 32:     VecGetArrayRead(vec,&x);
 33:     for (i=0; i<n; i++) {
 34:       if (begin) {
 35:         if (PetscIsInfOrNanScalar(x[i])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FP,"Vec entry at local location %D is not-a-number or infinite at beginning of function: Parameter number %D",i,argnum);
 36:       } else {
 37:         if (PetscIsInfOrNanScalar(x[i])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FP,"Vec entry at local location %D is not-a-number or infinite at end of function: Parameter number %D",i,argnum);
 38:       }
 39:     }
 40:     VecRestoreArrayRead(vec,&x);
 41:   }
 42: #else
 44: #endif
 45:   return(0);
 46: }

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

 51:    Logically Collective on Vec

 53:    Input Parameters:
 54: .  x, y  - the vectors

 56:    Output Parameter:
 57: .  max - the result

 59:    Level: advanced

 61:    Notes:
 62:     x and y may be the same vector
 63:           if a particular y_i is zero, it is treated as 1 in the above formula

 65: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs()
 66: @*/
 67: PetscErrorCode  VecMaxPointwiseDivide(Vec x,Vec y,PetscReal *max)
 68: {

 78:   VecCheckSameSize(x,1,y,2);
 79:   (*x->ops->maxpointwisedivide)(x,y,max);
 80:   return(0);
 81: }

 83: /*@
 84:    VecDot - Computes the vector dot product.

 86:    Collective on Vec

 88:    Input Parameters:
 89: .  x, y - the vectors

 91:    Output Parameter:
 92: .  val - the dot product

 94:    Performance Issues:
 95: $    per-processor memory bandwidth
 96: $    interprocessor latency
 97: $    work load inbalance that causes certain processes to arrive much earlier than others

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

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

110:    Level: intermediate

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

125:   VecCheckSameSize(x,1,y,2);

127:   PetscLogEventBegin(VEC_Dot,x,y,0,0);
128:   (*x->ops->dot)(x,y,val);
129:   PetscLogEventEnd(VEC_Dot,x,y,0,0);
130:   return(0);
131: }

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

136:    Collective on Vec

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

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

144:    Performance Issues:
145: $    per-processor memory bandwidth
146: $    interprocessor latency
147: $    work load inbalance that causes certain processes to arrive much earlier than others

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

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

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

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

159:    Level: intermediate

161: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDot(), VecDotNorm2()
162: @*/
163: PetscErrorCode  VecDotRealPart(Vec x,Vec y,PetscReal *val)
164: {
166:   PetscScalar    fdot;

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

174: /*@
175:    VecNorm  - Computes the vector norm.

177:    Collective on Vec

179:    Input Parameters:
180: +  x - the vector
181: -  type - one of NORM_1, NORM_2, NORM_INFINITY.  Also available
182:           NORM_1_AND_2, which computes both norms and stores them
183:           in a two element array.

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

188:    Notes:
189: $     NORM_1 denotes sum_i |x_i|
190: $     NORM_2 denotes sqrt(sum_i |x_i|^2)
191: $     NORM_INFINITY denotes max_i |x_i|

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

198:    Level: intermediate

200:    Performance Issues:
201: $    per-processor memory bandwidth
202: $    interprocessor latency
203: $    work load inbalance that causes certain processes to arrive much earlier than others

205: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNormAvailable(),
206:           VecNormBegin(), VecNormEnd()

208: @*/

210: PetscErrorCode  VecNorm(Vec x,NormType type,PetscReal *val)
211: {
212:   PetscBool      flg;


220:   /*
221:    * Cached data?
222:    */
223:   if (type!=NORM_1_AND_2) {
224:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,flg);
225:     if (flg) return(0);
226:   }
227:   PetscLogEventBegin(VEC_Norm,x,0,0,0);
228:   (*x->ops->norm)(x,type,val);
229:   PetscLogEventEnd(VEC_Norm,x,0,0,0);
230:   if (type!=NORM_1_AND_2) {
231:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[type],*val);
232:   }
233:   return(0);
234: }

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

239:    Not Collective

241:    Input Parameters:
242: +  x - the vector
243: -  type - one of NORM_1, NORM_2, NORM_INFINITY.  Also available
244:           NORM_1_AND_2, which computes both norms and stores them
245:           in a two element array.

247:    Output Parameter:
248: +  available - PETSC_TRUE if the val returned is valid
249: -  val - the norm

251:    Notes:
252: $     NORM_1 denotes sum_i |x_i|
253: $     NORM_2 denotes sqrt(sum_i (x_i)^2)
254: $     NORM_INFINITY denotes max_i |x_i|

256:    Level: intermediate

258:    Performance Issues:
259: $    per-processor memory bandwidth
260: $    interprocessor latency
261: $    work load inbalance that causes certain processes to arrive much earlier than others

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

268: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNorm()
269:           VecNormBegin(), VecNormEnd()

271: @*/
272: PetscErrorCode  VecNormAvailable(Vec x,NormType type,PetscBool  *available,PetscReal *val)
273: {


281:   *available = PETSC_FALSE;
282:   if (type!=NORM_1_AND_2) {
283:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,*available);
284:   }
285:   return(0);
286: }

288: /*@
289:    VecNormalize - Normalizes a vector by 2-norm.

291:    Collective on Vec

293:    Input Parameters:
294: +  x - the vector

296:    Output Parameter:
297: .  x - the normalized vector
298: -  val - the vector norm before normalization

300:    Level: intermediate

302: @*/
303: PetscErrorCode  VecNormalize(Vec x,PetscReal *val)
304: {
306:   PetscReal      norm;

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

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

327:    Collective on Vec

329:    Input Parameter:
330: .  x - the vector

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

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

339:    Returns the smallest index with the maximum value
340:    Level: intermediate

342: .seealso: VecNorm(), VecMin()
343: @*/
344: PetscErrorCode  VecMax(Vec x,PetscInt *p,PetscReal *val)
345: {

352:   PetscLogEventBegin(VEC_Max,x,0,0,0);
353:   (*x->ops->max)(x,p,val);
354:   PetscLogEventEnd(VEC_Max,x,0,0,0);
355:   return(0);
356: }

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

361:    Collective on Vec

363:    Input Parameters:
364: .  x - the vector

366:    Output Parameter:
367: +  p - the location of val (pass NULL if you don't want this location)
368: -  val - the minimum component

370:    Level: intermediate

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

375:    This returns the smallest index with the minumum value

377: .seealso: VecMax()
378: @*/
379: PetscErrorCode  VecMin(Vec x,PetscInt *p,PetscReal *val)
380: {

387:   PetscLogEventBegin(VEC_Min,x,0,0,0);
388:   (*x->ops->min)(x,p,val);
389:   PetscLogEventEnd(VEC_Min,x,0,0,0);
390:   return(0);
391: }

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

397:    Collective on Vec

399:    Input Parameters:
400: .  x, y - the vectors

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

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

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

414:    Level: intermediate

416: .seealso: VecDot(), VecMTDot()
417: @*/
418: PetscErrorCode  VecTDot(Vec x,Vec y,PetscScalar *val)
419: {

429:   VecCheckSameSize(x,1,y,2);

431:   PetscLogEventBegin(VEC_TDot,x,y,0,0);
432:   (*x->ops->tdot)(x,y,val);
433:   PetscLogEventEnd(VEC_TDot,x,y,0,0);
434:   return(0);
435: }

437: /*@
438:    VecScale - Scales a vector.

440:    Not collective on Vec

442:    Input Parameters:
443: +  x - the vector
444: -  alpha - the scalar

446:    Output Parameter:
447: .  x - the scaled vector

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

453:    Level: intermediate

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

466:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
467:   PetscLogEventBegin(VEC_Scale,x,0,0,0);
468:   if (alpha != (PetscScalar)1.0) {
469:     VecSetErrorIfLocked(x,1);
470:     /* get current stashed norms */
471:     for (i=0; i<4; i++) {
472:       PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
473:     }
474:     (*x->ops->scale)(x,alpha);
475:     PetscObjectStateIncrease((PetscObject)x);
476:     /* put the scaled stashed norms back into the Vec */
477:     for (i=0; i<4; i++) {
478:       if (flgs[i]) {
479:         PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],PetscAbsScalar(alpha)*norms[i]);
480:       }
481:     }
482:   }
483:   PetscLogEventEnd(VEC_Scale,x,0,0,0);
484:   return(0);
485: }

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

490:    Logically Collective on Vec

492:    Input Parameters:
493: +  x  - the vector
494: -  alpha - the scalar

496:    Output Parameter:
497: .  x  - the vector

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

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

509:    Level: beginner

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

513: @*/
514: PetscErrorCode  VecSet(Vec x,PetscScalar alpha)
515: {
516:   PetscReal      val;

522:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"You cannot call this after you have called VecSetValues() but\n before you have called VecAssemblyBegin/End()");
524:   VecSetErrorIfLocked(x,1);

526:   PetscLogEventBegin(VEC_Set,x,0,0,0);
527:   (*x->ops->set)(x,alpha);
528:   PetscLogEventEnd(VEC_Set,x,0,0,0);
529:   PetscObjectStateIncrease((PetscObject)x);

531:   /*  norms can be simply set (if |alpha|*N not too large) */
532:   val  = PetscAbsScalar(alpha);
533:   if (x->map->N == 0) {
534:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],0.0l);
535:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],0.0);
536:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],0.0);
537:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],0.0);
538:   } else if (val > PETSC_MAX_REAL/x->map->N) {
539:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
540:   } else {
541:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],x->map->N * val);
542:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
543:     val  = PetscSqrtReal((PetscReal)x->map->N) * val;
544:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],val);
545:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],val);
546:   }
547:   return(0);
548: }

550: /*@
551:    VecAXPY - Computes y = alpha x + y.

553:    Logically Collective on Vec

555:    Input Parameters:
556: +  alpha - the scalar
557: -  x, y  - the vectors

559:    Output Parameter:
560: .  y - output vector

562:    Level: intermediate

564:    Notes:
565:     x and y MUST be different vectors
566:     This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine

568: $    VecAXPY(y,alpha,x)                   y = alpha x           +      y
569: $    VecAYPX(y,beta,x)                    y =       x           + beta y
570: $    VecAXPBY(y,alpha,beta,x)             y = alpha x           + beta y
571: $    VecWAXPY(w,alpha,x,y)                w = alpha x           +      y
572: $    VecAXPBYPCZ(w,alpha,beta,gamma,x,y)  z = alpha x           + beta y + gamma z
573: $    VecMAXPY(y,nv,alpha[],x[])           y = sum alpha[i] x[i] +      y

575: .seealso:  VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPBYPCZ(), VecAXPBY()
576: @*/
577: PetscErrorCode  VecAXPY(Vec y,PetscScalar alpha,Vec x)
578: {

587:   VecCheckSameSize(x,3,y,1);
588:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
590:   if (alpha == (PetscScalar)0.0) return(0);
591:   VecSetErrorIfLocked(y,1);

593:   VecLockReadPush(x);
594:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
595:   (*y->ops->axpy)(y,alpha,x);
596:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
597:   VecLockReadPop(x);
598:   PetscObjectStateIncrease((PetscObject)y);
599:   return(0);
600: }

602: /*@
603:    VecAXPBY - Computes y = alpha x + beta y.

605:    Logically Collective on Vec

607:    Input Parameters:
608: +  alpha,beta - the scalars
609: -  x, y  - the vectors

611:    Output Parameter:
612: .  y - output vector

614:    Level: intermediate

616:    Notes:
617:     x and y MUST be different vectors
618:     The implementation is optimized for alpha and/or beta values of 0.0 and 1.0

620: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ()
621: @*/
622: PetscErrorCode  VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
623: {

632:   VecCheckSameSize(y,1,x,4);
633:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
636:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) return(0);
637:   VecSetErrorIfLocked(y,1);
638:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
639:   (*y->ops->axpby)(y,alpha,beta,x);
640:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
641:   PetscObjectStateIncrease((PetscObject)y);
642:   return(0);
643: }

645: /*@
646:    VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z

648:    Logically Collective on Vec

650:    Input Parameters:
651: +  alpha,beta, gamma - the scalars
652: -  x, y, z  - the vectors

654:    Output Parameter:
655: .  z - output vector

657:    Level: intermediate

659:    Notes:
660:     x, y and z must be different vectors
661:     The implementation is optimized for alpha of 1.0 and gamma of 1.0 or 0.0

663: .seealso:  VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBY()
664: @*/
665: PetscErrorCode  VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
666: {

678:   VecCheckSameSize(x,1,y,5);
679:   VecCheckSameSize(x,1,z,6);
680:   if (x == y || x == z) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
681:   if (y == z) SETERRQ(PetscObjectComm((PetscObject)y),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
685:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) return(0);
686:   VecSetErrorIfLocked(z,1);

688:   PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
689:   (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
690:   PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
691:   PetscObjectStateIncrease((PetscObject)z);
692:   return(0);
693: }

695: /*@
696:    VecAYPX - Computes y = x + beta y.

698:    Logically Collective on Vec

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

704:    Output Parameter:
705: .  y - output vector

707:    Level: intermediate

709:    Notes:
710:     x and y MUST be different vectors
711:     The implementation is optimized for beta of -1.0, 0.0, and 1.0

713: .seealso:  VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
714: @*/
715: PetscErrorCode  VecAYPX(Vec y,PetscScalar beta,Vec x)
716: {

725:   VecCheckSameSize(x,1,y,3);
726:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
728:   VecSetErrorIfLocked(y,1);

730:   PetscLogEventBegin(VEC_AYPX,x,y,0,0);
731:    (*y->ops->aypx)(y,beta,x);
732:   PetscLogEventEnd(VEC_AYPX,x,y,0,0);
733:   PetscObjectStateIncrease((PetscObject)y);
734:   return(0);
735: }

737: /*@
738:    VecWAXPY - Computes w = alpha x + y.

740:    Logically Collective on Vec

742:    Input Parameters:
743: +  alpha - the scalar
744: -  x, y  - the vectors

746:    Output Parameter:
747: .  w - the result

749:    Level: intermediate

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

755: .seealso: VecAXPY(), VecAYPX(), VecAXPBY(), VecMAXPY(), VecAXPBYPCZ()
756: @*/
757: PetscErrorCode  VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
758: {

770:   VecCheckSameSize(x,3,y,4);
771:   VecCheckSameSize(x,3,w,1);
772:   if (w == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector y, suggest VecAXPY()");
773:   if (w == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector x, suggest VecAYPX()");
775:   VecSetErrorIfLocked(w,1);

777:   PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
778:    (*w->ops->waxpy)(w,alpha,x,y);
779:   PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
780:   PetscObjectStateIncrease((PetscObject)w);
781:   return(0);
782: }

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

787:    Not Collective

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

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

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

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

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

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

816:    Level: beginner

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

827:   if (!ni) return(0);

832:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
833:   (*x->ops->setvalues)(x,ni,ix,y,iora);
834:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
835:   PetscObjectStateIncrease((PetscObject)x);
836:   return(0);
837: }

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

843:     Not Collective

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

850:    Output Parameter:
851: .   y - array of values

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

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

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

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

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

866:    Level: beginner

868: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues()
869: @*/
870: PetscErrorCode  VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
871: {

876:   if (!ni) return(0);
880:   (*x->ops->getvalues)(x,ni,ix,y);
881:   return(0);
882: }

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

887:    Not Collective

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

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

902:    Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
903:    options cannot be mixed without intervening calls to the assembly
904:    routines.

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

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

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

916:    Level: intermediate

918: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
919:            VecSetValues()
920: @*/
921: PetscErrorCode  VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
922: {

927:   if (!ni) return(0);

932:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
933:   (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
934:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
935:   PetscObjectStateIncrease((PetscObject)x);
936:   return(0);
937: }

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

943:    Not Collective

945:    Input Parameters:
946: +  x - vector to insert in
947: .  ni - number of elements to add
948: .  ix - indices where to add
949: .  y - array of values
950: -  iora - either INSERT_VALUES or ADD_VALUES, where
951:    ADD_VALUES adds values to any existing entries, and
952:    INSERT_VALUES replaces existing entries with new values

954:    Level: intermediate

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

959:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
960:    options cannot be mixed without intervening calls to the assembly
961:    routines.

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

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

968: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetLocalToGlobalMapping(),
969:            VecSetValuesBlockedLocal()
970: @*/
971: PetscErrorCode  VecSetValuesLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
972: {
974:   PetscInt       lixp[128],*lix = lixp;

978:   if (!ni) return(0);

983:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
984:   if (!x->ops->setvalueslocal) {
985:     if (!x->map->mapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMapping()");
986:     if (ni > 128) {
987:       PetscMalloc1(ni,&lix);
988:     }
989:     ISLocalToGlobalMappingApply(x->map->mapping,ni,(PetscInt*)ix,lix);
990:     (*x->ops->setvalues)(x,ni,lix,y,iora);
991:     if (ni > 128) {
992:       PetscFree(lix);
993:     }
994:   } else {
995:     (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
996:   }
997:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
998:   PetscObjectStateIncrease((PetscObject)x);
999:   return(0);
1000: }

1002: /*@
1003:    VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1004:    using a local ordering of the nodes.

1006:    Not Collective

1008:    Input Parameters:
1009: +  x - vector to insert in
1010: .  ni - number of blocks to add
1011: .  ix - indices where to add in block count, not element count
1012: .  y - array of values
1013: -  iora - either INSERT_VALUES or ADD_VALUES, where
1014:    ADD_VALUES adds values to any existing entries, and
1015:    INSERT_VALUES replaces existing entries with new values

1017:    Level: intermediate

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

1023:    Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1024:    options cannot be mixed without intervening calls to the assembly
1025:    routines.

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

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

1032: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
1033:            VecSetLocalToGlobalMapping()
1034: @*/
1035: PetscErrorCode  VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1036: {
1038:   PetscInt       lixp[128],*lix = lixp;

1042:   if (!ni) return(0);
1046:   if (ni > 128) {
1047:     PetscMalloc1(ni,&lix);
1048:   }

1050:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1051:   ISLocalToGlobalMappingApplyBlock(x->map->mapping,ni,(PetscInt*)ix,lix);
1052:   (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1053:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1054:   if (ni > 128) {
1055:     PetscFree(lix);
1056:   }
1057:   PetscObjectStateIncrease((PetscObject)x);
1058:   return(0);
1059: }

1061: /*@
1062:    VecMTDot - Computes indefinite vector multiple dot products.
1063:    That is, it does NOT use the complex conjugate.

1065:    Collective on Vec

1067:    Input Parameters:
1068: +  x - one vector
1069: .  nv - number of vectors
1070: -  y - array of vectors.  Note that vectors are pointers

1072:    Output Parameter:
1073: .  val - array of the dot products

1075:    Notes for Users of Complex Numbers:
1076:    For complex vectors, VecMTDot() computes the indefinite form
1077: $      val = (x,y) = y^T x,
1078:    where y^T denotes the transpose of y.

1080:    Use VecMDot() for the inner product
1081: $      val = (x,y) = y^H x,
1082:    where y^H denotes the conjugate transpose of y.

1084:    Level: intermediate

1086: .seealso: VecMDot(), VecTDot()
1087: @*/
1088: PetscErrorCode  VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1089: {

1095:   if (!nv) return(0);
1102:   VecCheckSameSize(x,1,*y,3);

1104:   PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1105:   (*x->ops->mtdot)(x,nv,y,val);
1106:   PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1107:   return(0);
1108: }

1110: /*@
1111:    VecMDot - Computes vector multiple dot products.

1113:    Collective on Vec

1115:    Input Parameters:
1116: +  x - one vector
1117: .  nv - number of vectors
1118: -  y - array of vectors.

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

1123:    Notes for Users of Complex Numbers:
1124:    For complex vectors, VecMDot() computes
1125: $     val = (x,y) = y^H x,
1126:    where y^H denotes the conjugate transpose of y.

1128:    Use VecMTDot() for the indefinite form
1129: $     val = (x,y) = y^T x,
1130:    where y^T denotes the transpose of y.

1132:    Level: intermediate

1134: .seealso: VecMTDot(), VecDot()
1135: @*/
1136: PetscErrorCode  VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1137: {

1143:   if (!nv) return(0);
1144:   if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1151:   VecCheckSameSize(x,1,*y,3);

1153:   PetscLogEventBegin(VEC_MDot,x,*y,0,0);
1154:   (*x->ops->mdot)(x,nv,y,val);
1155:   PetscLogEventEnd(VEC_MDot,x,*y,0,0);
1156:   return(0);
1157: }

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

1162:    Logically Collective on Vec

1164:    Input Parameters:
1165: +  nv - number of scalars and x-vectors
1166: .  alpha - array of scalars
1167: .  y - one vector
1168: -  x - array of vectors

1170:    Level: intermediate

1172:    Notes:
1173:     y cannot be any of the x vectors

1175: .seealso:  VecAYPX(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
1176: @*/
1177: PetscErrorCode  VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1178: {
1180:   PetscInt       i;
1181:   PetscBool      nonzero;

1186:   if (!nv) return(0);
1187:   if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1194:   VecCheckSameSize(y,1,*x,4);
1196:   for (i=0, nonzero = PETSC_FALSE; i<nv && !nonzero; i++) nonzero = (PetscBool)(nonzero || alpha[i] != (PetscScalar)0.0);
1197:   if (!nonzero) return(0);
1198:   VecSetErrorIfLocked(y,1);
1199:   PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1200:   (*y->ops->maxpy)(y,nv,alpha,x);
1201:   PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1202:   PetscObjectStateIncrease((PetscObject)y);
1203:   return(0);
1204: }

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

1211:    Collective on X

1213:    Input Arguments:
1214: +  nx   - number of vectors to be concatenated
1215: -  X    - array containing the vectors to be concatenated in the order of concatenation

1217:    Output Arguments:
1218: +  Y    - concatenated vector
1219: -  x_is - array of index sets corresponding to the concatenated components of Y (NULL if not needed)

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

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

1231:    Level: advanced

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


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

1289: /*@
1290:    VecGetSubVector - Gets a vector representing part of another vector

1292:    Collective on X and IS

1294:    Input Arguments:
1295: + X - vector from which to extract a subvector
1296: - is - index set representing portion of X to extract

1298:    Output Arguments:
1299: . Y - subvector corresponding to is

1301:    Level: advanced

1303:    Notes:
1304:    The subvector Y should be returned with VecRestoreSubVector().
1305:    X and is must be defined on the same communicator

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

1311: .seealso: MatCreateSubMatrix()
1312: @*/
1313: PetscErrorCode  VecGetSubVector(Vec X,IS is,Vec *Y)
1314: {
1315:   PetscErrorCode   ierr;
1316:   Vec              Z;

1323:   if (X->ops->getsubvector) {
1324:     (*X->ops->getsubvector)(X,is,&Z);
1325:   } else { /* Default implementation currently does no caching */
1326:     PetscInt  gstart,gend,start;
1327:     PetscBool red[2] = { PETSC_TRUE, PETSC_TRUE };
1328:     PetscInt  n,N,ibs,vbs,bs = -1;

1330:     ISGetLocalSize(is,&n);
1331:     ISGetSize(is,&N);
1332:     ISGetBlockSize(is,&ibs);
1333:     VecGetBlockSize(X,&vbs);
1334:     VecGetOwnershipRange(X,&gstart,&gend);
1335:     ISContiguousLocal(is,gstart,gend,&start,&red[0]);
1336:     /* block size is given by IS if ibs > 1; otherwise, check the vector */
1337:     if (ibs > 1) {
1338:       MPIU_Allreduce(MPI_IN_PLACE,red,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1339:       bs   = ibs;
1340:     } else {
1341:       if (n%vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1342:       MPIU_Allreduce(MPI_IN_PLACE,red,2,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1343:       if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1344:     }
1345:     if (red[0]) { /* We can do a no-copy implementation */
1346:       const PetscScalar *x;
1347:       PetscInt          state = 0;
1348:       PetscBool         isstd,iscuda,iship;

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

1359:         VecCUDAGetArrays_Private(X,&x,&x_d,&flg);
1360:         if (flg == PETSC_OFFLOAD_UNALLOCATED) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for PETSC_OFFLOAD_UNALLOCATED");
1361:         if (n && !x && !x_d) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Missing vector data");
1362:         if (x) x += start;
1363:         if (x_d) x_d += start;
1364:         MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1365:         if (size == 1) {
1366:           VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X),bs,n,x,x_d,&Z);
1367:         } else {
1368:           VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X),bs,n,N,x,x_d,&Z);
1369:         }
1370:         Z->offloadmask = flg;
1371: #endif
1372:       } else if (iship) {
1373: #if defined(PETSC_HAVE_HIP)
1374:         const PetscScalar *x_d;
1375:         PetscMPIInt       size;
1376:         PetscOffloadMask  flg;

1378:         VecHIPGetArrays_Private(X,&x,&x_d,&flg);
1379:         if (flg == PETSC_OFFLOAD_UNALLOCATED) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for PETSC_OFFLOAD_UNALLOCATED");
1380:         if (n && !x && !x_d) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Missing vector data");
1381:         if (x) x += start;
1382:         if (x_d) x_d += start;
1383:         MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1384:         if (size == 1) {
1385:           VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X),bs,n,x,x_d,&Z);
1386:         } else {
1387:           VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X),bs,n,N,x,x_d,&Z);
1388:         }
1389:         Z->offloadmask = flg;
1390: #endif
1391:       } else if (isstd) {
1392:         PetscMPIInt size;

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

1413:       /* this is relevant only in debug mode */
1414:       VecLockGet(X,&state);
1415:       if (state) {
1416:         VecLockReadPush(Z);
1417:       }
1418:       Z->ops->placearray = NULL;
1419:       Z->ops->replacearray = NULL;
1420:     } else { /* Have to create a scatter and do a copy */
1421:       VecScatter scatter;

1423:       VecCreate(PetscObjectComm((PetscObject)is),&Z);
1424:       VecSetSizes(Z,n,N);
1425:       VecSetBlockSize(Z,bs);
1426:       VecSetType(Z,((PetscObject)X)->type_name);
1427:       VecScatterCreate(X,is,Z,NULL,&scatter);
1428:       VecScatterBegin(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1429:       VecScatterEnd(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1430:       PetscObjectCompose((PetscObject)Z,"VecGetSubVector_Scatter",(PetscObject)scatter);
1431:       VecScatterDestroy(&scatter);
1432:     }
1433:   }
1434:   /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1435:   if (VecGetSubVectorSavedStateId < 0) {PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);}
1436:   PetscObjectComposedDataSetInt((PetscObject)Z,VecGetSubVectorSavedStateId,1);
1437:   *Y   = Z;
1438:   return(0);
1439: }

1441: /*@
1442:    VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()

1444:    Collective on IS

1446:    Input Arguments:
1447: + X - vector from which subvector was obtained
1448: . is - index set representing the subset of X
1449: - Y - subvector being restored

1451:    Level: advanced

1453: .seealso: VecGetSubVector()
1454: @*/
1455: PetscErrorCode  VecRestoreSubVector(Vec X,IS is,Vec *Y)
1456: {

1465:   if (X->ops->restoresubvector) {
1466:     (*X->ops->restoresubvector)(X,is,Y);
1467:   } else {
1468:     PETSC_UNUSED PetscObjectState dummystate = 0;
1469:     PetscBool valid;

1471:     PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,dummystate,valid);
1472:     if (!valid) {
1473:       VecScatter scatter;
1474:       PetscInt   state;

1476:       VecLockGet(X,&state);
1477:       if (state != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vec X is locked for read-only or read/write access");

1479:       PetscObjectQuery((PetscObject)*Y,"VecGetSubVector_Scatter",(PetscObject*)&scatter);
1480:       if (scatter) {
1481:         VecScatterBegin(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1482:         VecScatterEnd(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1483:       } else {
1484:         PetscBool         iscuda,iship;
1485:         PetscObjectTypeCompareAny((PetscObject)X,&iscuda,VECSEQCUDA,VECMPICUDA,"");
1486:         PetscObjectTypeCompareAny((PetscObject)X,&iship,VECSEQHIP,VECMPIHIP,"");

1488:         if (iscuda) {
1489: #if defined(PETSC_HAVE_CUDA)
1490:           PetscOffloadMask ymask = (*Y)->offloadmask;

1492:           /* The offloadmask of X dictates where to move memory
1493:              If X GPU data is valid, then move Y data on GPU if needed
1494:              Otherwise, move back to the CPU */
1495:           switch (X->offloadmask) {
1496:           case PETSC_OFFLOAD_BOTH:
1497:             if (ymask == PETSC_OFFLOAD_CPU) {
1498:               VecCUDAResetArray(*Y);
1499:             } else if (ymask == PETSC_OFFLOAD_GPU) {
1500:               X->offloadmask = PETSC_OFFLOAD_GPU;
1501:             }
1502:             break;
1503:           case PETSC_OFFLOAD_GPU:
1504:             if (ymask == PETSC_OFFLOAD_CPU) {
1505:               VecCUDAResetArray(*Y);
1506:             }
1507:             break;
1508:           case PETSC_OFFLOAD_CPU:
1509:             if (ymask == PETSC_OFFLOAD_GPU) {
1510:               VecResetArray(*Y);
1511:             }
1512:             break;
1513:           case PETSC_OFFLOAD_UNALLOCATED:
1514:           case PETSC_OFFLOAD_VECKOKKOS:
1515:             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"This should not happen");
1516:             break;
1517:           }
1518: #endif
1519:         } else if (iship) {
1520: #if defined(PETSC_HAVE_HIP)
1521:           PetscOffloadMask ymask = (*Y)->offloadmask;

1523:           /* The offloadmask of X dictates where to move memory
1524:              If X GPU data is valid, then move Y data on GPU if needed
1525:              Otherwise, move back to the CPU */
1526:           switch (X->offloadmask) {
1527:           case PETSC_OFFLOAD_BOTH:
1528:             if (ymask == PETSC_OFFLOAD_CPU) {
1529:               VecHIPResetArray(*Y);
1530:             } else if (ymask == PETSC_OFFLOAD_GPU) {
1531:               X->offloadmask = PETSC_OFFLOAD_GPU;
1532:             }
1533:             break;
1534:           case PETSC_OFFLOAD_GPU:
1535:             if (ymask == PETSC_OFFLOAD_CPU) {
1536:               VecHIPResetArray(*Y);
1537:             }
1538:             break;
1539:           case PETSC_OFFLOAD_CPU:
1540:             if (ymask == PETSC_OFFLOAD_GPU) {
1541:               VecResetArray(*Y);
1542:             }
1543:             break;
1544:           case PETSC_OFFLOAD_UNALLOCATED:
1545:           case PETSC_OFFLOAD_VECKOKKOS:
1546:             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"This should not happen");
1547:             break;
1548:           }
1549: #endif
1550:         } else {
1551:           /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1552:           VecResetArray(*Y);
1553:         }
1554:         PetscObjectStateIncrease((PetscObject)X);
1555:       }
1556:     }
1557:     VecDestroy(Y);
1558:   }
1559:   return(0);
1560: }

1562: /*@
1563:    VecGetLocalVectorRead - Maps the local portion of a vector into a
1564:    vector.  You must call VecRestoreLocalVectorRead() when the local
1565:    vector is no longer needed.

1567:    Not collective.

1569:    Input parameter:
1570: .  v - The vector for which the local vector is desired.

1572:    Output parameter:
1573: .  w - Upon exit this contains the local vector.

1575:    Level: beginner

1577:    Notes:
1578:    This function is similar to VecGetArrayRead() which maps the local
1579:    portion into a raw pointer.  VecGetLocalVectorRead() is usually
1580:    almost as efficient as VecGetArrayRead() but in certain circumstances
1581:    VecGetLocalVectorRead() can be much more efficient than
1582:    VecGetArrayRead().  This is because the construction of a contiguous
1583:    array representing the vector data required by VecGetArrayRead() can
1584:    be an expensive operation for certain vector types.  For example, for
1585:    GPU vectors VecGetArrayRead() requires that the data between device
1586:    and host is synchronized.

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

1591: .seealso: VecRestoreLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1592: @*/
1593: PetscErrorCode VecGetLocalVectorRead(Vec v,Vec w)
1594: {
1596:   PetscScalar    *a;

1601:   VecCheckSameLocalSize(v,1,w,2);
1602:   if (v->ops->getlocalvectorread) {
1603:     (*v->ops->getlocalvectorread)(v,w);
1604:   } else {
1605:     VecGetArrayRead(v,(const PetscScalar**)&a);
1606:     VecPlaceArray(w,a);
1607:   }
1608:   PetscObjectStateIncrease((PetscObject)w);
1609:   VecLockReadPush(v);
1610:   VecLockReadPush(w);
1611:   return(0);
1612: }

1614: /*@
1615:    VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1616:    previously mapped into a vector using VecGetLocalVectorRead().

1618:    Not collective.

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

1624:    Level: beginner

1626: .seealso: VecGetLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1627: @*/
1628: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec w)
1629: {
1631:   PetscScalar    *a;

1636:   if (v->ops->restorelocalvectorread) {
1637:     (*v->ops->restorelocalvectorread)(v,w);
1638:   } else {
1639:     VecGetArrayRead(w,(const PetscScalar**)&a);
1640:     VecRestoreArrayRead(v,(const PetscScalar**)&a);
1641:     VecResetArray(w);
1642:   }
1643:   VecLockReadPop(v);
1644:   VecLockReadPop(w);
1645:   PetscObjectStateIncrease((PetscObject)w);
1646:   return(0);
1647: }

1649: /*@
1650:    VecGetLocalVector - Maps the local portion of a vector into a
1651:    vector.

1653:    Collective on v, not collective on w.

1655:    Input parameter:
1656: .  v - The vector for which the local vector is desired.

1658:    Output parameter:
1659: .  w - Upon exit this contains the local vector.

1661:    Level: beginner

1663:    Notes:
1664:    This function is similar to VecGetArray() which maps the local
1665:    portion into a raw pointer.  VecGetLocalVector() is usually about as
1666:    efficient as VecGetArray() but in certain circumstances
1667:    VecGetLocalVector() can be much more efficient than VecGetArray().
1668:    This is because the construction of a contiguous array representing
1669:    the vector data required by VecGetArray() can be an expensive
1670:    operation for certain vector types.  For example, for GPU vectors
1671:    VecGetArray() requires that the data between device and host is
1672:    synchronized.

1674: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1675: @*/
1676: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1677: {
1679:   PetscScalar    *a;

1684:   VecCheckSameLocalSize(v,1,w,2);
1685:   if (v->ops->getlocalvector) {
1686:     (*v->ops->getlocalvector)(v,w);
1687:   } else {
1688:     VecGetArray(v,&a);
1689:     VecPlaceArray(w,a);
1690:   }
1691:   PetscObjectStateIncrease((PetscObject)w);
1692:   return(0);
1693: }

1695: /*@
1696:    VecRestoreLocalVector - Unmaps the local portion of a vector
1697:    previously mapped into a vector using VecGetLocalVector().

1699:    Logically collective.

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

1705:    Level: beginner

1707: .seealso: VecGetLocalVector(), VecGetLocalVectorRead(), VecRestoreLocalVectorRead(), LocalVectorRead(), VecGetArrayRead(), VecGetArray()
1708: @*/
1709: PetscErrorCode VecRestoreLocalVector(Vec v,Vec w)
1710: {
1712:   PetscScalar    *a;

1717:   if (v->ops->restorelocalvector) {
1718:     (*v->ops->restorelocalvector)(v,w);
1719:   } else {
1720:     VecGetArray(w,&a);
1721:     VecRestoreArray(v,&a);
1722:     VecResetArray(w);
1723:   }
1724:   PetscObjectStateIncrease((PetscObject)w);
1725:   PetscObjectStateIncrease((PetscObject)v);
1726:   return(0);
1727: }

1729: /*@C
1730:    VecGetArray - Returns a pointer to a contiguous array that contains this
1731:    processor's portion of the vector data. For the standard PETSc
1732:    vectors, VecGetArray() returns a pointer to the local data array and
1733:    does not use any copies. If the underlying vector data is not stored
1734:    in a contiguous array this routine will copy the data to a contiguous
1735:    array and return a pointer to that. You MUST call VecRestoreArray()
1736:    when you no longer need access to the array.

1738:    Logically Collective on Vec

1740:    Input Parameter:
1741: .  x - the vector

1743:    Output Parameter:
1744: .  a - location to put pointer to the array

1746:    Fortran Note:
1747:    This routine is used differently from Fortran 77
1748: $    Vec         x
1749: $    PetscScalar x_array(1)
1750: $    PetscOffset i_x
1751: $    PetscErrorCode ierr
1752: $       call VecGetArray(x,x_array,i_x,ierr)
1753: $
1754: $   Access first local entry in vector with
1755: $      value = x_array(i_x + 1)
1756: $
1757: $      ...... other code
1758: $       call VecRestoreArray(x,x_array,i_x,ierr)
1759:    For Fortran 90 see VecGetArrayF90()

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

1764:    Level: beginner

1766: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1767:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1768: @*/
1769: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1770: {

1775:   VecSetErrorIfLocked(x,1);
1776:   if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
1777:     (*x->ops->getarray)(x,a);
1778:   } else if (x->petscnative) { /* VECSTANDARD */
1779:     *a = *((PetscScalar**)x->data);
1780:   } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array for vector type \"%s\"",((PetscObject)x)->type_name);
1781:   return(0);
1782: }

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

1787:    Logically Collective on Vec

1789:    Input Parameters:
1790: +  x - the vector
1791: -  a - location of pointer to array obtained from VecGetArray()

1793:    Level: beginner

1795: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1796:           VecGetArrayPair(), VecRestoreArrayPair()
1797: @*/
1798: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1799: {

1804:   if (x->ops->restorearray) { /* VECNEST, VECCUDA etc */
1805:     (*x->ops->restorearray)(x,a);
1806:   } else if (x->petscnative) { /* VECSTANDARD */
1807:     /* nothing */
1808:   } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot restore array for vector type \"%s\"",((PetscObject)x)->type_name);
1809:   if (a) *a = NULL;
1810:   PetscObjectStateIncrease((PetscObject)x);
1811:   return(0);
1812: }
1813: /*@C
1814:    VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.

1816:    Not Collective

1818:    Input Parameter:
1819: .  x - the vector

1821:    Output Parameter:
1822: .  a - the array

1824:    Level: beginner

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

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

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

1835: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1836: @*/
1837: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1838: {

1843:   if (x->ops->getarray) { /* VECNEST, VECCUDA, VECKOKKOS etc */
1844:     (*x->ops->getarray)(x,(PetscScalar**)a);
1845:   } else if (x->petscnative) { /* VECSTANDARD */
1846:     *a = *((PetscScalar**)x->data);
1847:   } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array read for vector type \"%s\"",((PetscObject)x)->type_name);
1848:   return(0);
1849: }

1851: /*@C
1852:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

1854:    Not Collective

1856:    Input Parameters:
1857: +  vec - the vector
1858: -  array - the array

1860:    Level: beginner

1862: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1863: @*/
1864: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1865: {

1870:   if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
1871:     /* nothing */
1872:   } else if (x->ops->restorearrayread) { /* VECNEST */
1873:     (*x->ops->restorearrayread)(x,a);
1874:   } else { /* No one? */
1875:     (*x->ops->restorearray)(x,(PetscScalar**)a);
1876:   }
1877:   if (a) *a = NULL;
1878:   return(0);
1879: }

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

1886:    Logically Collective on Vec

1888:    Input Parameter:
1889: .  x - the vector

1891:    Output Parameter:
1892: .  a - location to put pointer to the array

1894:    Level: intermediate

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

1899:    Concepts: vector^accessing local values

1901: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1902:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArray(), VecRestoreArrayWrite()
1903: @*/
1904: PetscErrorCode VecGetArrayWrite(Vec x,PetscScalar **a)
1905: {

1910:   VecSetErrorIfLocked(x,1);
1911:   if (x->ops->getarraywrite) {
1912:     (*x->ops->getarraywrite)(x,a);
1913:   } else {
1914:     VecGetArray(x,a);
1915:   }
1916:   return(0);
1917: }

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

1922:    Logically Collective on Vec

1924:    Input Parameters:
1925: +  x - the vector
1926: -  a - location of pointer to array obtained from VecGetArray()

1928:    Level: beginner

1930: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1931:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite()
1932: @*/
1933: PetscErrorCode VecRestoreArrayWrite(Vec x,PetscScalar **a)
1934: {

1939:   if (x->ops->restorearraywrite) {
1940:     (*x->ops->restorearraywrite)(x,a);
1941:   } else if (x->ops->restorearray) {
1942:     (*x->ops->restorearray)(x,a);
1943:   }
1944:   if (a) *a = NULL;
1945:   PetscObjectStateIncrease((PetscObject)x);
1946:   return(0);
1947: }

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

1954:    Logically Collective on Vec

1956:    Input Parameters:
1957: +  x - the vectors
1958: -  n - the number of vectors

1960:    Output Parameter:
1961: .  a - location to put pointer to the array

1963:    Fortran Note:
1964:    This routine is not supported in Fortran.

1966:    Level: intermediate

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

1980:   if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1981:   PetscMalloc1(n,&q);
1982:   for (i=0; i<n; ++i) {
1983:     VecGetArray(x[i],&q[i]);
1984:   }
1985:   *a = q;
1986:   return(0);
1987: }

1989: /*@C
1990:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1991:    has been called.

1993:    Logically Collective on Vec

1995:    Input Parameters:
1996: +  x - the vector
1997: .  n - the number of vectors
1998: -  a - location of pointer to arrays obtained from VecGetArrays()

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

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

2009:    Level: intermediate

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


2024:   for (i=0; i<n; ++i) {
2025:     VecRestoreArray(x[i],&q[i]);
2026:   }
2027:   PetscFree(q);
2028:   return(0);
2029: }

2031: /*@C
2032:    VecGetArrayAndMemType - Like VecGetArray(), but if this is a GPU vector and it is currently offloaded to GPU,
2033:    the returned pointer will be a GPU pointer to the GPU memory that contains this processor's portion of the
2034:    vector data. Otherwise, it functions as VecGetArray().

2036:    Logically Collective on Vec

2038:    Input Parameter:
2039: .  x - the vector

2041:    Output Parameters:
2042: +  a - location to put pointer to the array
2043: -  mtype - memory type of the array

2045:    Level: beginner

2047: .seealso: VecRestoreArrayAndMemType(), VecRestoreArrayAndMemType(), VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(),
2048:           VecPlaceArray(), VecGetArray2d(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
2049: @*/
2050: PetscErrorCode VecGetArrayAndMemType(Vec x,PetscScalar **a,PetscMemType *mtype)
2051: {

2057:   VecSetErrorIfLocked(x,1);
2058:   if (x->ops->getarrayandmemtype) { /* VECCUDA, VECKOKKOS etc */
2059:     (*x->ops->getarrayandmemtype)(x,a,mtype);
2060:   } else { /* VECSTANDARD, VECNEST, VECVIENNACL */
2061:     VecGetArray(x,a);
2062:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2063:   }
2064:   return(0);
2065: }

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

2070:    Logically Collective on Vec

2072:    Input Parameters:
2073: +  x - the vector
2074: -  a - location of pointer to array obtained from VecGetArrayAndMemType()

2076:    Level: beginner

2078: .seealso: VecGetArrayAndMemType(), VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(),
2079:           VecPlaceArray(), VecRestoreArray2d(), VecGetArrayPair(), VecRestoreArrayPair()
2080: @*/
2081: PetscErrorCode VecRestoreArrayAndMemType(Vec x,PetscScalar **a)
2082: {

2088:   if (x->ops->restorearrayandmemtype) { /* VECCUDA, VECKOKKOS etc */
2089:     (*x->ops->restorearrayandmemtype)(x,a);
2090:   } else if (x->ops->restorearray) { /* VECNEST, VECVIENNACL */
2091:     (*x->ops->restorearray)(x,a);
2092:   } /* VECSTANDARD does nothing */
2093:   if (a) *a = NULL;
2094:   PetscObjectStateIncrease((PetscObject)x);
2095:   return(0);
2096: }

2098: /*@C
2099:    VecGetArrayReadAndMemType - Like VecGetArrayRead(), but if this is a CUDA vector and it is currently offloaded to GPU,
2100:    the returned pointer will be a GPU pointer to the GPU memory that contains this processor's portion of the
2101:    vector data. Otherwise, it functions as VecGetArrayRead().

2103:    Not Collective

2105:    Input Parameter:
2106: .  x - the vector

2108:    Output Parameters:
2109: +  a - the array
2110: -  mtype - memory type of the array

2112:    Level: beginner

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

2117: .seealso: VecRestoreArrayReadAndMemType(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayAndMemType()
2118: @*/
2119: PetscErrorCode VecGetArrayReadAndMemType(Vec x,const PetscScalar **a,PetscMemType *mtype)
2120: {

2126:  #if defined(PETSC_USE_DEBUG)
2127:   if (x->ops->getarrayreadandmemtype) SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Not expected vector type \"%s\" has ops->getarrayreadandmemtype",((PetscObject)x)->type_name);
2128:  #endif

2130:   if (x->ops->getarrayandmemtype) { /* VECCUDA, VECKOKKOS etc, though they are also petscnative */
2131:     (*x->ops->getarrayandmemtype)(x,(PetscScalar**)a,mtype);
2132:   } else if (x->ops->getarray) { /* VECNEST, VECVIENNACL */
2133:     (*x->ops->getarray)(x,(PetscScalar**)a);
2134:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2135:   } else if (x->petscnative) { /* VECSTANDARD */
2136:     *a = *((PetscScalar**)x->data);
2137:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2138:   } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array read in place for vector type \"%s\"",((PetscObject)x)->type_name);
2139:   return(0);
2140: }

2142: /*@C
2143:    VecRestoreArrayReadAndMemType - Restore array obtained with VecGetArrayReadAndMemType()

2145:    Not Collective

2147:    Input Parameters:
2148: +  vec - the vector
2149: -  array - the array

2151:    Level: beginner

2153: .seealso: VecGetArrayReadAndMemType(), VecRestoreArrayAndMemType(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
2154: @*/
2155: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x,const PetscScalar **a)
2156: {

2162:   if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS, VECVIENNACL etc */
2163:     /* nothing */
2164:   } else if (x->ops->restorearrayread) { /* VECNEST */
2165:     (*x->ops->restorearrayread)(x,a);
2166:   } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot restore array read in place for vector type \"%s\"",((PetscObject)x)->type_name);
2167:   if (a) *a = NULL;
2168:   return(0);
2169: }

2171: /*@
2172:    VecPlaceArray - Allows one to replace the array in a vector with an
2173:    array provided by the user. This is useful to avoid copying an array
2174:    into a vector.

2176:    Not Collective

2178:    Input Parameters:
2179: +  vec - the vector
2180: -  array - the array

2182:    Notes:
2183:    You can return to the original array with a call to VecResetArray()

2185:    Level: developer

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

2189: @*/
2190: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
2191: {

2198:   if (vec->ops->placearray) {
2199:     (*vec->ops->placearray)(vec,array);
2200:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
2201:   PetscObjectStateIncrease((PetscObject)vec);
2202:   return(0);
2203: }

2205: /*@C
2206:    VecReplaceArray - Allows one to replace the array in a vector with an
2207:    array provided by the user. This is useful to avoid copying an array
2208:    into a vector.

2210:    Not Collective

2212:    Input Parameters:
2213: +  vec - the vector
2214: -  array - the array

2216:    Notes:
2217:    This permanently replaces the array and frees the memory associated
2218:    with the old array.

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

2223:    Not supported from Fortran

2225:    Level: developer

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

2229: @*/
2230: PetscErrorCode  VecReplaceArray(Vec vec,const PetscScalar array[])
2231: {

2237:   if (vec->ops->replacearray) {
2238:     (*vec->ops->replacearray)(vec,array);
2239:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
2240:   PetscObjectStateIncrease((PetscObject)vec);
2241:   return(0);
2242: }

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

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

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

2262:    Input Parameter:
2263: .  v - the vector

2265:    Output Parameter:
2266: .  a - the CUDA device pointer

2268:    Fortran note:
2269:    This function is not currently available from Fortran.

2271:    Level: intermediate

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

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

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

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

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

2304:    Level: intermediate

2306: .seealso: VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2307: @*/
2308: PETSC_EXTERN PetscErrorCode VecCUDARestoreArray(Vec v, PetscScalar **a)
2309: {

2314: #if defined(PETSC_HAVE_CUDA)
2315:   v->offloadmask = PETSC_OFFLOAD_GPU;
2316: #endif
2317:   PetscObjectStateIncrease((PetscObject)v);
2318:   return(0);
2319: }

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

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

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

2339:    Input Parameter:
2340: .  v - the vector

2342:    Output Parameter:
2343: .  a - the CUDA pointer.

2345:    Fortran note:
2346:    This function is not currently available from Fortran.

2348:    Level: intermediate

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

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

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

2368:    Input Parameter:
2369: +  v - the vector
2370: -  a - the CUDA device pointer.  This pointer is invalid after
2371:        VecCUDARestoreArrayRead() returns.

2373:    Fortran note:
2374:    This function is not currently available from Fortran.

2376:    Level: intermediate

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

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

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

2397:    The device pointer needs to be released with
2398:    VecCUDARestoreArrayWrite().  When the pointer is released the
2399:    host data of the vector is marked as out of data.  Subsequent access
2400:    of the host data with e.g. VecGetArray() incurs a device to host data
2401:    transfer.

2403:    Input Parameter:
2404: .  v - the vector

2406:    Output Parameter:
2407: .  a - the CUDA pointer

2409:    Fortran note:
2410:    This function is not currently available from Fortran.

2412:    Level: advanced

2414: .seealso: VecCUDARestoreArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2415: @*/
2416: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayWrite(Vec v, PetscScalar **a)
2417: {
2420:  #if defined(PETSC_HAVE_CUDA)
2421:   {
2423:     VecCUDAAllocateCheck(v);
2424:     *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2425:   }
2426:  #endif
2427:   return(0);
2428: }

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

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

2437:    Input Parameter:
2438: +  v - the vector
2439: -  a - the CUDA device pointer.  This pointer is invalid after
2440:        VecCUDARestoreArrayWrite() returns.

2442:    Fortran note:
2443:    This function is not currently available from Fortran.

2445:    Level: intermediate

2447: .seealso: VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2448: @*/
2449: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayWrite(Vec v, PetscScalar **a)
2450: {

2455:  #if defined(PETSC_HAVE_CUDA)
2456:   v->offloadmask = PETSC_OFFLOAD_GPU;
2457:   if (a) *a = NULL;
2458:  #endif
2459:   PetscObjectStateIncrease((PetscObject)v);
2460:   return(0);
2461: }

2463: /*@C
2464:    VecCUDAPlaceArray - Allows one to replace the GPU array in a vector with a
2465:    GPU array provided by the user. This is useful to avoid copying an
2466:    array into a vector.

2468:    Not Collective

2470:    Input Parameters:
2471: +  vec - the vector
2472: -  array - the GPU array

2474:    Notes:
2475:    You can return to the original GPU array with a call to VecCUDAResetArray()
2476:    It is not possible to use VecCUDAPlaceArray() and VecPlaceArray() at the
2477:    same time on the same vector.

2479:    Level: developer

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

2483: @*/
2484: PetscErrorCode VecCUDAPlaceArray(Vec vin,const PetscScalar a[])
2485: {

2490: #if defined(PETSC_HAVE_CUDA)
2491:   VecCUDACopyToGPU(vin);
2492:   if (((Vec_Seq*)vin->data)->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"VecCUDAPlaceArray()/VecPlaceArray() was already called on this vector, without a call to VecCUDAResetArray()/VecResetArray()");
2493:   ((Vec_Seq*)vin->data)->unplacedarray  = (PetscScalar *) ((Vec_CUDA*)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2494:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2495:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2496: #endif
2497:   PetscObjectStateIncrease((PetscObject)vin);
2498:   return(0);
2499: }

2501: /*@C
2502:    VecCUDAReplaceArray - Allows one to replace the GPU array in a vector
2503:    with a GPU array provided by the user. This is useful to avoid copying
2504:    a GPU array into a vector.

2506:    Not Collective

2508:    Input Parameters:
2509: +  vec - the vector
2510: -  array - the GPU array

2512:    Notes:
2513:    This permanently replaces the GPU array and frees the memory associated
2514:    with the old GPU array.

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

2519:    Not supported from Fortran

2521:    Level: developer

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

2525: @*/
2526: PetscErrorCode VecCUDAReplaceArray(Vec vin,const PetscScalar a[])
2527: {
2528: #if defined(PETSC_HAVE_CUDA)
2529:   cudaError_t err;
2530: #endif

2535: #if defined(PETSC_HAVE_CUDA)
2536:   if (((Vec_CUDA*)vin->spptr)->GPUarray_allocated) {
2537:     err = cudaFree(((Vec_CUDA*)vin->spptr)->GPUarray_allocated);CHKERRCUDA(err);
2538:   }
2539:   ((Vec_CUDA*)vin->spptr)->GPUarray_allocated = ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2540:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2541: #endif
2542:   PetscObjectStateIncrease((PetscObject)vin);
2543:   return(0);
2544: }

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

2550:    Not Collective

2552:    Input Parameters:
2553: .  vec - the vector

2555:    Level: developer

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

2559: @*/
2560: PetscErrorCode VecCUDAResetArray(Vec vin)
2561: {

2566: #if defined(PETSC_HAVE_CUDA)
2567:   VecCUDACopyToGPU(vin);
2568:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar *) ((Vec_Seq*)vin->data)->unplacedarray;
2569:   ((Vec_Seq*)vin->data)->unplacedarray = 0;
2570:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2571: #endif
2572:   PetscObjectStateIncrease((PetscObject)vin);
2573:   return(0);
2574: }

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

2579:    This function has semantics similar to VecGetArray():  the pointer
2580:    returned by this function points to a consistent view of the vector
2581:    data.  This may involve a copy operation of data from the host to the
2582:    device if the data on the device is out of date.  If the device
2583:    memory hasn't been allocated previously it will be allocated as part
2584:    of this function call.  VecHIPGetArray() assumes that
2585:    the user will modify the vector data.  This is similar to
2586:    intent(inout) in fortran.

2588:    The HIP device pointer has to be released by calling
2589:    VecHIPRestoreArray().  Upon restoring the vector data
2590:    the data on the host will be marked as out of date.  A subsequent
2591:    access of the host data will thus incur a data transfer from the
2592:    device to the host.

2594:    Input Parameter:
2595: .  v - the vector

2597:    Output Parameter:
2598: .  a - the HIP device pointer

2600:    Fortran note:
2601:    This function is not currently available from Fortran.

2603:    Level: intermediate

2605: .seealso: VecHIPRestoreArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2606: @*/
2607: PETSC_EXTERN PetscErrorCode VecHIPGetArray(Vec v, PetscScalar **a)
2608: {
2609: #if defined(PETSC_HAVE_HIP)
2611: #endif

2615: #if defined(PETSC_HAVE_HIP)
2616:   *a   = 0;
2617:   VecHIPCopyToGPU(v);
2618:   *a   = ((Vec_HIP*)v->spptr)->GPUarray;
2619: #endif
2620:   return(0);
2621: }

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

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

2630:    Input Parameter:
2631: +  v - the vector
2632: -  a - the HIP device pointer.  This pointer is invalid after
2633:        VecHIPRestoreArray() returns.

2635:    Fortran note:
2636:    This function is not currently available from Fortran.

2638:    Level: intermediate

2640: .seealso: VecHIPGetArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2641: @*/
2642: PETSC_EXTERN PetscErrorCode VecHIPRestoreArray(Vec v, PetscScalar **a)
2643: {

2648: #if defined(PETSC_HAVE_HIP)
2649:   v->offloadmask = PETSC_OFFLOAD_GPU;
2650: #endif

2652:   PetscObjectStateIncrease((PetscObject)v);
2653:   return(0);
2654: }

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

2659:    This function is analogous to VecGetArrayRead():  The pointer
2660:    returned by this function points to a consistent view of the vector
2661:    data.  This may involve a copy operation of data from the host to the
2662:    device if the data on the device is out of date.  If the device
2663:    memory hasn't been allocated previously it will be allocated as part
2664:    of this function call.  VecHIPGetArrayRead() assumes that the
2665:    user will not modify the vector data.  This is analgogous to
2666:    intent(in) in Fortran.

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

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

2677:    Output Parameter:
2678: .  a - the HIP pointer.

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

2683:    Level: intermediate

2685: .seealso: VecHIPRestoreArrayRead(), VecHIPGetArray(), VecHIPGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2686: @*/
2687: PETSC_EXTERN PetscErrorCode VecHIPGetArrayRead(Vec v, const PetscScalar **a)
2688: {
2689: #if defined(PETSC_HAVE_HIP)
2691: #endif

2695: #if defined(PETSC_HAVE_HIP)
2696:   *a   = 0;
2697:   VecHIPCopyToGPU(v);
2698:   *a   = ((Vec_HIP*)v->spptr)->GPUarray;
2699: #endif
2700:   return(0);
2701: }

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

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

2711:    Input Parameter:
2712: +  v - the vector
2713: -  a - the HIP device pointer.  This pointer is invalid after
2714:        VecHIPRestoreArrayRead() returns.

2716:    Fortran note:
2717:    This function is not currently available from Fortran.

2719:    Level: intermediate

2721: .seealso: VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecHIPGetArray(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2722: @*/
2723: PETSC_EXTERN PetscErrorCode VecHIPRestoreArrayRead(Vec v, const PetscScalar **a)
2724: {
2727:   *a = NULL;
2728:   return(0);
2729: }

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

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

2740:    The device pointer needs to be released with
2741:    VecHIPRestoreArrayWrite().  When the pointer is released the
2742:    host data of the vector is marked as out of data.  Subsequent access
2743:    of the host data with e.g. VecGetArray() incurs a device to host data
2744:    transfer.

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

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

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

2755:    Level: advanced

2757: .seealso: VecHIPRestoreArrayWrite(), VecHIPGetArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2758: @*/
2759: PETSC_EXTERN PetscErrorCode VecHIPGetArrayWrite(Vec v, PetscScalar **a)
2760: {
2761: #if defined(PETSC_HAVE_HIP)
2763: #endif

2767: #if defined(PETSC_HAVE_HIP)
2768:   *a   = 0;
2769:   VecHIPAllocateCheck(v);
2770:   *a   = ((Vec_HIP*)v->spptr)->GPUarray;
2771: #endif
2772:   return(0);
2773: }

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

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

2782:    Input Parameter:
2783: +  v - the vector
2784: -  a - the HIP device pointer.  This pointer is invalid after
2785:        VecHIPRestoreArrayWrite() returns.

2787:    Fortran note:
2788:    This function is not currently available from Fortran.

2790:    Level: intermediate

2792: .seealso: VecHIPGetArrayWrite(), VecHIPGetArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2793: @*/
2794: PETSC_EXTERN PetscErrorCode VecHIPRestoreArrayWrite(Vec v, PetscScalar **a)
2795: {

2800: #if defined(PETSC_HAVE_HIP)
2801:   v->offloadmask = PETSC_OFFLOAD_GPU;
2802: #endif

2804:   PetscObjectStateIncrease((PetscObject)v);
2805:   return(0);
2806: }

2808: /*@C
2809:    VecHIPPlaceArray - Allows one to replace the GPU array in a vector with a
2810:    GPU array provided by the user. This is useful to avoid copying an
2811:    array into a vector.

2813:    Not Collective

2815:    Input Parameters:
2816: +  vec - the vector
2817: -  array - the GPU array

2819:    Notes:
2820:    You can return to the original GPU array with a call to VecHIPResetArray()
2821:    It is not possible to use VecHIPPlaceArray() and VecPlaceArray() at the
2822:    same time on the same vector.

2824:    Level: developer

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

2828: @*/
2829: PetscErrorCode VecHIPPlaceArray(Vec vin,const PetscScalar a[])
2830: {

2835: #if defined(PETSC_HAVE_HIP)
2836:   VecHIPCopyToGPU(vin);
2837:   if (((Vec_Seq*)vin->data)->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"VecHIPPlaceArray()/VecPlaceArray() was already called on this vector, without a call to VecHIPResetArray()/VecResetArray()");
2838:   ((Vec_Seq*)vin->data)->unplacedarray  = (PetscScalar *) ((Vec_HIP*)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2839:   ((Vec_HIP*)vin->spptr)->GPUarray = (PetscScalar*)a;
2840:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2841: #endif
2842:   PetscObjectStateIncrease((PetscObject)vin);
2843:   return(0);
2844: }

2846: /*@C
2847:    VecHIPReplaceArray - Allows one to replace the GPU array in a vector
2848:    with a GPU array provided by the user. This is useful to avoid copying
2849:    a GPU array into a vector.

2851:    Not Collective

2853:    Input Parameters:
2854: +  vec - the vector
2855: -  array - the GPU array

2857:    Notes:
2858:    This permanently replaces the GPU array and frees the memory associated
2859:    with the old GPU array.

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

2864:    Not supported from Fortran

2866:    Level: developer

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

2870: @*/
2871: PetscErrorCode VecHIPReplaceArray(Vec vin,const PetscScalar a[])
2872: {
2873: #if defined(PETSC_HAVE_HIP)
2874:   hipError_t err;
2875: #endif

2880: #if defined(PETSC_HAVE_HIP)
2881:   err = hipFree(((Vec_HIP*)vin->spptr)->GPUarray);CHKERRHIP(err);
2882:   ((Vec_HIP*)vin->spptr)->GPUarray = (PetscScalar*)a;
2883:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2884: #endif
2885:   PetscObjectStateIncrease((PetscObject)vin);
2886:   return(0);
2887: }

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

2893:    Not Collective

2895:    Input Parameters:
2896: .  vec - the vector

2898:    Level: developer

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

2902: @*/
2903: PetscErrorCode VecHIPResetArray(Vec vin)
2904: {

2909: #if defined(PETSC_HAVE_HIP)
2910:   VecHIPCopyToGPU(vin);
2911:   ((Vec_HIP*)vin->spptr)->GPUarray = (PetscScalar *) ((Vec_Seq*)vin->data)->unplacedarray;
2912:   ((Vec_Seq*)vin->data)->unplacedarray = 0;
2913:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2914: #endif
2915:   PetscObjectStateIncrease((PetscObject)vin);
2916:   return(0);
2917: }

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

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

2926:     Collective on Vec

2928:     Input Parameters:
2929: +   x - a vector to mimic
2930: -   n - the number of vectors to obtain

2932:     Output Parameters:
2933: +   y - Fortran90 pointer to the array of vectors
2934: -   ierr - error code

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

2941:     Vec x
2942:     Vec, pointer :: y(:)
2943:     ....
2944:     call VecDuplicateVecsF90(x,2,y,ierr)
2945:     call VecSet(y(2),alpha,ierr)
2946:     call VecSet(y(2),alpha,ierr)
2947:     ....
2948:     call VecDestroyVecsF90(2,y,ierr)
2949: .ve

2951:     Notes:
2952:     Not yet supported for all F90 compilers

2954:     Use VecDestroyVecsF90() to free the space.

2956:     Level: beginner

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

2960: M*/

2962: /*MC
2963:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2964:     VecGetArrayF90().

2966:     Synopsis:
2967:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2969:     Logically Collective on Vec

2971:     Input Parameters:
2972: +   x - vector
2973: -   xx_v - the Fortran90 pointer to the array

2975:     Output Parameter:
2976: .   ierr - error code

2978:     Example of Usage:
2979: .vb
2980: #include <petsc/finclude/petscvec.h>
2981:     use petscvec

2983:     PetscScalar, pointer :: xx_v(:)
2984:     ....
2985:     call VecGetArrayF90(x,xx_v,ierr)
2986:     xx_v(3) = a
2987:     call VecRestoreArrayF90(x,xx_v,ierr)
2988: .ve

2990:     Level: beginner

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

2994: M*/

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

2999:     Synopsis:
3000:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

3002:     Collective on Vec

3004:     Input Parameters:
3005: +   n - the number of vectors previously obtained
3006: -   x - pointer to array of vector pointers

3008:     Output Parameter:
3009: .   ierr - error code

3011:     Notes:
3012:     Not yet supported for all F90 compilers

3014:     Level: beginner

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

3018: M*/

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

3026:     Synopsis:
3027:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

3029:     Logically Collective on Vec

3031:     Input Parameter:
3032: .   x - vector

3034:     Output Parameters:
3035: +   xx_v - the Fortran90 pointer to the array
3036: -   ierr - error code

3038:     Example of Usage:
3039: .vb
3040: #include <petsc/finclude/petscvec.h>
3041:     use petscvec

3043:     PetscScalar, pointer :: xx_v(:)
3044:     ....
3045:     call VecGetArrayF90(x,xx_v,ierr)
3046:     xx_v(3) = a
3047:     call VecRestoreArrayF90(x,xx_v,ierr)
3048: .ve

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

3052:     Level: beginner

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

3056: M*/

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

3064:     Synopsis:
3065:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

3067:     Logically Collective on Vec

3069:     Input Parameter:
3070: .   x - vector

3072:     Output Parameters:
3073: +   xx_v - the Fortran90 pointer to the array
3074: -   ierr - error code

3076:     Example of Usage:
3077: .vb
3078: #include <petsc/finclude/petscvec.h>
3079:     use petscvec

3081:     PetscScalar, pointer :: xx_v(:)
3082:     ....
3083:     call VecGetArrayReadF90(x,xx_v,ierr)
3084:     a = xx_v(3)
3085:     call VecRestoreArrayReadF90(x,xx_v,ierr)
3086: .ve

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

3090:     Level: beginner

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

3094: M*/

3096: /*MC
3097:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
3098:     VecGetArrayReadF90().

3100:     Synopsis:
3101:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

3103:     Logically Collective on Vec

3105:     Input Parameters:
3106: +   x - vector
3107: -   xx_v - the Fortran90 pointer to the array

3109:     Output Parameter:
3110: .   ierr - error code

3112:     Example of Usage:
3113: .vb
3114: #include <petsc/finclude/petscvec.h>
3115:     use petscvec

3117:     PetscScalar, pointer :: xx_v(:)
3118:     ....
3119:     call VecGetArrayReadF90(x,xx_v,ierr)
3120:     a = xx_v(3)
3121:     call VecRestoreArrayReadF90(x,xx_v,ierr)
3122: .ve

3124:     Level: beginner

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

3128: M*/

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

3135:    Logically Collective

3137:    Input Parameter:
3138: +  x - the vector
3139: .  m - first dimension of two dimensional array
3140: .  n - second dimension of two dimensional array
3141: .  mstart - first index you will use in first coordinate direction (often 0)
3142: -  nstart - first index in the second coordinate direction (often 0)

3144:    Output Parameter:
3145: .  a - location to put pointer to the array

3147:    Level: developer

3149:   Notes:
3150:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
3151:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3152:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3153:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

3157: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3158:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3159:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3160: @*/
3161: PetscErrorCode  VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3162: {
3164:   PetscInt       i,N;
3165:   PetscScalar    *aa;

3171:   VecGetLocalSize(x,&N);
3172:   if (m*n != N) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 2d array dimensions %D by %D",N,m,n);
3173:   VecGetArray(x,&aa);

3175:   PetscMalloc1(m,a);
3176:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
3177:   *a -= mstart;
3178:   return(0);
3179: }

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

3186:    Logically Collective

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

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

3198:    Level: developer

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

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

3208:    Concepts: vector^accessing local values as 2d array

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

3224:   VecGetLocalSize(x,&N);
3225:   if (m*n != N) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 2d array dimensions %D by %D",N,m,n);
3226:   VecGetArrayWrite(x,&aa);

3228:   PetscMalloc1(m,a);
3229:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
3230:   *a -= mstart;
3231:   return(0);
3232: }

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

3237:    Logically Collective

3239:    Input Parameters:
3240: +  x - the vector
3241: .  m - first dimension of two dimensional array
3242: .  n - second dimension of the two dimensional array
3243: .  mstart - first index you will use in first coordinate direction (often 0)
3244: .  nstart - first index in the second coordinate direction (often 0)
3245: -  a - location of pointer to array obtained from VecGetArray2d()

3247:    Level: developer

3249:    Notes:
3250:    For regular PETSc vectors this routine does not involve any copies. For
3251:    any special vectors that do not store local vector data in a contiguous
3252:    array, this routine will copy the data back into the underlying
3253:    vector data structure from the array obtained with VecGetArray().

3255:    This routine actually zeros out the a pointer.

3257: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3258:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3259:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3260: @*/
3261: PetscErrorCode  VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3262: {
3264:   void           *dummy;

3270:   dummy = (void*)(*a + mstart);
3271:   PetscFree(dummy);
3272:   VecRestoreArray(x,NULL);
3273:   return(0);
3274: }

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

3279:    Logically Collective

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

3289:    Level: developer

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

3297:    This routine actually zeros out the a pointer.

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

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

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

3323:    Logically Collective

3325:    Input Parameter:
3326: +  x - the vector
3327: .  m - first dimension of two dimensional array
3328: -  mstart - first index you will use in first coordinate direction (often 0)

3330:    Output Parameter:
3331: .  a - location to put pointer to the array

3333:    Level: developer

3335:   Notes:
3336:    For a vector obtained from DMCreateLocalVector() mstart are likely
3337:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3338:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3342: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3343:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3344:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3345: @*/
3346: PetscErrorCode  VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3347: {
3349:   PetscInt       N;

3355:   VecGetLocalSize(x,&N);
3356:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3357:   VecGetArray(x,a);
3358:   *a  -= mstart;
3359:   return(0);
3360: }

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

3367:    Logically Collective

3369:    Input Parameter:
3370: +  x - the vector
3371: .  m - first dimension of two dimensional array
3372: -  mstart - first index you will use in first coordinate direction (often 0)

3374:    Output Parameter:
3375: .  a - location to put pointer to the array

3377:    Level: developer

3379:   Notes:
3380:    For a vector obtained from DMCreateLocalVector() mstart are likely
3381:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3382:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3386: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3387:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3388:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3389: @*/
3390: PetscErrorCode  VecGetArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3391: {
3393:   PetscInt       N;

3399:   VecGetLocalSize(x,&N);
3400:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3401:   VecGetArrayWrite(x,a);
3402:   *a  -= mstart;
3403:   return(0);
3404: }

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

3409:    Logically Collective

3411:    Input Parameters:
3412: +  x - the vector
3413: .  m - first dimension of two dimensional array
3414: .  mstart - first index you will use in first coordinate direction (often 0)
3415: -  a - location of pointer to array obtained from VecGetArray21()

3417:    Level: developer

3419:    Notes:
3420:    For regular PETSc vectors this routine does not involve any copies. For
3421:    any special vectors that do not store local vector data in a contiguous
3422:    array, this routine will copy the data back into the underlying
3423:    vector data structure from the array obtained with VecGetArray1d().

3425:    This routine actually zeros out the a pointer.

3427: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3428:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3429:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3430: @*/
3431: PetscErrorCode  VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3432: {

3438:   VecRestoreArray(x,NULL);
3439:   return(0);
3440: }

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

3445:    Logically Collective

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

3453:    Level: developer

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

3461:    This routine actually zeros out the a pointer.

3463:    Concepts: vector^accessing local values as 1d array

3465: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3466:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3467:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3468: @*/
3469: PetscErrorCode  VecRestoreArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3470: {

3476:   VecRestoreArrayWrite(x,NULL);
3477:   return(0);
3478: }

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

3485:    Logically Collective

3487:    Input Parameter:
3488: +  x - the vector
3489: .  m - first dimension of three dimensional array
3490: .  n - second dimension of three dimensional array
3491: .  p - third dimension of three dimensional array
3492: .  mstart - first index you will use in first coordinate direction (often 0)
3493: .  nstart - first index in the second coordinate direction (often 0)
3494: -  pstart - first index in the third coordinate direction (often 0)

3496:    Output Parameter:
3497: .  a - location to put pointer to the array

3499:    Level: developer

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

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

3509: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3510:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3511:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3512: @*/
3513: PetscErrorCode  VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3514: {
3516:   PetscInt       i,N,j;
3517:   PetscScalar    *aa,**b;

3523:   VecGetLocalSize(x,&N);
3524:   if (m*n*p != N) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 3d array dimensions %D by %D by %D",N,m,n,p);
3525:   VecGetArray(x,&aa);

3527:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3528:   b    = (PetscScalar**)((*a) + m);
3529:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3530:   for (i=0; i<m; i++)
3531:     for (j=0; j<n; j++)
3532:       b[i*n+j] = aa + i*n*p + j*p - pstart;

3534:   *a -= mstart;
3535:   return(0);
3536: }

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

3543:    Logically Collective

3545:    Input Parameter:
3546: +  x - the vector
3547: .  m - first dimension of three dimensional array
3548: .  n - second dimension of three dimensional array
3549: .  p - third dimension of three dimensional array
3550: .  mstart - first index you will use in first coordinate direction (often 0)
3551: .  nstart - first index in the second coordinate direction (often 0)
3552: -  pstart - first index in the third coordinate direction (often 0)

3554:    Output Parameter:
3555: .  a - location to put pointer to the array

3557:    Level: developer

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

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

3567:    Concepts: vector^accessing local values as 3d array

3569: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3570:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3571:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3572: @*/
3573: PetscErrorCode  VecGetArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3574: {
3576:   PetscInt       i,N,j;
3577:   PetscScalar    *aa,**b;

3583:   VecGetLocalSize(x,&N);
3584:   if (m*n*p != N) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 3d array dimensions %D by %D by %D",N,m,n,p);
3585:   VecGetArrayWrite(x,&aa);

3587:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3588:   b    = (PetscScalar**)((*a) + m);
3589:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3590:   for (i=0; i<m; i++)
3591:     for (j=0; j<n; j++)
3592:       b[i*n+j] = aa + i*n*p + j*p - pstart;

3594:   *a -= mstart;
3595:   return(0);
3596: }

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

3601:    Logically Collective

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

3613:    Level: developer

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

3621:    This routine actually zeros out the a pointer.

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

3636:   dummy = (void*)(*a + mstart);
3637:   PetscFree(dummy);
3638:   VecRestoreArray(x,NULL);
3639:   return(0);
3640: }

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

3645:    Logically Collective

3647:    Input Parameters:
3648: +  x - the vector
3649: .  m - first dimension of three dimensional array
3650: .  n - second dimension of the three dimensional array
3651: .  p - third dimension of the three dimensional array
3652: .  mstart - first index you will use in first coordinate direction (often 0)
3653: .  nstart - first index in the second coordinate direction (often 0)
3654: .  pstart - first index in the third coordinate direction (often 0)
3655: -  a - location of pointer to array obtained from VecGetArray3d()

3657:    Level: developer

3659:    Notes:
3660:    For regular PETSc vectors this routine does not involve any copies. For
3661:    any special vectors that do not store local vector data in a contiguous
3662:    array, this routine will copy the data back into the underlying
3663:    vector data structure from the array obtained with VecGetArray().

3665:    This routine actually zeros out the a pointer.

3667: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3668:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3669:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3670: @*/
3671: PetscErrorCode  VecRestoreArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3672: {
3674:   void           *dummy;

3680:   dummy = (void*)(*a + mstart);
3681:   PetscFree(dummy);
3682:   VecRestoreArrayWrite(x,NULL);
3683:   return(0);
3684: }

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

3691:    Logically Collective

3693:    Input Parameter:
3694: +  x - the vector
3695: .  m - first dimension of four dimensional array
3696: .  n - second dimension of four dimensional array
3697: .  p - third dimension of four dimensional array
3698: .  q - fourth dimension of four dimensional array
3699: .  mstart - first index you will use in first coordinate direction (often 0)
3700: .  nstart - first index in the second coordinate direction (often 0)
3701: .  pstart - first index in the third coordinate direction (often 0)
3702: -  qstart - first index in the fourth coordinate direction (often 0)

3704:    Output Parameter:
3705: .  a - location to put pointer to the array

3707:    Level: beginner

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

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

3717: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3718:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3719:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3720: @*/
3721: PetscErrorCode  VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3722: {
3724:   PetscInt       i,N,j,k;
3725:   PetscScalar    *aa,***b,**c;

3731:   VecGetLocalSize(x,&N);
3732:   if (m*n*p*q != N) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 4d array dimensions %D by %D by %D by %D",N,m,n,p,q);
3733:   VecGetArray(x,&aa);

3735:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3736:   b    = (PetscScalar***)((*a) + m);
3737:   c    = (PetscScalar**)(b + m*n);
3738:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3739:   for (i=0; i<m; i++)
3740:     for (j=0; j<n; j++)
3741:       b[i*n+j] = c + i*n*p + j*p - pstart;
3742:   for (i=0; i<m; i++)
3743:     for (j=0; j<n; j++)
3744:       for (k=0; k<p; k++)
3745:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3746:   *a -= mstart;
3747:   return(0);
3748: }

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

3755:    Logically Collective

3757:    Input Parameter:
3758: +  x - the vector
3759: .  m - first dimension of four dimensional array
3760: .  n - second dimension of four dimensional array
3761: .  p - third dimension of four dimensional array
3762: .  q - fourth dimension of four dimensional array
3763: .  mstart - first index you will use in first coordinate direction (often 0)
3764: .  nstart - first index in the second coordinate direction (often 0)
3765: .  pstart - first index in the third coordinate direction (often 0)
3766: -  qstart - first index in the fourth coordinate direction (often 0)

3768:    Output Parameter:
3769: .  a - location to put pointer to the array

3771:    Level: beginner

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

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

3781:    Concepts: vector^accessing local values as 3d array

3783: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3784:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3785:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3786: @*/
3787: PetscErrorCode  VecGetArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3788: {
3790:   PetscInt       i,N,j,k;
3791:   PetscScalar    *aa,***b,**c;

3797:   VecGetLocalSize(x,&N);
3798:   if (m*n*p*q != N) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 4d array dimensions %D by %D by %D by %D",N,m,n,p,q);
3799:   VecGetArrayWrite(x,&aa);

3801:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3802:   b    = (PetscScalar***)((*a) + m);
3803:   c    = (PetscScalar**)(b + m*n);
3804:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3805:   for (i=0; i<m; i++)
3806:     for (j=0; j<n; j++)
3807:       b[i*n+j] = c + i*n*p + j*p - pstart;
3808:   for (i=0; i<m; i++)
3809:     for (j=0; j<n; j++)
3810:       for (k=0; k<p; k++)
3811:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3812:   *a -= mstart;
3813:   return(0);
3814: }

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

3819:    Logically Collective

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

3833:    Level: beginner

3835:    Notes:
3836:    For regular PETSc vectors this routine does not involve any copies. For
3837:    any special vectors that do not store local vector data in a contiguous
3838:    array, this routine will copy the data back into the underlying
3839:    vector data structure from the array obtained with VecGetArray().

3841:    This routine actually zeros out the a pointer.

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

3856:   dummy = (void*)(*a + mstart);
3857:   PetscFree(dummy);
3858:   VecRestoreArray(x,NULL);
3859:   return(0);
3860: }

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

3865:    Logically Collective

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

3879:    Level: beginner

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

3887:    This routine actually zeros out the a pointer.

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

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

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

3913:    Logically Collective

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

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

3925:    Level: developer

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

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:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3938: @*/
3939: PetscErrorCode  VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3940: {
3941:   PetscErrorCode    ierr;
3942:   PetscInt          i,N;
3943:   const PetscScalar *aa;

3949:   VecGetLocalSize(x,&N);
3950:   if (m*n != N) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 2d array dimensions %D by %D",N,m,n);
3951:   VecGetArrayRead(x,&aa);

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

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

3962:    Logically Collective

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

3972:    Level: developer

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

3980:    This routine actually zeros out the a pointer.

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

3995:   dummy = (void*)(*a + mstart);
3996:   PetscFree(dummy);
3997:   VecRestoreArrayRead(x,NULL);
3998:   return(0);
3999: }

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

4006:    Logically Collective

4008:    Input Parameter:
4009: +  x - the vector
4010: .  m - first dimension of two dimensional array
4011: -  mstart - first index you will use in first coordinate direction (often 0)

4013:    Output Parameter:
4014: .  a - location to put pointer to the array

4016:    Level: developer

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

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

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

4038:   VecGetLocalSize(x,&N);
4039:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
4040:   VecGetArrayRead(x,(const PetscScalar**)a);
4041:   *a  -= mstart;
4042:   return(0);
4043: }

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

4048:    Logically Collective

4050:    Input Parameters:
4051: +  x - the vector
4052: .  m - first dimension of two dimensional array
4053: .  mstart - first index you will use in first coordinate direction (often 0)
4054: -  a - location of pointer to array obtained from VecGetArray21()

4056:    Level: developer

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

4064:    This routine actually zeros out the a pointer.

4066: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
4067:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
4068:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
4069: @*/
4070: PetscErrorCode  VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
4071: {

4077:   VecRestoreArrayRead(x,NULL);
4078:   return(0);
4079: }

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

4086:    Logically Collective

4088:    Input Parameter:
4089: +  x - the vector
4090: .  m - first dimension of three dimensional array
4091: .  n - second dimension of three dimensional array
4092: .  p - third dimension of three dimensional array
4093: .  mstart - first index you will use in first coordinate direction (often 0)
4094: .  nstart - first index in the second coordinate direction (often 0)
4095: -  pstart - first index in the third coordinate direction (often 0)

4097:    Output Parameter:
4098: .  a - location to put pointer to the array

4100:    Level: developer

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

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

4110: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
4111:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
4112:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
4113: @*/
4114: PetscErrorCode  VecGetArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
4115: {
4116:   PetscErrorCode    ierr;
4117:   PetscInt          i,N,j;
4118:   const PetscScalar *aa;
4119:   PetscScalar       **b;

4125:   VecGetLocalSize(x,&N);
4126:   if (m*n*p != N) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 3d array dimensions %D by %D by %D",N,m,n,p);
4127:   VecGetArrayRead(x,&aa);

4129:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
4130:   b    = (PetscScalar**)((*a) + m);
4131:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
4132:   for (i=0; i<m; i++)
4133:     for (j=0; j<n; j++)
4134:       b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;

4136:   *a -= mstart;
4137:   return(0);
4138: }

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

4143:    Logically Collective

4145:    Input Parameters:
4146: +  x - the vector
4147: .  m - first dimension of three dimensional array
4148: .  n - second dimension of the three dimensional array
4149: .  p - third dimension of the three dimensional array
4150: .  mstart - first index you will use in first coordinate direction (often 0)
4151: .  nstart - first index in the second coordinate direction (often 0)
4152: .  pstart - first index in the third coordinate direction (often 0)
4153: -  a - location of pointer to array obtained from VecGetArray3dRead()

4155:    Level: developer

4157:    Notes:
4158:    For regular PETSc vectors this routine does not involve any copies. For
4159:    any special vectors that do not store local vector data in a contiguous
4160:    array, this routine will copy the data back into the underlying
4161:    vector data structure from the array obtained with VecGetArray().

4163:    This routine actually zeros out the a pointer.

4165: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
4166:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
4167:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
4168: @*/
4169: PetscErrorCode  VecRestoreArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
4170: {
4172:   void           *dummy;

4178:   dummy = (void*)(*a + mstart);
4179:   PetscFree(dummy);
4180:   VecRestoreArrayRead(x,NULL);
4181:   return(0);
4182: }

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

4189:    Logically Collective

4191:    Input Parameter:
4192: +  x - the vector
4193: .  m - first dimension of four dimensional array
4194: .  n - second dimension of four dimensional array
4195: .  p - third dimension of four dimensional array
4196: .  q - fourth dimension of four dimensional array
4197: .  mstart - first index you will use in first coordinate direction (often 0)
4198: .  nstart - first index in the second coordinate direction (often 0)
4199: .  pstart - first index in the third coordinate direction (often 0)
4200: -  qstart - first index in the fourth coordinate direction (often 0)

4202:    Output Parameter:
4203: .  a - location to put pointer to the array

4205:    Level: beginner

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

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

4215: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
4216:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
4217:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
4218: @*/
4219: PetscErrorCode  VecGetArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
4220: {
4221:   PetscErrorCode    ierr;
4222:   PetscInt          i,N,j,k;
4223:   const PetscScalar *aa;
4224:   PetscScalar       ***b,**c;

4230:   VecGetLocalSize(x,&N);
4231:   if (m*n*p*q != N) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 4d array dimensions %D by %D by %D by %D",N,m,n,p,q);
4232:   VecGetArrayRead(x,&aa);

4234:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
4235:   b    = (PetscScalar***)((*a) + m);
4236:   c    = (PetscScalar**)(b + m*n);
4237:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
4238:   for (i=0; i<m; i++)
4239:     for (j=0; j<n; j++)
4240:       b[i*n+j] = c + i*n*p + j*p - pstart;
4241:   for (i=0; i<m; i++)
4242:     for (j=0; j<n; j++)
4243:       for (k=0; k<p; k++)
4244:         c[i*n*p+j*p+k] = (PetscScalar*) aa + i*n*p*q + j*p*q + k*q - qstart;
4245:   *a -= mstart;
4246:   return(0);
4247: }

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

4252:    Logically Collective

4254:    Input Parameters:
4255: +  x - the vector
4256: .  m - first dimension of four dimensional array
4257: .  n - second dimension of the four dimensional array
4258: .  p - third dimension of the four dimensional array
4259: .  q - fourth dimension of the four dimensional array
4260: .  mstart - first index you will use in first coordinate direction (often 0)
4261: .  nstart - first index in the second coordinate direction (often 0)
4262: .  pstart - first index in the third coordinate direction (often 0)
4263: .  qstart - first index in the fourth coordinate direction (often 0)
4264: -  a - location of pointer to array obtained from VecGetArray4dRead()

4266:    Level: beginner

4268:    Notes:
4269:    For regular PETSc vectors this routine does not involve any copies. For
4270:    any special vectors that do not store local vector data in a contiguous
4271:    array, this routine will copy the data back into the underlying
4272:    vector data structure from the array obtained with VecGetArray().

4274:    This routine actually zeros out the a pointer.

4276: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
4277:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
4278:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
4279: @*/
4280: PetscErrorCode  VecRestoreArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
4281: {
4283:   void           *dummy;

4289:   dummy = (void*)(*a + mstart);
4290:   PetscFree(dummy);
4291:   VecRestoreArrayRead(x,NULL);
4292:   return(0);
4293: }

4295: #if defined(PETSC_USE_DEBUG)

4297: /*@
4298:    VecLockGet  - Gets the current lock status of a vector

4300:    Logically Collective on Vec

4302:    Input Parameter:
4303: .  x - the vector

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

4309:    Level: beginner

4311: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop()
4312: @*/
4313: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
4314: {
4317:   *state = x->lock;
4318:   return(0);
4319: }

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

4324:    Logically Collective on Vec

4326:    Input Parameter:
4327: .  x - the vector

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

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

4335:    Level: beginner

4337: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPop(), VecLockGet()
4338: @*/
4339: PetscErrorCode VecLockReadPush(Vec x)
4340: {
4343:   if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector is already locked for exclusive write access but you want to read it");
4344:   x->lock++;
4345:   return(0);
4346: }

4348: /*@
4349:    VecLockReadPop  - Pops a read-only lock from a vector

4351:    Logically Collective on Vec

4353:    Input Parameter:
4354: .  x - the vector

4356:    Level: beginner

4358: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockGet()
4359: @*/
4360: PetscErrorCode VecLockReadPop(Vec x)
4361: {
4364:   x->lock--;
4365:   if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector has been unlocked from read-only access too many times");
4366:   return(0);
4367: }

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

4372:    Logically Collective on Vec

4374:    Input Parameter:
4375: +  x   - the vector
4376: -  flg - PETSC_TRUE to lock the vector for writing; PETSC_FALSE to unlock it.

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

4384:        VecGetArray(x,&xdata); // begin phase
4385:        VecLockWriteSet_Private(v,PETSC_TRUE);

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

4389:        VecRestoreArray(x,&vdata); // end phase
4390:        VecLockWriteSet_Private(v,PETSC_FALSE);

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

4395:    Level: beginner

4397: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop(), VecLockGet()
4398: @*/
4399: PetscErrorCode VecLockWriteSet_Private(Vec x,PetscBool flg)
4400: {
4403:   if (flg) {
4404:     if (x->lock > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector is already locked for read-only access but you want to write it");
4405:     else if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector is already locked for exclusive write access but you want to write it");
4406:     else x->lock = -1;
4407:   } else {
4408:     if (x->lock != -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector is not locked for exclusive write access but you want to unlock it from that");
4409:     x->lock = 0;
4410:   }
4411:   return(0);
4412: }

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

4417:    Level: deprecated

4419: .seealso: VecLockReadPush()
4420: @*/
4421: PetscErrorCode VecLockPush(Vec x)
4422: {
4425:   VecLockReadPush(x);
4426:   return(0);
4427: }

4429: /*@
4430:    VecLockPop  - Pops a read-only lock from a vector

4432:    Level: deprecated

4434: .seealso: VecLockReadPop()
4435: @*/
4436: PetscErrorCode VecLockPop(Vec x)
4437: {
4440:   VecLockReadPop(x);
4441:   return(0);
4442: }

4444: #endif