Actual source code: rvector.c

petsc-master 2020-11-28
Report Typos and Errors
  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: static PetscInt VecGetSubVectorSavedStateId = -1;

 14: PETSC_EXTERN PetscErrorCode VecValidValues(Vec vec,PetscInt argnum,PetscBool begin)
 15: {
 16: #if defined(PETSC_USE_DEBUG)
 17:   PetscErrorCode    ierr;
 18:   PetscInt          n,i;
 19:   const PetscScalar *x;

 22: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL)
 23:   if ((vec->petscnative || vec->ops->getarray) && (vec->offloadmask & PETSC_OFFLOAD_CPU)) {
 24: #else
 25:   if (vec->petscnative || vec->ops->getarray) {
 26: #endif
 27:     VecGetLocalSize(vec,&n);
 28:     VecGetArrayRead(vec,&x);
 29:     for (i=0; i<n; i++) {
 30:       if (begin) {
 31:         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);
 32:       } else {
 33:         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);
 34:       }
 35:     }
 36:     VecRestoreArrayRead(vec,&x);
 37:   }
 38: #else
 40: #endif
 41:   return(0);
 42: }

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

 47:    Logically Collective on Vec

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

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

 55:    Level: advanced

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

 61: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs()
 62: @*/
 63: PetscErrorCode  VecMaxPointwiseDivide(Vec x,Vec y,PetscReal *max)
 64: {

 74:   VecCheckSameSize(x,1,y,2);
 75:   (*x->ops->maxpointwisedivide)(x,y,max);
 76:   return(0);
 77: }

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

 82:    Collective on Vec

 84:    Input Parameters:
 85: .  x, y - the vectors

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

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

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

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

106:    Level: intermediate


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

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

124:   PetscLogEventBegin(VEC_Dot,x,y,0,0);
125:   (*x->ops->dot)(x,y,val);
126:   PetscLogEventEnd(VEC_Dot,x,y,0,0);
127:   return(0);
128: }

130: /*@
131:    VecDotRealPart - Computes the real part of the vector dot product.

133:    Collective on Vec

135:    Input Parameters:
136: .  x, y - the vectors

138:    Output Parameter:
139: .  val - the real part of the dot product;

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

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

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

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

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

156:    Level: intermediate


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

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

172: /*@
173:    VecNorm  - Computes the vector norm.

175:    Collective on Vec

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

183:    Output Parameter:
184: .  val - the norm

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

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

196:    Level: intermediate

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


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

207: @*/

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


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

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

238:    Not Collective

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

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

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

255:    Level: intermediate

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

262:    Compile Option:
263:    PETSC_HAVE_SLOW_BLAS_NORM2 will cause a C (loop unrolled) version of the norm to be used, rather
264:  than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
265:  (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


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

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

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

328:    Collective on Vec

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

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

337:    Notes:
338:    Returns the value PETSC_MIN_REAL and p = -1 if the vector is of length 0.

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


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

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

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

363:    Collective on Vec

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

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

372:    Level: intermediate

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

377:    This returns the smallest index with the minumum value


380: .seealso: VecMax()
381: @*/
382: PetscErrorCode  VecMin(Vec x,PetscInt *p,PetscReal *val)
383: {

390:   PetscLogEventBegin(VEC_Min,x,0,0,0);
391:   (*x->ops->min)(x,p,val);
392:   PetscLogEventEnd(VEC_Min,x,0,0,0);
393:   return(0);
394: }

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

400:    Collective on Vec

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

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

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

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

417:    Level: intermediate

419: .seealso: VecDot(), VecMTDot()
420: @*/
421: PetscErrorCode  VecTDot(Vec x,Vec y,PetscScalar *val)
422: {

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

434:   PetscLogEventBegin(VEC_TDot,x,y,0,0);
435:   (*x->ops->tdot)(x,y,val);
436:   PetscLogEventEnd(VEC_TDot,x,y,0,0);
437:   return(0);
438: }

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

443:    Not collective on Vec

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

449:    Output Parameter:
450: .  x - the scaled vector

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

456:    Level: intermediate


459: @*/
460: PetscErrorCode  VecScale(Vec x, PetscScalar alpha)
461: {
462:   PetscReal      norms[4] = {0.0,0.0,0.0, 0.0};
463:   PetscBool      flgs[4];
465:   PetscInt       i;

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

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

494:    Logically Collective on Vec

496:    Input Parameters:
497: +  x  - the vector
498: -  alpha - the scalar

500:    Output Parameter:
501: .  x  - the vector

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

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

513:    Level: beginner

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

517: @*/
518: PetscErrorCode  VecSet(Vec x,PetscScalar alpha)
519: {
520:   PetscReal      val;

526:   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()");
528:   VecSetErrorIfLocked(x,1);

530:   PetscLogEventBegin(VEC_Set,x,0,0,0);
531:   (*x->ops->set)(x,alpha);
532:   PetscLogEventEnd(VEC_Set,x,0,0,0);
533:   PetscObjectStateIncrease((PetscObject)x);

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


555: /*@
556:    VecAXPY - Computes y = alpha x + y.

558:    Logically Collective on Vec

560:    Input Parameters:
561: +  alpha - the scalar
562: -  x, y  - the vectors

564:    Output Parameter:
565: .  y - output vector

567:    Level: intermediate

569:    Notes:
570:     x and y MUST be different vectors
571:     This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine

573: $    VecAXPY(y,alpha,x)                   y = alpha x           +      y
574: $    VecAYPX(y,beta,x)                    y =       x           + beta y
575: $    VecAXPBY(y,alpha,beta,x)             y = alpha x           + beta y
576: $    VecWAXPY(w,alpha,x,y)                w = alpha x           +      y
577: $    VecAXPBYPCZ(w,alpha,beta,gamma,x,y)  z = alpha x           + beta y + gamma z
578: $    VecMAXPY(y,nv,alpha[],x[])           y = sum alpha[i] x[i] +      y


581: .seealso:  VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPBYPCZ(), VecAXPBY()
582: @*/
583: PetscErrorCode  VecAXPY(Vec y,PetscScalar alpha,Vec x)
584: {

593:   VecCheckSameSize(x,1,y,3);
594:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
596:   if (alpha == (PetscScalar)0.0) return(0);
597:   VecSetErrorIfLocked(y,1);

599:   VecLockReadPush(x);
600:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
601:   (*y->ops->axpy)(y,alpha,x);
602:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
603:   VecLockReadPop(x);
604:   PetscObjectStateIncrease((PetscObject)y);
605:   return(0);
606: }

608: /*@
609:    VecAXPBY - Computes y = alpha x + beta y.

611:    Logically Collective on Vec

613:    Input Parameters:
614: +  alpha,beta - the scalars
615: -  x, y  - the vectors

617:    Output Parameter:
618: .  y - output vector

620:    Level: intermediate

622:    Notes:
623:     x and y MUST be different vectors
624:     The implementation is optimized for alpha and/or beta values of 0.0 and 1.0


627: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ()
628: @*/
629: PetscErrorCode  VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
630: {

639:   VecCheckSameSize(y,1,x,4);
640:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
643:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) return(0);
644:   VecSetErrorIfLocked(y,1);
645:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
646:   (*y->ops->axpby)(y,alpha,beta,x);
647:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
648:   PetscObjectStateIncrease((PetscObject)y);
649:   return(0);
650: }

652: /*@
653:    VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z

655:    Logically Collective on Vec

657:    Input Parameters:
658: +  alpha,beta, gamma - the scalars
659: -  x, y, z  - the vectors

661:    Output Parameter:
662: .  z - output vector

664:    Level: intermediate

666:    Notes:
667:     x, y and z must be different vectors
668:     The implementation is optimized for alpha of 1.0 and gamma of 1.0 or 0.0


671: .seealso:  VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBY()
672: @*/
673: PetscErrorCode  VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
674: {

686:   VecCheckSameSize(x,1,y,5);
687:   VecCheckSameSize(x,1,z,6);
688:   if (x == y || x == z) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
689:   if (y == z) SETERRQ(PetscObjectComm((PetscObject)y),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
693:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) return(0);
694:   VecSetErrorIfLocked(z,1);

696:   PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
697:   (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
698:   PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
699:   PetscObjectStateIncrease((PetscObject)z);
700:   return(0);
701: }

703: /*@
704:    VecAYPX - Computes y = x + beta y.

706:    Logically Collective on Vec

708:    Input Parameters:
709: +  beta - the scalar
710: -  x, y  - the vectors

712:    Output Parameter:
713: .  y - output vector

715:    Level: intermediate

717:    Notes:
718:     x and y MUST be different vectors
719:     The implementation is optimized for beta of -1.0, 0.0, and 1.0


722: .seealso:  VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
723: @*/
724: PetscErrorCode  VecAYPX(Vec y,PetscScalar beta,Vec x)
725: {

734:   VecCheckSameSize(x,1,y,3);
735:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
737:   VecSetErrorIfLocked(y,1);

739:   PetscLogEventBegin(VEC_AYPX,x,y,0,0);
740:    (*y->ops->aypx)(y,beta,x);
741:   PetscLogEventEnd(VEC_AYPX,x,y,0,0);
742:   PetscObjectStateIncrease((PetscObject)y);
743:   return(0);
744: }


747: /*@
748:    VecWAXPY - Computes w = alpha x + y.

750:    Logically Collective on Vec

752:    Input Parameters:
753: +  alpha - the scalar
754: -  x, y  - the vectors

756:    Output Parameter:
757: .  w - the result

759:    Level: intermediate

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


766: .seealso: VecAXPY(), VecAYPX(), VecAXPBY(), VecMAXPY(), VecAXPBYPCZ()
767: @*/
768: PetscErrorCode  VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
769: {

781:   VecCheckSameSize(x,3,y,4);
782:   VecCheckSameSize(x,3,w,1);
783:   if (w == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector y, suggest VecAXPY()");
784:   if (w == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector x, suggest VecAYPX()");
786:   VecSetErrorIfLocked(w,1);

788:   PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
789:    (*w->ops->waxpy)(w,alpha,x,y);
790:   PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
791:   PetscObjectStateIncrease((PetscObject)w);
792:   return(0);
793: }


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

799:    Not Collective

801:    Input Parameters:
802: +  x - vector to insert in
803: .  ni - number of elements to add
804: .  ix - indices where to add
805: .  y - array of values
806: -  iora - either INSERT_VALUES or ADD_VALUES, where
807:    ADD_VALUES adds values to any existing entries, and
808:    INSERT_VALUES replaces existing entries with new values

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

813:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
814:    options cannot be mixed without intervening calls to the assembly
815:    routines.

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

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

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

828:    Level: beginner

830: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
831:            VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
832: @*/
833: PetscErrorCode  VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
834: {

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

844:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
845:   (*x->ops->setvalues)(x,ni,ix,y,iora);
846:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
847:   PetscObjectStateIncrease((PetscObject)x);
848:   return(0);
849: }

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

855:     Not Collective

857:    Input Parameters:
858: +  x - vector to get values from
859: .  ni - number of elements to get
860: -  ix - indices where to get them from (in global 1d numbering)

862:    Output Parameter:
863: .   y - array of values

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

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

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

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

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

878:    Level: beginner

880: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues()
881: @*/
882: PetscErrorCode  VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
883: {

888:   if (!ni) return(0);
892:   (*x->ops->getvalues)(x,ni,ix,y);
893:   return(0);
894: }

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

899:    Not Collective

901:    Input Parameters:
902: +  x - vector to insert in
903: .  ni - number of blocks to add
904: .  ix - indices where to add in block count, rather than element count
905: .  y - array of values
906: -  iora - either INSERT_VALUES or ADD_VALUES, where
907:    ADD_VALUES adds values to any existing entries, and
908:    INSERT_VALUES replaces existing entries with new values

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

914:    Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
915:    options cannot be mixed without intervening calls to the assembly
916:    routines.

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

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

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

928:    Level: intermediate

930: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
931:            VecSetValues()
932: @*/
933: PetscErrorCode  VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
934: {

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

944:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
945:   (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
946:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
947:   PetscObjectStateIncrease((PetscObject)x);
948:   return(0);
949: }


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

956:    Not Collective

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

967:    Level: intermediate

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

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

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

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

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

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

996:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
997:   if (!x->ops->setvalueslocal) {
998:     if (!x->map->mapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMapping()");
999:     if (ni > 128) {
1000:       PetscMalloc1(ni,&lix);
1001:     }
1002:     ISLocalToGlobalMappingApply(x->map->mapping,ni,(PetscInt*)ix,lix);
1003:     (*x->ops->setvalues)(x,ni,lix,y,iora);
1004:     if (ni > 128) {
1005:       PetscFree(lix);
1006:     }
1007:   } else {
1008:     (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
1009:   }
1010:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1011:   PetscObjectStateIncrease((PetscObject)x);
1012:   return(0);
1013: }

1015: /*@
1016:    VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1017:    using a local ordering of the nodes.

1019:    Not Collective

1021:    Input Parameters:
1022: +  x - vector to insert in
1023: .  ni - number of blocks to add
1024: .  ix - indices where to add in block count, not element count
1025: .  y - array of values
1026: -  iora - either INSERT_VALUES or ADD_VALUES, where
1027:    ADD_VALUES adds values to any existing entries, and
1028:    INSERT_VALUES replaces existing entries with new values

1030:    Level: intermediate

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

1036:    Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1037:    options cannot be mixed without intervening calls to the assembly
1038:    routines.

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

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


1046: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
1047:            VecSetLocalToGlobalMapping()
1048: @*/
1049: PetscErrorCode  VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1050: {
1052:   PetscInt       lixp[128],*lix = lixp;

1056:   if (!ni) return(0);
1060:   if (ni > 128) {
1061:     PetscMalloc1(ni,&lix);
1062:   }

1064:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1065:   ISLocalToGlobalMappingApplyBlock(x->map->mapping,ni,(PetscInt*)ix,lix);
1066:   (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1067:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1068:   if (ni > 128) {
1069:     PetscFree(lix);
1070:   }
1071:   PetscObjectStateIncrease((PetscObject)x);
1072:   return(0);
1073: }

1075: /*@
1076:    VecMTDot - Computes indefinite vector multiple dot products.
1077:    That is, it does NOT use the complex conjugate.

1079:    Collective on Vec

1081:    Input Parameters:
1082: +  x - one vector
1083: .  nv - number of vectors
1084: -  y - array of vectors.  Note that vectors are pointers

1086:    Output Parameter:
1087: .  val - array of the dot products

1089:    Notes for Users of Complex Numbers:
1090:    For complex vectors, VecMTDot() computes the indefinite form
1091: $      val = (x,y) = y^T x,
1092:    where y^T denotes the transpose of y.

1094:    Use VecMDot() for the inner product
1095: $      val = (x,y) = y^H x,
1096:    where y^H denotes the conjugate transpose of y.

1098:    Level: intermediate


1101: .seealso: VecMDot(), VecTDot()
1102: @*/
1103: PetscErrorCode  VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1104: {

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

1119:   PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1120:   (*x->ops->mtdot)(x,nv,y,val);
1121:   PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1122:   return(0);
1123: }

1125: /*@
1126:    VecMDot - Computes vector multiple dot products.

1128:    Collective on Vec

1130:    Input Parameters:
1131: +  x - one vector
1132: .  nv - number of vectors
1133: -  y - array of vectors.

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

1138:    Notes for Users of Complex Numbers:
1139:    For complex vectors, VecMDot() computes
1140: $     val = (x,y) = y^H x,
1141:    where y^H denotes the conjugate transpose of y.

1143:    Use VecMTDot() for the indefinite form
1144: $     val = (x,y) = y^T x,
1145:    where y^T denotes the transpose of y.

1147:    Level: intermediate


1150: .seealso: VecMTDot(), VecDot()
1151: @*/
1152: PetscErrorCode  VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1153: {

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

1169:   PetscLogEventBegin(VEC_MDot,x,*y,0,0);
1170:   (*x->ops->mdot)(x,nv,y,val);
1171:   PetscLogEventEnd(VEC_MDot,x,*y,0,0);
1172:   return(0);
1173: }

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

1178:    Logically Collective on Vec

1180:    Input Parameters:
1181: +  nv - number of scalars and x-vectors
1182: .  alpha - array of scalars
1183: .  y - one vector
1184: -  x - array of vectors

1186:    Level: intermediate

1188:    Notes:
1189:     y cannot be any of the x vectors

1191: .seealso:  VecAYPX(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
1192: @*/
1193: PetscErrorCode  VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1194: {
1196:   PetscInt       i;
1197:   PetscBool      nonzero;

1202:   if (!nv) return(0);
1203:   if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1210:   VecCheckSameSize(y,1,*x,4);
1212:   for (i=0, nonzero = PETSC_FALSE; i<nv && !nonzero; i++) nonzero = (PetscBool)(nonzero || alpha[i] != (PetscScalar)0.0);
1213:   if (!nonzero) return(0);
1214:   VecSetErrorIfLocked(y,1);
1215:   PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1216:   (*y->ops->maxpy)(y,nv,alpha,x);
1217:   PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1218:   PetscObjectStateIncrease((PetscObject)y);
1219:   return(0);
1220: }

1222: /*@
1223:    VecGetSubVector - Gets a vector representing part of another vector

1225:    Collective on X and IS

1227:    Input Arguments:
1228: + X - vector from which to extract a subvector
1229: - is - index set representing portion of X to extract

1231:    Output Arguments:
1232: . Y - subvector corresponding to is

1234:    Level: advanced

1236:    Notes:
1237:    The subvector Y should be returned with VecRestoreSubVector().
1238:    X and is must be defined on the same communicator

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

1244: .seealso: MatCreateSubMatrix()
1245: @*/
1246: PetscErrorCode  VecGetSubVector(Vec X,IS is,Vec *Y)
1247: {
1248:   PetscErrorCode   ierr;
1249:   Vec              Z;

1256:   if (X->ops->getsubvector) {
1257:     (*X->ops->getsubvector)(X,is,&Z);
1258:   } else { /* Default implementation currently does no caching */
1259:     PetscInt  gstart,gend,start;
1260:     PetscBool red[2] = { PETSC_TRUE, PETSC_TRUE };
1261:     PetscInt  n,N,ibs,vbs,bs = -1;

1263:     ISGetLocalSize(is,&n);
1264:     ISGetSize(is,&N);
1265:     ISGetBlockSize(is,&ibs);
1266:     VecGetBlockSize(X,&vbs);
1267:     VecGetOwnershipRange(X,&gstart,&gend);
1268:     ISContiguousLocal(is,gstart,gend,&start,&red[0]);
1269:     /* block size is given by IS if ibs > 1; otherwise, check the vector */
1270:     if (ibs > 1) {
1271:       MPIU_Allreduce(MPI_IN_PLACE,red,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1272:       bs   = ibs;
1273:     } else {
1274:       if (n%vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1275:       MPIU_Allreduce(MPI_IN_PLACE,red,2,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1276:       if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1277:     }
1278:     if (red[0]) { /* We can do a no-copy implementation */
1279:       const PetscScalar *x;
1280:       PetscInt          state = 0;
1281:       PetscBool         isstd;
1282: #if defined(PETSC_HAVE_CUDA)
1283:       PetscBool         iscuda;
1284: #endif

1286:       PetscObjectTypeCompareAny((PetscObject)X,&isstd,VECSEQ,VECMPI,VECSTANDARD,"");
1287: #if defined(PETSC_HAVE_CUDA)
1288:       PetscObjectTypeCompareAny((PetscObject)X,&iscuda,VECSEQCUDA,VECMPICUDA,"");
1289:       if (iscuda) {
1290:         const PetscScalar *x_d;
1291:         PetscMPIInt       size;
1292:         PetscOffloadMask  flg;

1294:         VecCUDAGetArrays_Private(X,&x,&x_d,&flg);
1295:         if (flg == PETSC_OFFLOAD_UNALLOCATED) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for PETSC_OFFLOAD_UNALLOCATED");
1296:         if (n && !x && !x_d) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Missing vector data");
1297:         if (x) x += start;
1298:         if (x_d) x_d += start;
1299:         MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1300:         if (size == 1) {
1301:           VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X),bs,n,x,x_d,&Z);
1302:         } else {
1303:           VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X),bs,n,N,x,x_d,&Z);
1304:         }
1305:         Z->offloadmask = flg;
1306:       } else if (isstd) {
1307: #else
1308:       if (isstd) { /* standard CPU: use CreateWithArray pattern */
1309: #endif
1310:         PetscMPIInt size;

1312:         MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1313:         VecGetArrayRead(X,&x);
1314:         if (x) x += start;
1315:         if (size == 1) {
1316:           VecCreateSeqWithArray(PetscObjectComm((PetscObject)X),bs,n,x,&Z);
1317:         } else {
1318:           VecCreateMPIWithArray(PetscObjectComm((PetscObject)X),bs,n,N,x,&Z);
1319:         }
1320:         VecRestoreArrayRead(X,&x);
1321:       } else { /* default implementation: use place array */
1322:         VecGetArrayRead(X,&x);
1323:         VecCreate(PetscObjectComm((PetscObject)X),&Z);
1324:         VecSetType(Z,((PetscObject)X)->type_name);
1325:         VecSetSizes(Z,n,N);
1326:         VecSetBlockSize(Z,bs);
1327:         VecPlaceArray(Z,x ? x+start : NULL);
1328:         VecRestoreArrayRead(X,&x);
1329:       }

1331:       /* this is relevant only in debug mode */
1332:       VecLockGet(X,&state);
1333:       if (state) {
1334:         VecLockReadPush(Z);
1335:       }
1336:       Z->ops->placearray = NULL;
1337:       Z->ops->replacearray = NULL;
1338:     } else { /* Have to create a scatter and do a copy */
1339:       VecScatter scatter;

1341:       VecCreate(PetscObjectComm((PetscObject)is),&Z);
1342:       VecSetSizes(Z,n,N);
1343:       VecSetBlockSize(Z,bs);
1344:       VecSetType(Z,((PetscObject)X)->type_name);
1345:       VecScatterCreate(X,is,Z,NULL,&scatter);
1346:       VecScatterBegin(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1347:       VecScatterEnd(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1348:       PetscObjectCompose((PetscObject)Z,"VecGetSubVector_Scatter",(PetscObject)scatter);
1349:       VecScatterDestroy(&scatter);
1350:     }
1351:   }
1352:   /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1353:   if (VecGetSubVectorSavedStateId < 0) {PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);}
1354:   PetscObjectComposedDataSetInt((PetscObject)Z,VecGetSubVectorSavedStateId,1);
1355:   *Y   = Z;
1356:   return(0);
1357: }

1359: /*@
1360:    VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()

1362:    Collective on IS

1364:    Input Arguments:
1365: + X - vector from which subvector was obtained
1366: . is - index set representing the subset of X
1367: - Y - subvector being restored

1369:    Level: advanced

1371: .seealso: VecGetSubVector()
1372: @*/
1373: PetscErrorCode  VecRestoreSubVector(Vec X,IS is,Vec *Y)
1374: {

1383:   if (X->ops->restoresubvector) {
1384:     (*X->ops->restoresubvector)(X,is,Y);
1385:   } else {
1386:     PETSC_UNUSED PetscObjectState dummystate = 0;
1387:     PetscBool valid;

1389:     PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,dummystate,valid);
1390:     if (!valid) {
1391:       VecScatter scatter;
1392:       PetscInt   state;

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

1397:       PetscObjectQuery((PetscObject)*Y,"VecGetSubVector_Scatter",(PetscObject*)&scatter);
1398:       if (scatter) {
1399:         VecScatterBegin(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1400:         VecScatterEnd(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1401:       } else {
1402: #if defined(PETSC_HAVE_CUDA)
1403:         PetscBool iscuda;

1405:         PetscObjectTypeCompareAny((PetscObject)*Y,&iscuda,VECSEQCUDA,VECMPICUDA,"");
1406:         if (iscuda) {
1407:           PetscOffloadMask ymask = (*Y)->offloadmask;

1409:           /* The offloadmask of X dictates where to move memory
1410:              If X GPU data is valid, then move Y data on GPU if needed
1411:              Otherwise, move back to the CPU */
1412:           switch (X->offloadmask) {
1413:           case PETSC_OFFLOAD_BOTH:
1414:             if (ymask == PETSC_OFFLOAD_CPU) {
1415:               VecCUDAResetArray(*Y);
1416:             } else if (ymask == PETSC_OFFLOAD_GPU) {
1417:               X->offloadmask = PETSC_OFFLOAD_GPU;
1418:             }
1419:             break;
1420:           case PETSC_OFFLOAD_GPU:
1421:             if (ymask == PETSC_OFFLOAD_CPU) {
1422:               VecCUDAResetArray(*Y);
1423:             }
1424:             break;
1425:           case PETSC_OFFLOAD_CPU:
1426:             if (ymask == PETSC_OFFLOAD_GPU) {
1427:               VecResetArray(*Y);
1428:             }
1429:             break;
1430:           case PETSC_OFFLOAD_UNALLOCATED:
1431:           case PETSC_OFFLOAD_VECKOKKOS:
1432:             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"This should not happen");
1433:             break;
1434:           }
1435:         } else {
1436: #else
1437:         {
1438: #endif
1439:           /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1440:           VecResetArray(*Y);
1441:         }
1442:         PetscObjectStateIncrease((PetscObject)X);
1443:       }
1444:     }
1445:     VecDestroy(Y);
1446:   }
1447:   return(0);
1448: }

1450: /*@
1451:    VecGetLocalVectorRead - Maps the local portion of a vector into a
1452:    vector.  You must call VecRestoreLocalVectorRead() when the local
1453:    vector is no longer needed.

1455:    Not collective.

1457:    Input parameter:
1458: .  v - The vector for which the local vector is desired.

1460:    Output parameter:
1461: .  w - Upon exit this contains the local vector.

1463:    Level: beginner

1465:    Notes:
1466:    This function is similar to VecGetArrayRead() which maps the local
1467:    portion into a raw pointer.  VecGetLocalVectorRead() is usually
1468:    almost as efficient as VecGetArrayRead() but in certain circumstances
1469:    VecGetLocalVectorRead() can be much more efficient than
1470:    VecGetArrayRead().  This is because the construction of a contiguous
1471:    array representing the vector data required by VecGetArrayRead() can
1472:    be an expensive operation for certain vector types.  For example, for
1473:    GPU vectors VecGetArrayRead() requires that the data between device
1474:    and host is synchronized.

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

1479: .seealso: VecRestoreLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1480: @*/
1481: PetscErrorCode VecGetLocalVectorRead(Vec v,Vec w)
1482: {
1484:   PetscScalar    *a;

1489:   VecCheckSameLocalSize(v,1,w,2);
1490:   if (v->ops->getlocalvectorread) {
1491:     (*v->ops->getlocalvectorread)(v,w);
1492:   } else {
1493:     VecGetArrayRead(v,(const PetscScalar**)&a);
1494:     VecPlaceArray(w,a);
1495:   }
1496:   PetscObjectStateIncrease((PetscObject)w);
1497:   VecLockReadPush(v);
1498:   VecLockReadPush(w);
1499:   return(0);
1500: }

1502: /*@
1503:    VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1504:    previously mapped into a vector using VecGetLocalVectorRead().

1506:    Not collective.

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

1512:    Level: beginner

1514: .seealso: VecGetLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1515: @*/
1516: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec w)
1517: {
1519:   PetscScalar    *a;

1524:   if (v->ops->restorelocalvectorread) {
1525:     (*v->ops->restorelocalvectorread)(v,w);
1526:   } else {
1527:     VecGetArrayRead(w,(const PetscScalar**)&a);
1528:     VecRestoreArrayRead(v,(const PetscScalar**)&a);
1529:     VecResetArray(w);
1530:   }
1531:   VecLockReadPop(v);
1532:   VecLockReadPop(w);
1533:   PetscObjectStateIncrease((PetscObject)w);
1534:   return(0);
1535: }

1537: /*@
1538:    VecGetLocalVector - Maps the local portion of a vector into a
1539:    vector.

1541:    Collective on v, not collective on w.

1543:    Input parameter:
1544: .  v - The vector for which the local vector is desired.

1546:    Output parameter:
1547: .  w - Upon exit this contains the local vector.

1549:    Level: beginner

1551:    Notes:
1552:    This function is similar to VecGetArray() which maps the local
1553:    portion into a raw pointer.  VecGetLocalVector() is usually about as
1554:    efficient as VecGetArray() but in certain circumstances
1555:    VecGetLocalVector() can be much more efficient than VecGetArray().
1556:    This is because the construction of a contiguous array representing
1557:    the vector data required by VecGetArray() can be an expensive
1558:    operation for certain vector types.  For example, for GPU vectors
1559:    VecGetArray() requires that the data between device and host is
1560:    synchronized.

1562: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1563: @*/
1564: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1565: {
1567:   PetscScalar    *a;

1572:   VecCheckSameLocalSize(v,1,w,2);
1573:   if (v->ops->getlocalvector) {
1574:     (*v->ops->getlocalvector)(v,w);
1575:   } else {
1576:     VecGetArray(v,&a);
1577:     VecPlaceArray(w,a);
1578:   }
1579:   PetscObjectStateIncrease((PetscObject)w);
1580:   return(0);
1581: }

1583: /*@
1584:    VecRestoreLocalVector - Unmaps the local portion of a vector
1585:    previously mapped into a vector using VecGetLocalVector().

1587:    Logically collective.

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

1593:    Level: beginner

1595: .seealso: VecGetLocalVector(), VecGetLocalVectorRead(), VecRestoreLocalVectorRead(), LocalVectorRead(), VecGetArrayRead(), VecGetArray()
1596: @*/
1597: PetscErrorCode VecRestoreLocalVector(Vec v,Vec w)
1598: {
1600:   PetscScalar    *a;

1605:   if (v->ops->restorelocalvector) {
1606:     (*v->ops->restorelocalvector)(v,w);
1607:   } else {
1608:     VecGetArray(w,&a);
1609:     VecRestoreArray(v,&a);
1610:     VecResetArray(w);
1611:   }
1612:   PetscObjectStateIncrease((PetscObject)w);
1613:   PetscObjectStateIncrease((PetscObject)v);
1614:   return(0);
1615: }

1617: /*@C
1618:    VecGetArray - Returns a pointer to a contiguous array that contains this
1619:    processor's portion of the vector data. For the standard PETSc
1620:    vectors, VecGetArray() returns a pointer to the local data array and
1621:    does not use any copies. If the underlying vector data is not stored
1622:    in a contiguous array this routine will copy the data to a contiguous
1623:    array and return a pointer to that. You MUST call VecRestoreArray()
1624:    when you no longer need access to the array.

1626:    Logically Collective on Vec

1628:    Input Parameter:
1629: .  x - the vector

1631:    Output Parameter:
1632: .  a - location to put pointer to the array

1634:    Fortran Note:
1635:    This routine is used differently from Fortran 77
1636: $    Vec         x
1637: $    PetscScalar x_array(1)
1638: $    PetscOffset i_x
1639: $    PetscErrorCode ierr
1640: $       call VecGetArray(x,x_array,i_x,ierr)
1641: $
1642: $   Access first local entry in vector with
1643: $      value = x_array(i_x + 1)
1644: $
1645: $      ...... other code
1646: $       call VecRestoreArray(x,x_array,i_x,ierr)
1647:    For Fortran 90 see VecGetArrayF90()

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

1652:    Level: beginner

1654: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1655:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1656: @*/
1657: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1658: {

1663:   VecSetErrorIfLocked(x,1);
1664:   if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
1665:     (*x->ops->getarray)(x,a);
1666:   } else if (x->petscnative) { /* VECSTANDARD */
1667:     *a = *((PetscScalar**)x->data);
1668:   } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array for vector type \"%s\"",((PetscObject)x)->type_name);
1669:   return(0);
1670: }

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

1675:    Logically Collective on Vec

1677:    Input Parameters:
1678: +  x - the vector
1679: -  a - location of pointer to array obtained from VecGetArray()

1681:    Level: beginner

1683: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1684:           VecGetArrayPair(), VecRestoreArrayPair()
1685: @*/
1686: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1687: {

1692:   if (x->ops->restorearray) { /* VECNEST, VECCUDA etc */
1693:     (*x->ops->restorearray)(x,a);
1694:   } else if (x->petscnative) { /* VECSTANDARD */
1695:     /* nothing */
1696:   } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot restore array for vector type \"%s\"",((PetscObject)x)->type_name);
1697:   if (a) *a = NULL;
1698:   PetscObjectStateIncrease((PetscObject)x);
1699:   return(0);
1700: }
1701: /*@C
1702:    VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.

1704:    Not Collective

1706:    Input Parameter:
1707: .  x - the vector

1709:    Output Parameter:
1710: .  a - the array

1712:    Level: beginner

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

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

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

1723: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1724: @*/
1725: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1726: {

1731:   if (x->ops->getarray) { /* VECNEST, VECCUDA, VECKOKKOS etc */
1732:     (*x->ops->getarray)(x,(PetscScalar**)a);
1733:   } else if (x->petscnative) { /* VECSTANDARD */
1734:     *a = *((PetscScalar**)x->data);
1735:   } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array read for vector type \"%s\"",((PetscObject)x)->type_name);
1736:   return(0);
1737: }

1739: /*@C
1740:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

1742:    Not Collective

1744:    Input Parameters:
1745: +  vec - the vector
1746: -  array - the array

1748:    Level: beginner

1750: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1751: @*/
1752: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1753: {

1758:   if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
1759:     /* nothing */
1760:   } else if (x->ops->restorearrayread) { /* VECNEST */
1761:     (*x->ops->restorearrayread)(x,a);
1762:   } else { /* No one? */
1763:     (*x->ops->restorearray)(x,(PetscScalar**)a);
1764:   }
1765:   if (a) *a = NULL;
1766:   return(0);
1767: }

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

1774:    Logically Collective on Vec

1776:    Input Parameter:
1777: .  x - the vector

1779:    Output Parameter:
1780: .  a - location to put pointer to the array

1782:    Level: intermediate

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

1787:    Concepts: vector^accessing local values

1789: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1790:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArray(), VecRestoreArrayWrite()
1791: @*/
1792: PetscErrorCode VecGetArrayWrite(Vec x,PetscScalar **a)
1793: {

1798:   VecSetErrorIfLocked(x,1);
1799:   if (x->ops->getarraywrite) {
1800:     (*x->ops->getarraywrite)(x,a);
1801:   } else {
1802:     VecGetArray(x,a);
1803:   }
1804:   return(0);
1805: }

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

1810:    Logically Collective on Vec

1812:    Input Parameters:
1813: +  x - the vector
1814: -  a - location of pointer to array obtained from VecGetArray()

1816:    Level: beginner

1818: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1819:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite()
1820: @*/
1821: PetscErrorCode VecRestoreArrayWrite(Vec x,PetscScalar **a)
1822: {

1827:   if (x->ops->restorearraywrite) {
1828:     (*x->ops->restorearraywrite)(x,a);
1829:   } else if (x->ops->restorearray) {
1830:     (*x->ops->restorearray)(x,a);
1831:   }
1832:   if (a) *a = NULL;
1833:   PetscObjectStateIncrease((PetscObject)x);
1834:   return(0);
1835: }

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

1842:    Logically Collective on Vec

1844:    Input Parameters:
1845: +  x - the vectors
1846: -  n - the number of vectors

1848:    Output Parameter:
1849: .  a - location to put pointer to the array

1851:    Fortran Note:
1852:    This routine is not supported in Fortran.

1854:    Level: intermediate

1856: .seealso: VecGetArray(), VecRestoreArrays()
1857: @*/
1858: PetscErrorCode  VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1859: {
1861:   PetscInt       i;
1862:   PetscScalar    **q;

1868:   if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1869:   PetscMalloc1(n,&q);
1870:   for (i=0; i<n; ++i) {
1871:     VecGetArray(x[i],&q[i]);
1872:   }
1873:   *a = q;
1874:   return(0);
1875: }

1877: /*@C
1878:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1879:    has been called.

1881:    Logically Collective on Vec

1883:    Input Parameters:
1884: +  x - the vector
1885: .  n - the number of vectors
1886: -  a - location of pointer to arrays obtained from VecGetArrays()

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

1894:    Fortran Note:
1895:    This routine is not supported in Fortran.

1897:    Level: intermediate

1899: .seealso: VecGetArrays(), VecRestoreArray()
1900: @*/
1901: PetscErrorCode  VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1902: {
1904:   PetscInt       i;
1905:   PetscScalar    **q = *a;


1912:   for (i=0; i<n; ++i) {
1913:     VecRestoreArray(x[i],&q[i]);
1914:   }
1915:   PetscFree(q);
1916:   return(0);
1917: }

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

1924:    Logically Collective on Vec

1926:    Input Parameter:
1927: .  x - the vector

1929:    Output Parameters:
1930: +  a - location to put pointer to the array
1931: -  mtype - memory type of the array

1933:    Level: beginner

1935: .seealso: VecRestoreArrayAndMemType(), VecRestoreArrayAndMemType(), VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(),
1936:           VecPlaceArray(), VecGetArray2d(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1937: @*/
1938: PetscErrorCode VecGetArrayAndMemType(Vec x,PetscScalar **a,PetscMemType *mtype)
1939: {

1945:   VecSetErrorIfLocked(x,1);
1946:   if (x->ops->getarrayandmemtype) { /* VECCUDA, VECKOKKOS etc */
1947:     (*x->ops->getarrayandmemtype)(x,a,mtype);
1948:   } else { /* VECSTANDARD, VECNEST, VECVIENNACL */
1949:     VecGetArray(x,a);
1950:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
1951:   }
1952:   return(0);
1953: }

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

1958:    Logically Collective on Vec

1960:    Input Parameters:
1961: +  x - the vector
1962: -  a - location of pointer to array obtained from VecGetArrayAndMemType()

1964:    Level: beginner

1966: .seealso: VecGetArrayAndMemType(), VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(),
1967:           VecPlaceArray(), VecRestoreArray2d(), VecGetArrayPair(), VecRestoreArrayPair()
1968: @*/
1969: PetscErrorCode VecRestoreArrayAndMemType(Vec x,PetscScalar **a)
1970: {

1976:   if (x->ops->restorearrayandmemtype) { /* VECCUDA, VECKOKKOS etc */
1977:     (*x->ops->restorearrayandmemtype)(x,a);
1978:   } else if (x->ops->restorearray) { /* VECNEST, VECVIENNACL */
1979:     (*x->ops->restorearray)(x,a);
1980:   } /* VECSTANDARD does nothing */
1981:   if (a) *a = NULL;
1982:   PetscObjectStateIncrease((PetscObject)x);
1983:   return(0);
1984: }

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

1991:    Not Collective

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

1996:    Output Parameters:
1997: +  a - the array
1998: -  mtype - memory type of the array

2000:    Level: beginner

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


2006: .seealso: VecRestoreArrayReadAndMemType(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayAndMemType()
2007: @*/
2008: PetscErrorCode VecGetArrayReadAndMemType(Vec x,const PetscScalar **a,PetscMemType *mtype)
2009: {

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

2019:   if (x->ops->getarrayandmemtype) { /* VECCUDA, VECKOKKOS etc, though they are also petscnative */
2020:     (*x->ops->getarrayandmemtype)(x,(PetscScalar**)a,mtype);
2021:   } else if (x->ops->getarray) { /* VECNEST, VECVIENNACL */
2022:     (*x->ops->getarray)(x,(PetscScalar**)a);
2023:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2024:   } else if (x->petscnative) { /* VECSTANDARD */
2025:     *a = *((PetscScalar**)x->data);
2026:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2027:   } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array read in place for vector type \"%s\"",((PetscObject)x)->type_name);
2028:   return(0);
2029: }

2031: /*@C
2032:    VecRestoreArrayReadAndMemType - Restore array obtained with VecGetArrayReadAndMemType()

2034:    Not Collective

2036:    Input Parameters:
2037: +  vec - the vector
2038: -  array - the array

2040:    Level: beginner

2042: .seealso: VecGetArrayReadAndMemType(), VecRestoreArrayAndMemType(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
2043: @*/
2044: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x,const PetscScalar **a)
2045: {

2051:   if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS, VECVIENNACL etc */
2052:     /* nothing */
2053:   } else if (x->ops->restorearrayread) { /* VECNEST */
2054:     (*x->ops->restorearrayread)(x,a);
2055:   } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot restore array read in place for vector type \"%s\"",((PetscObject)x)->type_name);
2056:   if (a) *a = NULL;
2057:   return(0);
2058: }

2060: /*@
2061:    VecPlaceArray - Allows one to replace the array in a vector with an
2062:    array provided by the user. This is useful to avoid copying an array
2063:    into a vector.

2065:    Not Collective

2067:    Input Parameters:
2068: +  vec - the vector
2069: -  array - the array

2071:    Notes:
2072:    You can return to the original array with a call to VecResetArray()

2074:    Level: developer

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

2078: @*/
2079: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
2080: {

2087:   if (vec->ops->placearray) {
2088:     (*vec->ops->placearray)(vec,array);
2089:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
2090:   PetscObjectStateIncrease((PetscObject)vec);
2091:   return(0);
2092: }

2094: /*@C
2095:    VecReplaceArray - Allows one to replace the array in a vector with an
2096:    array provided by the user. This is useful to avoid copying an array
2097:    into a vector.

2099:    Not Collective

2101:    Input Parameters:
2102: +  vec - the vector
2103: -  array - the array

2105:    Notes:
2106:    This permanently replaces the array and frees the memory associated
2107:    with the old array.

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

2112:    Not supported from Fortran

2114:    Level: developer

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

2118: @*/
2119: PetscErrorCode  VecReplaceArray(Vec vec,const PetscScalar array[])
2120: {

2126:   if (vec->ops->replacearray) {
2127:     (*vec->ops->replacearray)(vec,array);
2128:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
2129:   PetscObjectStateIncrease((PetscObject)vec);
2130:   return(0);
2131: }


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

2137:    This function has semantics similar to VecGetArray():  the pointer
2138:    returned by this function points to a consistent view of the vector
2139:    data.  This may involve a copy operation of data from the host to the
2140:    device if the data on the device is out of date.  If the device
2141:    memory hasn't been allocated previously it will be allocated as part
2142:    of this function call.  VecCUDAGetArray() assumes that
2143:    the user will modify the vector data.  This is similar to
2144:    intent(inout) in fortran.

2146:    The CUDA device pointer has to be released by calling
2147:    VecCUDARestoreArray().  Upon restoring the vector data
2148:    the data on the host will be marked as out of date.  A subsequent
2149:    access of the host data will thus incur a data transfer from the
2150:    device to the host.


2153:    Input Parameter:
2154: .  v - the vector

2156:    Output Parameter:
2157: .  a - the CUDA device pointer

2159:    Fortran note:
2160:    This function is not currently available from Fortran.

2162:    Level: intermediate

2164: .seealso: VecCUDARestoreArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2165: @*/
2166: PETSC_EXTERN PetscErrorCode VecCUDAGetArray(Vec v, PetscScalar **a)
2167: {
2170:  #if defined(PETSC_HAVE_CUDA)
2171:   {
2173:     VecCUDACopyToGPU(v);
2174:     *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2175:   }
2176:  #endif
2177:   return(0);
2178: }

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

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

2187:    Input Parameter:
2188: +  v - the vector
2189: -  a - the CUDA device pointer.  This pointer is invalid after
2190:        VecCUDARestoreArray() returns.

2192:    Fortran note:
2193:    This function is not currently available from Fortran.

2195:    Level: intermediate

2197: .seealso: VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2198: @*/
2199: PETSC_EXTERN PetscErrorCode VecCUDARestoreArray(Vec v, PetscScalar **a)
2200: {

2205: #if defined(PETSC_HAVE_CUDA)
2206:   v->offloadmask = PETSC_OFFLOAD_GPU;
2207: #endif
2208:   PetscObjectStateIncrease((PetscObject)v);
2209:   return(0);
2210: }

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

2215:    This function is analogous to VecGetArrayRead():  The pointer
2216:    returned by this function points to a consistent view of the vector
2217:    data.  This may involve a copy operation of data from the host to the
2218:    device if the data on the device is out of date.  If the device
2219:    memory hasn't been allocated previously it will be allocated as part
2220:    of this function call.  VecCUDAGetArrayRead() assumes that the
2221:    user will not modify the vector data.  This is analgogous to
2222:    intent(in) in Fortran.

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

2230:    Input Parameter:
2231: .  v - the vector

2233:    Output Parameter:
2234: .  a - the CUDA pointer.

2236:    Fortran note:
2237:    This function is not currently available from Fortran.

2239:    Level: intermediate

2241: .seealso: VecCUDARestoreArrayRead(), VecCUDAGetArray(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2242: @*/
2243: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayRead(Vec v,const PetscScalar** a)
2244: {
2247:    VecCUDAGetArray(v,(PetscScalar**)a);
2248:    return(0);
2249: }

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

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

2259:    Input Parameter:
2260: +  v - the vector
2261: -  a - the CUDA device pointer.  This pointer is invalid after
2262:        VecCUDARestoreArrayRead() returns.

2264:    Fortran note:
2265:    This function is not currently available from Fortran.

2267:    Level: intermediate

2269: .seealso: VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2270: @*/
2271: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayRead(Vec v, const PetscScalar **a)
2272: {
2275:   *a = NULL;
2276:   return(0);
2277: }

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

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

2288:    The device pointer needs to be released with
2289:    VecCUDARestoreArrayWrite().  When the pointer is released the
2290:    host data of the vector is marked as out of data.  Subsequent access
2291:    of the host data with e.g. VecGetArray() incurs a device to host data
2292:    transfer.


2295:    Input Parameter:
2296: .  v - the vector

2298:    Output Parameter:
2299: .  a - the CUDA pointer

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

2304:    Level: advanced

2306: .seealso: VecCUDARestoreArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2307: @*/
2308: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayWrite(Vec v, PetscScalar **a)
2309: {
2312:  #if defined(PETSC_HAVE_CUDA)
2313:   {
2315:     VecCUDAAllocateCheck(v);
2316:     *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2317:   }
2318:  #endif
2319:   return(0);
2320: }

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

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

2329:    Input Parameter:
2330: +  v - the vector
2331: -  a - the CUDA device pointer.  This pointer is invalid after
2332:        VecCUDARestoreArrayWrite() returns.

2334:    Fortran note:
2335:    This function is not currently available from Fortran.

2337:    Level: intermediate

2339: .seealso: VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2340: @*/
2341: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayWrite(Vec v, PetscScalar **a)
2342: {

2347:  #if defined(PETSC_HAVE_CUDA)
2348:   v->offloadmask = PETSC_OFFLOAD_GPU;
2349:   a              = NULL;
2350:  #endif
2351:   PetscObjectStateIncrease((PetscObject)v);
2352:   return(0);
2353: }

2355: /*@C
2356:    VecCUDAPlaceArray - Allows one to replace the GPU array in a vector with a
2357:    GPU array provided by the user. This is useful to avoid copying an
2358:    array into a vector.

2360:    Not Collective

2362:    Input Parameters:
2363: +  vec - the vector
2364: -  array - the GPU array

2366:    Notes:
2367:    You can return to the original GPU array with a call to VecCUDAResetArray()
2368:    It is not possible to use VecCUDAPlaceArray() and VecPlaceArray() at the
2369:    same time on the same vector.

2371:    Level: developer

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

2375: @*/
2376: PetscErrorCode VecCUDAPlaceArray(Vec vin,const PetscScalar a[])
2377: {

2382: #if defined(PETSC_HAVE_CUDA)
2383:   VecCUDACopyToGPU(vin);
2384:   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()");
2385:   ((Vec_Seq*)vin->data)->unplacedarray  = (PetscScalar *) ((Vec_CUDA*)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2386:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2387:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2388: #endif
2389:   PetscObjectStateIncrease((PetscObject)vin);
2390:   return(0);
2391: }

2393: /*@C
2394:    VecCUDAReplaceArray - Allows one to replace the GPU array in a vector
2395:    with a GPU array provided by the user. This is useful to avoid copying
2396:    a GPU array into a vector.

2398:    Not Collective

2400:    Input Parameters:
2401: +  vec - the vector
2402: -  array - the GPU array

2404:    Notes:
2405:    This permanently replaces the GPU array and frees the memory associated
2406:    with the old GPU array.

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

2411:    Not supported from Fortran

2413:    Level: developer

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

2417: @*/
2418: PetscErrorCode VecCUDAReplaceArray(Vec vin,const PetscScalar a[])
2419: {
2420: #if defined(PETSC_HAVE_CUDA)
2421:   cudaError_t err;
2422: #endif

2427: #if defined(PETSC_HAVE_CUDA)
2428:   err = cudaFree(((Vec_CUDA*)vin->spptr)->GPUarray);CHKERRCUDA(err);
2429:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2430:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2431: #endif
2432:   PetscObjectStateIncrease((PetscObject)vin);
2433:   return(0);
2434: }

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

2440:    Not Collective

2442:    Input Parameters:
2443: .  vec - the vector

2445:    Level: developer

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

2449: @*/
2450: PetscErrorCode VecCUDAResetArray(Vec vin)
2451: {

2456: #if defined(PETSC_HAVE_CUDA)
2457:   VecCUDACopyToGPU(vin);
2458:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar *) ((Vec_Seq*)vin->data)->unplacedarray;
2459:   ((Vec_Seq*)vin->data)->unplacedarray = 0;
2460:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2461: #endif
2462:   PetscObjectStateIncrease((PetscObject)vin);
2463:   return(0);
2464: }

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

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

2473:     Collective on Vec

2475:     Input Parameters:
2476: +   x - a vector to mimic
2477: -   n - the number of vectors to obtain

2479:     Output Parameters:
2480: +   y - Fortran90 pointer to the array of vectors
2481: -   ierr - error code

2483:     Example of Usage:
2484: .vb
2485: #include <petsc/finclude/petscvec.h>
2486:     use petscvec

2488:     Vec x
2489:     Vec, pointer :: y(:)
2490:     ....
2491:     call VecDuplicateVecsF90(x,2,y,ierr)
2492:     call VecSet(y(2),alpha,ierr)
2493:     call VecSet(y(2),alpha,ierr)
2494:     ....
2495:     call VecDestroyVecsF90(2,y,ierr)
2496: .ve

2498:     Notes:
2499:     Not yet supported for all F90 compilers

2501:     Use VecDestroyVecsF90() to free the space.

2503:     Level: beginner

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

2507: M*/

2509: /*MC
2510:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2511:     VecGetArrayF90().

2513:     Synopsis:
2514:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2516:     Logically Collective on Vec

2518:     Input Parameters:
2519: +   x - vector
2520: -   xx_v - the Fortran90 pointer to the array

2522:     Output Parameter:
2523: .   ierr - error code

2525:     Example of Usage:
2526: .vb
2527: #include <petsc/finclude/petscvec.h>
2528:     use petscvec

2530:     PetscScalar, pointer :: xx_v(:)
2531:     ....
2532:     call VecGetArrayF90(x,xx_v,ierr)
2533:     xx_v(3) = a
2534:     call VecRestoreArrayF90(x,xx_v,ierr)
2535: .ve

2537:     Level: beginner

2539: .seealso:  VecGetArrayF90(), VecGetArray(), VecRestoreArray(), UsingFortran, VecRestoreArrayReadF90()

2541: M*/

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

2546:     Synopsis:
2547:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

2549:     Collective on Vec

2551:     Input Parameters:
2552: +   n - the number of vectors previously obtained
2553: -   x - pointer to array of vector pointers

2555:     Output Parameter:
2556: .   ierr - error code

2558:     Notes:
2559:     Not yet supported for all F90 compilers

2561:     Level: beginner

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

2565: M*/

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

2573:     Synopsis:
2574:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2576:     Logically Collective on Vec

2578:     Input Parameter:
2579: .   x - vector

2581:     Output Parameters:
2582: +   xx_v - the Fortran90 pointer to the array
2583: -   ierr - error code

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

2590:     PetscScalar, pointer :: xx_v(:)
2591:     ....
2592:     call VecGetArrayF90(x,xx_v,ierr)
2593:     xx_v(3) = a
2594:     call VecRestoreArrayF90(x,xx_v,ierr)
2595: .ve

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

2599:     Level: beginner

2601: .seealso:  VecRestoreArrayF90(), VecGetArray(), VecRestoreArray(), VecGetArrayReadF90(), UsingFortran

2603: M*/

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

2611:     Synopsis:
2612:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2614:     Logically Collective on Vec

2616:     Input Parameter:
2617: .   x - vector

2619:     Output Parameters:
2620: +   xx_v - the Fortran90 pointer to the array
2621: -   ierr - error code

2623:     Example of Usage:
2624: .vb
2625: #include <petsc/finclude/petscvec.h>
2626:     use petscvec

2628:     PetscScalar, pointer :: xx_v(:)
2629:     ....
2630:     call VecGetArrayReadF90(x,xx_v,ierr)
2631:     a = xx_v(3)
2632:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2633: .ve

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

2637:     Level: beginner

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

2641: M*/

2643: /*MC
2644:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2645:     VecGetArrayReadF90().

2647:     Synopsis:
2648:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2650:     Logically Collective on Vec

2652:     Input Parameters:
2653: +   x - vector
2654: -   xx_v - the Fortran90 pointer to the array

2656:     Output Parameter:
2657: .   ierr - error code

2659:     Example of Usage:
2660: .vb
2661: #include <petsc/finclude/petscvec.h>
2662:     use petscvec

2664:     PetscScalar, pointer :: xx_v(:)
2665:     ....
2666:     call VecGetArrayReadF90(x,xx_v,ierr)
2667:     a = xx_v(3)
2668:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2669: .ve

2671:     Level: beginner

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

2675: M*/

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

2682:    Logically Collective

2684:    Input Parameter:
2685: +  x - the vector
2686: .  m - first dimension of two dimensional array
2687: .  n - second dimension of two dimensional array
2688: .  mstart - first index you will use in first coordinate direction (often 0)
2689: -  nstart - first index in the second coordinate direction (often 0)

2691:    Output Parameter:
2692: .  a - location to put pointer to the array

2694:    Level: developer

2696:   Notes:
2697:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2698:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2699:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2700:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

2704: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2705:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2706:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2707: @*/
2708: PetscErrorCode  VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2709: {
2711:   PetscInt       i,N;
2712:   PetscScalar    *aa;

2718:   VecGetLocalSize(x,&N);
2719:   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);
2720:   VecGetArray(x,&aa);

2722:   PetscMalloc1(m,a);
2723:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2724:   *a -= mstart;
2725:   return(0);
2726: }

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

2733:    Logically Collective

2735:    Input Parameter:
2736: +  x - the vector
2737: .  m - first dimension of two dimensional array
2738: .  n - second dimension of two dimensional array
2739: .  mstart - first index you will use in first coordinate direction (often 0)
2740: -  nstart - first index in the second coordinate direction (often 0)

2742:    Output Parameter:
2743: .  a - location to put pointer to the array

2745:    Level: developer

2747:   Notes:
2748:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2749:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2750:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2751:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

2757: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2758:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2759:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2760: @*/
2761: PetscErrorCode  VecGetArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2762: {
2764:   PetscInt       i,N;
2765:   PetscScalar    *aa;

2771:   VecGetLocalSize(x,&N);
2772:   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);
2773:   VecGetArrayWrite(x,&aa);

2775:   PetscMalloc1(m,a);
2776:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2777:   *a -= mstart;
2778:   return(0);
2779: }

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

2784:    Logically Collective

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

2794:    Level: developer

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

2802:    This routine actually zeros out the a pointer.

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

2817:   dummy = (void*)(*a + mstart);
2818:   PetscFree(dummy);
2819:   VecRestoreArray(x,NULL);
2820:   return(0);
2821: }

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

2826:    Logically Collective

2828:    Input Parameters:
2829: +  x - the vector
2830: .  m - first dimension of two dimensional array
2831: .  n - second dimension of the two dimensional array
2832: .  mstart - first index you will use in first coordinate direction (often 0)
2833: .  nstart - first index in the second coordinate direction (often 0)
2834: -  a - location of pointer to array obtained from VecGetArray2d()

2836:    Level: developer

2838:    Notes:
2839:    For regular PETSc vectors this routine does not involve any copies. For
2840:    any special vectors that do not store local vector data in a contiguous
2841:    array, this routine will copy the data back into the underlying
2842:    vector data structure from the array obtained with VecGetArray().

2844:    This routine actually zeros out the a pointer.

2846: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2847:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2848:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2849: @*/
2850: PetscErrorCode  VecRestoreArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2851: {
2853:   void           *dummy;

2859:   dummy = (void*)(*a + mstart);
2860:   PetscFree(dummy);
2861:   VecRestoreArrayWrite(x,NULL);
2862:   return(0);
2863: }

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

2870:    Logically Collective

2872:    Input Parameter:
2873: +  x - the vector
2874: .  m - first dimension of two dimensional array
2875: -  mstart - first index you will use in first coordinate direction (often 0)

2877:    Output Parameter:
2878: .  a - location to put pointer to the array

2880:    Level: developer

2882:   Notes:
2883:    For a vector obtained from DMCreateLocalVector() mstart are likely
2884:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2885:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

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

2902:   VecGetLocalSize(x,&N);
2903:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2904:   VecGetArray(x,a);
2905:   *a  -= mstart;
2906:   return(0);
2907: }

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

2914:    Logically Collective

2916:    Input Parameter:
2917: +  x - the vector
2918: .  m - first dimension of two dimensional array
2919: -  mstart - first index you will use in first coordinate direction (often 0)

2921:    Output Parameter:
2922: .  a - location to put pointer to the array

2924:    Level: developer

2926:   Notes:
2927:    For a vector obtained from DMCreateLocalVector() mstart are likely
2928:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2929:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2933: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2934:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2935:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2936: @*/
2937: PetscErrorCode  VecGetArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2938: {
2940:   PetscInt       N;

2946:   VecGetLocalSize(x,&N);
2947:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2948:   VecGetArrayWrite(x,a);
2949:   *a  -= mstart;
2950:   return(0);
2951: }

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

2956:    Logically Collective

2958:    Input Parameters:
2959: +  x - the vector
2960: .  m - first dimension of two dimensional array
2961: .  mstart - first index you will use in first coordinate direction (often 0)
2962: -  a - location of pointer to array obtained from VecGetArray21()

2964:    Level: developer

2966:    Notes:
2967:    For regular PETSc vectors this routine does not involve any copies. For
2968:    any special vectors that do not store local vector data in a contiguous
2969:    array, this routine will copy the data back into the underlying
2970:    vector data structure from the array obtained with VecGetArray1d().

2972:    This routine actually zeros out the a pointer.

2974: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2975:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2976:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2977: @*/
2978: PetscErrorCode  VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2979: {

2985:   VecRestoreArray(x,NULL);
2986:   return(0);
2987: }

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

2992:    Logically Collective

2994:    Input Parameters:
2995: +  x - the vector
2996: .  m - first dimension of two dimensional array
2997: .  mstart - first index you will use in first coordinate direction (often 0)
2998: -  a - location of pointer to array obtained from VecGetArray21()

3000:    Level: developer

3002:    Notes:
3003:    For regular PETSc vectors this routine does not involve any copies. For
3004:    any special vectors that do not store local vector data in a contiguous
3005:    array, this routine will copy the data back into the underlying
3006:    vector data structure from the array obtained with VecGetArray1d().

3008:    This routine actually zeros out the a pointer.

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

3012: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3013:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3014:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3015: @*/
3016: PetscErrorCode  VecRestoreArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3017: {

3023:   VecRestoreArrayWrite(x,NULL);
3024:   return(0);
3025: }

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

3032:    Logically Collective

3034:    Input Parameter:
3035: +  x - the vector
3036: .  m - first dimension of three dimensional array
3037: .  n - second dimension of three dimensional array
3038: .  p - third dimension of three dimensional array
3039: .  mstart - first index you will use in first coordinate direction (often 0)
3040: .  nstart - first index in the second coordinate direction (often 0)
3041: -  pstart - first index in the third coordinate direction (often 0)

3043:    Output Parameter:
3044: .  a - location to put pointer to the array

3046:    Level: developer

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

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

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

3070:   VecGetLocalSize(x,&N);
3071:   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);
3072:   VecGetArray(x,&aa);

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

3081:   *a -= mstart;
3082:   return(0);
3083: }

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

3090:    Logically Collective

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

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

3104:    Level: developer

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

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

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

3116: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3117:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3118:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3119: @*/
3120: PetscErrorCode  VecGetArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3121: {
3123:   PetscInt       i,N,j;
3124:   PetscScalar    *aa,**b;

3130:   VecGetLocalSize(x,&N);
3131:   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);
3132:   VecGetArrayWrite(x,&aa);

3134:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3135:   b    = (PetscScalar**)((*a) + m);
3136:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3137:   for (i=0; i<m; i++)
3138:     for (j=0; j<n; j++)
3139:       b[i*n+j] = aa + i*n*p + j*p - pstart;

3141:   *a -= mstart;
3142:   return(0);
3143: }

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

3148:    Logically Collective

3150:    Input Parameters:
3151: +  x - the vector
3152: .  m - first dimension of three dimensional array
3153: .  n - second dimension of the three dimensional array
3154: .  p - third dimension of the three dimensional array
3155: .  mstart - first index you will use in first coordinate direction (often 0)
3156: .  nstart - first index in the second coordinate direction (often 0)
3157: .  pstart - first index in the third coordinate direction (often 0)
3158: -  a - location of pointer to array obtained from VecGetArray3d()

3160:    Level: developer

3162:    Notes:
3163:    For regular PETSc vectors this routine does not involve any copies. For
3164:    any special vectors that do not store local vector data in a contiguous
3165:    array, this routine will copy the data back into the underlying
3166:    vector data structure from the array obtained with VecGetArray().

3168:    This routine actually zeros out the a pointer.

3170: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3171:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3172:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3173: @*/
3174: PetscErrorCode  VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3175: {
3177:   void           *dummy;

3183:   dummy = (void*)(*a + mstart);
3184:   PetscFree(dummy);
3185:   VecRestoreArray(x,NULL);
3186:   return(0);
3187: }

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

3192:    Logically Collective

3194:    Input Parameters:
3195: +  x - the vector
3196: .  m - first dimension of three dimensional array
3197: .  n - second dimension of the three dimensional array
3198: .  p - third dimension of the three dimensional array
3199: .  mstart - first index you will use in first coordinate direction (often 0)
3200: .  nstart - first index in the second coordinate direction (often 0)
3201: .  pstart - first index in the third coordinate direction (often 0)
3202: -  a - location of pointer to array obtained from VecGetArray3d()

3204:    Level: developer

3206:    Notes:
3207:    For regular PETSc vectors this routine does not involve any copies. For
3208:    any special vectors that do not store local vector data in a contiguous
3209:    array, this routine will copy the data back into the underlying
3210:    vector data structure from the array obtained with VecGetArray().

3212:    This routine actually zeros out the a pointer.

3214: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3215:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3216:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3217: @*/
3218: PetscErrorCode  VecRestoreArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3219: {
3221:   void           *dummy;

3227:   dummy = (void*)(*a + mstart);
3228:   PetscFree(dummy);
3229:   VecRestoreArrayWrite(x,NULL);
3230:   return(0);
3231: }

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

3238:    Logically Collective

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

3251:    Output Parameter:
3252: .  a - location to put pointer to the array

3254:    Level: beginner

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

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

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

3278:   VecGetLocalSize(x,&N);
3279:   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);
3280:   VecGetArray(x,&aa);

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

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

3302:    Logically Collective

3304:    Input Parameter:
3305: +  x - the vector
3306: .  m - first dimension of four dimensional array
3307: .  n - second dimension of four dimensional array
3308: .  p - third dimension of four dimensional array
3309: .  q - fourth dimension of four dimensional array
3310: .  mstart - first index you will use in first coordinate direction (often 0)
3311: .  nstart - first index in the second coordinate direction (often 0)
3312: .  pstart - first index in the third coordinate direction (often 0)
3313: -  qstart - first index in the fourth coordinate direction (often 0)

3315:    Output Parameter:
3316: .  a - location to put pointer to the array

3318:    Level: beginner

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

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

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

3330: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3331:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3332:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3333: @*/
3334: PetscErrorCode  VecGetArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3335: {
3337:   PetscInt       i,N,j,k;
3338:   PetscScalar    *aa,***b,**c;

3344:   VecGetLocalSize(x,&N);
3345:   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);
3346:   VecGetArrayWrite(x,&aa);

3348:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3349:   b    = (PetscScalar***)((*a) + m);
3350:   c    = (PetscScalar**)(b + m*n);
3351:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3352:   for (i=0; i<m; i++)
3353:     for (j=0; j<n; j++)
3354:       b[i*n+j] = c + i*n*p + j*p - pstart;
3355:   for (i=0; i<m; i++)
3356:     for (j=0; j<n; j++)
3357:       for (k=0; k<p; k++)
3358:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3359:   *a -= mstart;
3360:   return(0);
3361: }

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

3366:    Logically Collective

3368:    Input Parameters:
3369: +  x - the vector
3370: .  m - first dimension of four dimensional array
3371: .  n - second dimension of the four dimensional array
3372: .  p - third dimension of the four dimensional array
3373: .  q - fourth dimension of the four dimensional array
3374: .  mstart - first index you will use in first coordinate direction (often 0)
3375: .  nstart - first index in the second coordinate direction (often 0)
3376: .  pstart - first index in the third coordinate direction (often 0)
3377: .  qstart - first index in the fourth coordinate direction (often 0)
3378: -  a - location of pointer to array obtained from VecGetArray4d()

3380:    Level: beginner

3382:    Notes:
3383:    For regular PETSc vectors this routine does not involve any copies. For
3384:    any special vectors that do not store local vector data in a contiguous
3385:    array, this routine will copy the data back into the underlying
3386:    vector data structure from the array obtained with VecGetArray().

3388:    This routine actually zeros out the a pointer.

3390: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3391:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3392:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3393: @*/
3394: PetscErrorCode  VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3395: {
3397:   void           *dummy;

3403:   dummy = (void*)(*a + mstart);
3404:   PetscFree(dummy);
3405:   VecRestoreArray(x,NULL);
3406:   return(0);
3407: }

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

3412:    Logically Collective

3414:    Input Parameters:
3415: +  x - the vector
3416: .  m - first dimension of four dimensional array
3417: .  n - second dimension of the four dimensional array
3418: .  p - third dimension of the four dimensional array
3419: .  q - fourth dimension of the four dimensional array
3420: .  mstart - first index you will use in first coordinate direction (often 0)
3421: .  nstart - first index in the second coordinate direction (often 0)
3422: .  pstart - first index in the third coordinate direction (often 0)
3423: .  qstart - first index in the fourth coordinate direction (often 0)
3424: -  a - location of pointer to array obtained from VecGetArray4d()

3426:    Level: beginner

3428:    Notes:
3429:    For regular PETSc vectors this routine does not involve any copies. For
3430:    any special vectors that do not store local vector data in a contiguous
3431:    array, this routine will copy the data back into the underlying
3432:    vector data structure from the array obtained with VecGetArray().

3434:    This routine actually zeros out the a pointer.

3436: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3437:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3438:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3439: @*/
3440: PetscErrorCode  VecRestoreArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3441: {
3443:   void           *dummy;

3449:   dummy = (void*)(*a + mstart);
3450:   PetscFree(dummy);
3451:   VecRestoreArrayWrite(x,NULL);
3452:   return(0);
3453: }

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

3460:    Logically Collective

3462:    Input Parameter:
3463: +  x - the vector
3464: .  m - first dimension of two dimensional array
3465: .  n - second dimension of two dimensional array
3466: .  mstart - first index you will use in first coordinate direction (often 0)
3467: -  nstart - first index in the second coordinate direction (often 0)

3469:    Output Parameter:
3470: .  a - location to put pointer to the array

3472:    Level: developer

3474:   Notes:
3475:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
3476:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3477:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3478:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

3482: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3483:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3484:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3485: @*/
3486: PetscErrorCode  VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3487: {
3488:   PetscErrorCode    ierr;
3489:   PetscInt          i,N;
3490:   const PetscScalar *aa;

3496:   VecGetLocalSize(x,&N);
3497:   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);
3498:   VecGetArrayRead(x,&aa);

3500:   PetscMalloc1(m,a);
3501:   for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
3502:   *a -= mstart;
3503:   return(0);
3504: }

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

3509:    Logically Collective

3511:    Input Parameters:
3512: +  x - the vector
3513: .  m - first dimension of two dimensional array
3514: .  n - second dimension of the two dimensional array
3515: .  mstart - first index you will use in first coordinate direction (often 0)
3516: .  nstart - first index in the second coordinate direction (often 0)
3517: -  a - location of pointer to array obtained from VecGetArray2d()

3519:    Level: developer

3521:    Notes:
3522:    For regular PETSc vectors this routine does not involve any copies. For
3523:    any special vectors that do not store local vector data in a contiguous
3524:    array, this routine will copy the data back into the underlying
3525:    vector data structure from the array obtained with VecGetArray().

3527:    This routine actually zeros out the a pointer.

3529: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3530:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3531:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3532: @*/
3533: PetscErrorCode  VecRestoreArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3534: {
3536:   void           *dummy;

3542:   dummy = (void*)(*a + mstart);
3543:   PetscFree(dummy);
3544:   VecRestoreArrayRead(x,NULL);
3545:   return(0);
3546: }

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

3553:    Logically Collective

3555:    Input Parameter:
3556: +  x - the vector
3557: .  m - first dimension of two dimensional array
3558: -  mstart - first index you will use in first coordinate direction (often 0)

3560:    Output Parameter:
3561: .  a - location to put pointer to the array

3563:    Level: developer

3565:   Notes:
3566:    For a vector obtained from DMCreateLocalVector() mstart are likely
3567:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3568:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3572: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3573:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3574:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3575: @*/
3576: PetscErrorCode  VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3577: {
3579:   PetscInt       N;

3585:   VecGetLocalSize(x,&N);
3586:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3587:   VecGetArrayRead(x,(const PetscScalar**)a);
3588:   *a  -= mstart;
3589:   return(0);
3590: }

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

3595:    Logically Collective

3597:    Input Parameters:
3598: +  x - the vector
3599: .  m - first dimension of two dimensional array
3600: .  mstart - first index you will use in first coordinate direction (often 0)
3601: -  a - location of pointer to array obtained from VecGetArray21()

3603:    Level: developer

3605:    Notes:
3606:    For regular PETSc vectors this routine does not involve any copies. For
3607:    any special vectors that do not store local vector data in a contiguous
3608:    array, this routine will copy the data back into the underlying
3609:    vector data structure from the array obtained with VecGetArray1dRead().

3611:    This routine actually zeros out the a pointer.

3613: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3614:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3615:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3616: @*/
3617: PetscErrorCode  VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3618: {

3624:   VecRestoreArrayRead(x,NULL);
3625:   return(0);
3626: }


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

3634:    Logically Collective

3636:    Input Parameter:
3637: +  x - the vector
3638: .  m - first dimension of three dimensional array
3639: .  n - second dimension of three dimensional array
3640: .  p - third dimension of three dimensional array
3641: .  mstart - first index you will use in first coordinate direction (often 0)
3642: .  nstart - first index in the second coordinate direction (often 0)
3643: -  pstart - first index in the third coordinate direction (often 0)

3645:    Output Parameter:
3646: .  a - location to put pointer to the array

3648:    Level: developer

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

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

3658: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3659:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3660:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3661: @*/
3662: PetscErrorCode  VecGetArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3663: {
3664:   PetscErrorCode    ierr;
3665:   PetscInt          i,N,j;
3666:   const PetscScalar *aa;
3667:   PetscScalar       **b;

3673:   VecGetLocalSize(x,&N);
3674:   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);
3675:   VecGetArrayRead(x,&aa);

3677:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3678:   b    = (PetscScalar**)((*a) + m);
3679:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3680:   for (i=0; i<m; i++)
3681:     for (j=0; j<n; j++)
3682:       b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;

3684:   *a -= mstart;
3685:   return(0);
3686: }

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

3691:    Logically Collective

3693:    Input Parameters:
3694: +  x - the vector
3695: .  m - first dimension of three dimensional array
3696: .  n - second dimension of the three dimensional array
3697: .  p - third dimension of the three dimensional array
3698: .  mstart - first index you will use in first coordinate direction (often 0)
3699: .  nstart - first index in the second coordinate direction (often 0)
3700: .  pstart - first index in the third coordinate direction (often 0)
3701: -  a - location of pointer to array obtained from VecGetArray3dRead()

3703:    Level: developer

3705:    Notes:
3706:    For regular PETSc vectors this routine does not involve any copies. For
3707:    any special vectors that do not store local vector data in a contiguous
3708:    array, this routine will copy the data back into the underlying
3709:    vector data structure from the array obtained with VecGetArray().

3711:    This routine actually zeros out the a pointer.

3713: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3714:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3715:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3716: @*/
3717: PetscErrorCode  VecRestoreArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3718: {
3720:   void           *dummy;

3726:   dummy = (void*)(*a + mstart);
3727:   PetscFree(dummy);
3728:   VecRestoreArrayRead(x,NULL);
3729:   return(0);
3730: }

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

3737:    Logically Collective

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

3750:    Output Parameter:
3751: .  a - location to put pointer to the array

3753:    Level: beginner

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

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

3763: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3764:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3765:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3766: @*/
3767: PetscErrorCode  VecGetArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3768: {
3769:   PetscErrorCode    ierr;
3770:   PetscInt          i,N,j,k;
3771:   const PetscScalar *aa;
3772:   PetscScalar       ***b,**c;

3778:   VecGetLocalSize(x,&N);
3779:   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);
3780:   VecGetArrayRead(x,&aa);

3782:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3783:   b    = (PetscScalar***)((*a) + m);
3784:   c    = (PetscScalar**)(b + m*n);
3785:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3786:   for (i=0; i<m; i++)
3787:     for (j=0; j<n; j++)
3788:       b[i*n+j] = c + i*n*p + j*p - pstart;
3789:   for (i=0; i<m; i++)
3790:     for (j=0; j<n; j++)
3791:       for (k=0; k<p; k++)
3792:         c[i*n*p+j*p+k] = (PetscScalar*) aa + i*n*p*q + j*p*q + k*q - qstart;
3793:   *a -= mstart;
3794:   return(0);
3795: }

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

3800:    Logically Collective

3802:    Input Parameters:
3803: +  x - the vector
3804: .  m - first dimension of four dimensional array
3805: .  n - second dimension of the four dimensional array
3806: .  p - third dimension of the four dimensional array
3807: .  q - fourth dimension of the four dimensional array
3808: .  mstart - first index you will use in first coordinate direction (often 0)
3809: .  nstart - first index in the second coordinate direction (often 0)
3810: .  pstart - first index in the third coordinate direction (often 0)
3811: .  qstart - first index in the fourth coordinate direction (often 0)
3812: -  a - location of pointer to array obtained from VecGetArray4dRead()

3814:    Level: beginner

3816:    Notes:
3817:    For regular PETSc vectors this routine does not involve any copies. For
3818:    any special vectors that do not store local vector data in a contiguous
3819:    array, this routine will copy the data back into the underlying
3820:    vector data structure from the array obtained with VecGetArray().

3822:    This routine actually zeros out the a pointer.

3824: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3825:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3826:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3827: @*/
3828: PetscErrorCode  VecRestoreArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3829: {
3831:   void           *dummy;

3837:   dummy = (void*)(*a + mstart);
3838:   PetscFree(dummy);
3839:   VecRestoreArrayRead(x,NULL);
3840:   return(0);
3841: }

3843: #if defined(PETSC_USE_DEBUG)

3845: /*@
3846:    VecLockGet  - Gets the current lock status of a vector

3848:    Logically Collective on Vec

3850:    Input Parameter:
3851: .  x - the vector

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

3857:    Level: beginner

3859: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop()
3860: @*/
3861: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
3862: {
3865:   *state = x->lock;
3866:   return(0);
3867: }

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

3872:    Logically Collective on Vec

3874:    Input Parameter:
3875: .  x - the vector

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

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

3883:    Level: beginner

3885: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPop(), VecLockGet()
3886: @*/
3887: PetscErrorCode VecLockReadPush(Vec x)
3888: {
3891:   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");
3892:   x->lock++;
3893:   return(0);
3894: }

3896: /*@
3897:    VecLockReadPop  - Pops a read-only lock from a vector

3899:    Logically Collective on Vec

3901:    Input Parameter:
3902: .  x - the vector

3904:    Level: beginner

3906: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockGet()
3907: @*/
3908: PetscErrorCode VecLockReadPop(Vec x)
3909: {
3912:   x->lock--;
3913:   if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector has been unlocked from read-only access too many times");
3914:   return(0);
3915: }

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

3920:    Logically Collective on Vec

3922:    Input Parameter:
3923: +  x   - the vector
3924: -  flg - PETSC_TRUE to lock the vector for writing; PETSC_FALSE to unlock it.

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


3933:        VecGetArray(x,&xdata); // begin phase
3934:        VecLockWriteSet_Private(v,PETSC_TRUE);

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

3938:        VecRestoreArray(x,&vdata); // end phase
3939:        VecLockWriteSet_Private(v,PETSC_FALSE);

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

3944:    Level: beginner

3946: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop(), VecLockGet()
3947: @*/
3948: PetscErrorCode VecLockWriteSet_Private(Vec x,PetscBool flg)
3949: {
3952:   if (flg) {
3953:     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");
3954:     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");
3955:     else x->lock = -1;
3956:   } else {
3957:     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");
3958:     x->lock = 0;
3959:   }
3960:   return(0);
3961: }

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

3966:    Level: deprecated

3968: .seealso: VecLockReadPush()
3969: @*/
3970: PetscErrorCode VecLockPush(Vec x)
3971: {
3974:   VecLockReadPush(x);
3975:   return(0);
3976: }

3978: /*@
3979:    VecLockPop  - Pops a read-only lock from a vector

3981:    Level: deprecated

3983: .seealso: VecLockReadPop()
3984: @*/
3985: PetscErrorCode VecLockPop(Vec x)
3986: {
3989:   VecLockReadPop(x);
3990:   return(0);
3991: }

3993: #endif