Actual source code: rvector.c

petsc-3.14.0 2020-09-29
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/vecimpl.h>
  6: #if defined(PETSC_HAVE_CUDA)
  7: #include <../src/vec/vec/impls/dvecimpl.h>
  8: #include <petsc/private/cudavecimpl.h>
  9: #endif
 10: static PetscInt VecGetSubVectorSavedStateId = -1;

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

 20: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL)
 21:   if ((vec->petscnative || vec->ops->getarray) && (vec->offloadmask & PETSC_OFFLOAD_CPU)) {
 22: #else
 23:   if (vec->petscnative || vec->ops->getarray) {
 24: #endif
 25:     VecGetLocalSize(vec,&n);
 26:     VecGetArrayRead(vec,&x);
 27:     for (i=0; i<n; i++) {
 28:       if (begin) {
 29:         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);
 30:       } else {
 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 end of function: Parameter number %D",i,argnum);
 32:       }
 33:     }
 34:     VecRestoreArrayRead(vec,&x);
 35:   }
 36: #else
 38: #endif
 39:   return(0);
 40: }

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

 45:    Logically Collective on Vec

 47:    Input Parameters:
 48: .  x, y  - the vectors

 50:    Output Parameter:
 51: .  max - the result

 53:    Level: advanced

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

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

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

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

 80:    Collective on Vec

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

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

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

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

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

104:    Level: intermediate


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

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

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

128: /*@
129:    VecDotRealPart - Computes the real part of the vector dot product.

131:    Collective on Vec

133:    Input Parameters:
134: .  x, y - the vectors

136:    Output Parameter:
137: .  val - the real part of the dot product;

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

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

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

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

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

154:    Level: intermediate


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

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

170: /*@
171:    VecNorm  - Computes the vector norm.

173:    Collective on Vec

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

181:    Output Parameter:
182: .  val - the norm

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

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

194:    Level: intermediate

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


202: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNormAvailable(),
203:           VecNormBegin(), VecNormEnd()

205: @*/

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


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

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

236:    Not Collective

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

244:    Output Parameter:
245: +  available - PETSC_TRUE if the val returned is valid
246: -  val - the norm

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

253:    Level: intermediate

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

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


266: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNorm()
267:           VecNormBegin(), VecNormEnd()

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


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

286: /*@
287:    VecNormalize - Normalizes a vector by 2-norm.

289:    Collective on Vec

291:    Input Parameters:
292: +  x - the vector

294:    Output Parameter:
295: .  x - the normalized vector
296: -  val - the vector norm before normalization

298:    Level: intermediate


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

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

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

326:    Collective on Vec

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

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

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

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


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

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

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

361:    Collective on Vec

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

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

370:    Level: intermediate

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

375:    This returns the smallest index with the minumum value


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

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

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

398:    Collective on Vec

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

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

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

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

415:    Level: intermediate

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

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

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

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

441:    Not collective on Vec

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

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

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

454:    Level: intermediate


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

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

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

492:    Logically Collective on Vec

494:    Input Parameters:
495: +  x  - the vector
496: -  alpha - the scalar

498:    Output Parameter:
499: .  x  - the vector

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

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

511:    Level: beginner

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

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

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

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

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


553: /*@
554:    VecAXPY - Computes y = alpha x + y.

556:    Logically Collective on Vec

558:    Input Parameters:
559: +  alpha - the scalar
560: -  x, y  - the vectors

562:    Output Parameter:
563: .  y - output vector

565:    Level: intermediate

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

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


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

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

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

606: /*@
607:    VecAXPBY - Computes y = alpha x + beta y.

609:    Logically Collective on Vec

611:    Input Parameters:
612: +  alpha,beta - the scalars
613: -  x, y  - the vectors

615:    Output Parameter:
616: .  y - output vector

618:    Level: intermediate

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


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

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

650: /*@
651:    VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z

653:    Logically Collective on Vec

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

659:    Output Parameter:
660: .  z - output vector

662:    Level: intermediate

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


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

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

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

701: /*@
702:    VecAYPX - Computes y = x + beta y.

704:    Logically Collective on Vec

706:    Input Parameters:
707: +  beta - the scalar
708: -  x, y  - the vectors

710:    Output Parameter:
711: .  y - output vector

713:    Level: intermediate

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


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

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

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


745: /*@
746:    VecWAXPY - Computes w = alpha x + y.

748:    Logically Collective on Vec

750:    Input Parameters:
751: +  alpha - the scalar
752: -  x, y  - the vectors

754:    Output Parameter:
755: .  w - the result

757:    Level: intermediate

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


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

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

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


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

797:    Not Collective

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

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

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

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

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

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

826:    Level: beginner

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

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

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

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

853:     Not Collective

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

860:    Output Parameter:
861: .   y - array of values

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

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

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

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

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

876:    Level: beginner

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

886:   if (!ni) return(0);
890:   (*x->ops->getvalues)(x,ni,ix,y);
891:   return(0);
892: }

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

897:    Not Collective

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

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

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

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

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

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

926:    Level: intermediate

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

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

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


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

954:    Not Collective

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

965:    Level: intermediate

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

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

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

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

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

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

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

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

1017:    Not Collective

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

1028:    Level: intermediate

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

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

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

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


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

1054:   if (!ni) return(0);
1058:   if (ni > 128) {
1059:     PetscMalloc1(ni,&lix);
1060:   }

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

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

1077:    Collective on Vec

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

1084:    Output Parameter:
1085: .  val - array of the dot products

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

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

1096:    Level: intermediate


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

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

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

1123: /*@
1124:    VecMDot - Computes vector multiple dot products.

1126:    Collective on Vec

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

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

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

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

1145:    Level: intermediate


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

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

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

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

1176:    Logically Collective on Vec

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

1184:    Level: intermediate

1186:    Notes:
1187:     y cannot be any of the x vectors

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

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

1220: /*@
1221:    VecGetSubVector - Gets a vector representing part of another vector

1223:    Collective on IS

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

1229:    Output Arguments:
1230: . Y - subvector corresponding to is

1232:    Level: advanced

1234:    Notes:
1235:    The subvector Y should be returned with VecRestoreSubVector().

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

1240: .seealso: MatCreateSubMatrix()
1241: @*/
1242: PetscErrorCode  VecGetSubVector(Vec X,IS is,Vec *Y)
1243: {
1244:   PetscErrorCode   ierr;
1245:   Vec              Z;

1251:   if (X->ops->getsubvector) {
1252:     (*X->ops->getsubvector)(X,is,&Z);
1253:   } else { /* Default implementation currently does no caching */
1254:     PetscInt  gstart,gend,start;
1255:     PetscBool contiguous,gcontiguous;

1257:     VecGetOwnershipRange(X,&gstart,&gend);
1258:     ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1259:     MPIU_Allreduce(&contiguous,&gcontiguous,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1260:     if (gcontiguous) { /* We can do a no-copy implementation */
1261:       const PetscScalar *x;
1262:       PetscInt          n,N,bs;
1263:       PetscInt          state = 0;
1264: #if defined(PETSC_HAVE_CUDA)
1265:       PetscBool         iscuda;
1266: #endif

1268:       ISGetSize(is,&N);
1269:       ISGetLocalSize(is,&n);
1270:       VecGetBlockSize(X,&bs);
1271:       if (n%bs || bs == 1) bs = -1; /* Do not decide block size if we do not have to */
1272: #if defined(PETSC_HAVE_CUDA)
1273:       PetscObjectTypeCompareAny((PetscObject)X,&iscuda,VECSEQCUDA,VECMPICUDA,"");
1274:       if (iscuda) {
1275:         const PetscScalar *x_d;
1276:         PetscMPIInt       size;
1277:         PetscOffloadMask  flg;

1279:         VecCUDAGetArrays_Private(X,&x,&x_d,&flg);
1280:         if (flg == PETSC_OFFLOAD_UNALLOCATED) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for PETSC_OFFLOAD_UNALLOCATED");
1281:         if (n && !x && !x_d) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Missing vector data");
1282:         if (x) x += start;
1283:         if (x_d) x_d += start;
1284:         MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1285:         if (size == 1) {
1286:           VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X),bs,n,x,x_d,&Z);
1287:         } else {
1288:           VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X),bs,n,N,x,x_d,&Z);
1289:         }
1290:         Z->offloadmask = flg;
1291:       } else {
1292: #else
1293:       {
1294: #endif
1295:         VecGetArrayRead(X,&x);
1296:         VecCreate(PetscObjectComm((PetscObject)X),&Z);
1297:         VecSetType(Z,((PetscObject)X)->type_name);
1298:         VecSetSizes(Z,n,N);
1299:         VecSetBlockSize(Z,bs);
1300:         VecPlaceArray(Z,x+start);
1301:         VecRestoreArrayRead(X,&x);
1302:       }

1304:       /* this is relevant only in debug mode */
1305:       VecLockGet(X,&state);
1306:       if (state) {
1307:         VecLockReadPush(Z);
1308:       }
1309:       Z->ops->placearray = NULL;
1310:       Z->ops->replacearray = NULL;
1311:     } else { /* Have to create a scatter and do a copy */
1312:       VecScatter scatter;
1313:       PetscInt   n,N;

1315:       ISGetLocalSize(is,&n);
1316:       ISGetSize(is,&N);
1317:       VecCreate(PetscObjectComm((PetscObject)is),&Z);
1318:       VecSetSizes(Z,n,N);
1319:       VecSetType(Z,((PetscObject)X)->type_name);
1320:       VecScatterCreate(X,is,Z,NULL,&scatter);
1321:       VecScatterBegin(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1322:       VecScatterEnd(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1323:       PetscObjectCompose((PetscObject)Z,"VecGetSubVector_Scatter",(PetscObject)scatter);
1324:       VecScatterDestroy(&scatter);
1325:     }
1326:   }
1327:   /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1328:   if (VecGetSubVectorSavedStateId < 0) {PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);}
1329:   PetscObjectComposedDataSetInt((PetscObject)Z,VecGetSubVectorSavedStateId,1);
1330:   *Y   = Z;
1331:   return(0);
1332: }

1334: /*@
1335:    VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()

1337:    Collective on IS

1339:    Input Arguments:
1340: + X - vector from which subvector was obtained
1341: . is - index set representing the subset of X
1342: - Y - subvector being restored

1344:    Level: advanced

1346: .seealso: VecGetSubVector()
1347: @*/
1348: PetscErrorCode  VecRestoreSubVector(Vec X,IS is,Vec *Y)
1349: {

1357:   if (X->ops->restoresubvector) {
1358:     (*X->ops->restoresubvector)(X,is,Y);
1359:   } else {
1360:     PETSC_UNUSED PetscObjectState dummystate = 0;
1361:     PetscBool valid;

1363:     PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,dummystate,valid);
1364:     if (!valid) {
1365:       VecScatter scatter;
1366:       PetscInt   state;

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

1371:       PetscObjectQuery((PetscObject)*Y,"VecGetSubVector_Scatter",(PetscObject*)&scatter);
1372:       if (scatter) {
1373:         VecScatterBegin(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1374:         VecScatterEnd(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1375:       } else {
1376: #if defined(PETSC_HAVE_CUDA)
1377:         PetscBool iscuda;

1379:         PetscObjectTypeCompareAny((PetscObject)*Y,&iscuda,VECSEQCUDA,VECMPICUDA,"");
1380:         if (iscuda) {
1381:           PetscOffloadMask ymask = (*Y)->offloadmask;

1383:           /* The offloadmask of X dictates where to move memory
1384:              If X GPU data is valid, then move Y data on GPU if needed
1385:              Otherwise, move back to the CPU */
1386:           switch (X->offloadmask) {
1387:           case PETSC_OFFLOAD_BOTH:
1388:             if (ymask == PETSC_OFFLOAD_CPU) {
1389:               VecCUDAResetArray(*Y);
1390:             } else if (ymask == PETSC_OFFLOAD_GPU) {
1391:               X->offloadmask = PETSC_OFFLOAD_GPU;
1392:             }
1393:             break;
1394:           case PETSC_OFFLOAD_GPU:
1395:             if (ymask == PETSC_OFFLOAD_CPU) {
1396:               VecCUDAResetArray(*Y);
1397:             }
1398:             break;
1399:           case PETSC_OFFLOAD_CPU:
1400:             if (ymask == PETSC_OFFLOAD_GPU) {
1401:               VecResetArray(*Y);
1402:             }
1403:             break;
1404:           case PETSC_OFFLOAD_UNALLOCATED:
1405:           case PETSC_OFFLOAD_VECKOKKOS:
1406:             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"This should not happen");
1407:             break;
1408:           }
1409:         } else {
1410: #else
1411:         {
1412: #endif
1413:           /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1414:           VecResetArray(*Y);
1415:         }
1416:         PetscObjectStateIncrease((PetscObject)X);
1417:       }
1418:     }
1419:     VecDestroy(Y);
1420:   }
1421:   return(0);
1422: }

1424: /*@
1425:    VecGetLocalVectorRead - Maps the local portion of a vector into a
1426:    vector.  You must call VecRestoreLocalVectorRead() when the local
1427:    vector is no longer needed.

1429:    Not collective.

1431:    Input parameter:
1432: .  v - The vector for which the local vector is desired.

1434:    Output parameter:
1435: .  w - Upon exit this contains the local vector.

1437:    Level: beginner

1439:    Notes:
1440:    This function is similar to VecGetArrayRead() which maps the local
1441:    portion into a raw pointer.  VecGetLocalVectorRead() is usually
1442:    almost as efficient as VecGetArrayRead() but in certain circumstances
1443:    VecGetLocalVectorRead() can be much more efficient than
1444:    VecGetArrayRead().  This is because the construction of a contiguous
1445:    array representing the vector data required by VecGetArrayRead() can
1446:    be an expensive operation for certain vector types.  For example, for
1447:    GPU vectors VecGetArrayRead() requires that the data between device
1448:    and host is synchronized.

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

1453: .seealso: VecRestoreLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1454: @*/
1455: PetscErrorCode VecGetLocalVectorRead(Vec v,Vec w)
1456: {
1458:   PetscScalar    *a;

1463:   VecCheckSameLocalSize(v,1,w,2);
1464:   if (v->ops->getlocalvectorread) {
1465:     (*v->ops->getlocalvectorread)(v,w);
1466:   } else {
1467:     VecGetArrayRead(v,(const PetscScalar**)&a);
1468:     VecPlaceArray(w,a);
1469:   }
1470:   return(0);
1471: }

1473: /*@
1474:    VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1475:    previously mapped into a vector using VecGetLocalVectorRead().

1477:    Not collective.

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

1483:    Level: beginner

1485: .seealso: VecGetLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1486: @*/
1487: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec w)
1488: {
1490:   PetscScalar    *a;

1495:   if (v->ops->restorelocalvectorread) {
1496:     (*v->ops->restorelocalvectorread)(v,w);
1497:   } else {
1498:     VecGetArrayRead(w,(const PetscScalar**)&a);
1499:     VecRestoreArrayRead(v,(const PetscScalar**)&a);
1500:     VecResetArray(w);
1501:   }
1502:   return(0);
1503: }

1505: /*@
1506:    VecGetLocalVector - Maps the local portion of a vector into a
1507:    vector.

1509:    Collective on v, not collective on w.

1511:    Input parameter:
1512: .  v - The vector for which the local vector is desired.

1514:    Output parameter:
1515: .  w - Upon exit this contains the local vector.

1517:    Level: beginner

1519:    Notes:
1520:    This function is similar to VecGetArray() which maps the local
1521:    portion into a raw pointer.  VecGetLocalVector() is usually about as
1522:    efficient as VecGetArray() but in certain circumstances
1523:    VecGetLocalVector() can be much more efficient than VecGetArray().
1524:    This is because the construction of a contiguous array representing
1525:    the vector data required by VecGetArray() can be an expensive
1526:    operation for certain vector types.  For example, for GPU vectors
1527:    VecGetArray() requires that the data between device and host is
1528:    synchronized.

1530: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1531: @*/
1532: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1533: {
1535:   PetscScalar    *a;

1540:   VecCheckSameLocalSize(v,1,w,2);
1541:   if (v->ops->getlocalvector) {
1542:     (*v->ops->getlocalvector)(v,w);
1543:   } else {
1544:     VecGetArray(v,&a);
1545:     VecPlaceArray(w,a);
1546:   }
1547:   return(0);
1548: }

1550: /*@
1551:    VecRestoreLocalVector - Unmaps the local portion of a vector
1552:    previously mapped into a vector using VecGetLocalVector().

1554:    Logically collective.

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

1560:    Level: beginner

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

1572:   if (v->ops->restorelocalvector) {
1573:     (*v->ops->restorelocalvector)(v,w);
1574:   } else {
1575:     VecGetArray(w,&a);
1576:     VecRestoreArray(v,&a);
1577:     VecResetArray(w);
1578:   }
1579:   return(0);
1580: }

1582: /*@C
1583:    VecGetArray - Returns a pointer to a contiguous array that contains this
1584:    processor's portion of the vector data. For the standard PETSc
1585:    vectors, VecGetArray() returns a pointer to the local data array and
1586:    does not use any copies. If the underlying vector data is not stored
1587:    in a contiguous array this routine will copy the data to a contiguous
1588:    array and return a pointer to that. You MUST call VecRestoreArray()
1589:    when you no longer need access to the array.

1591:    Logically Collective on Vec

1593:    Input Parameter:
1594: .  x - the vector

1596:    Output Parameter:
1597: .  a - location to put pointer to the array

1599:    Fortran Note:
1600:    This routine is used differently from Fortran 77
1601: $    Vec         x
1602: $    PetscScalar x_array(1)
1603: $    PetscOffset i_x
1604: $    PetscErrorCode ierr
1605: $       call VecGetArray(x,x_array,i_x,ierr)
1606: $
1607: $   Access first local entry in vector with
1608: $      value = x_array(i_x + 1)
1609: $
1610: $      ...... other code
1611: $       call VecRestoreArray(x,x_array,i_x,ierr)
1612:    For Fortran 90 see VecGetArrayF90()

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

1617:    Level: beginner

1619: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1620:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1621: @*/
1622: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1623: {
1625: #if defined(PETSC_HAVE_VIENNACL)
1626:   PetscBool      is_viennacltype = PETSC_FALSE;
1627: #endif

1631:   VecSetErrorIfLocked(x,1);
1632:   if (x->petscnative) {
1633: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
1634:     if (x->offloadmask == PETSC_OFFLOAD_VECKOKKOS) { /* offloadmask here works as a tag quickly saying this is a VecKokkos */
1635:       VecKokkosSyncHost(x);
1636:       *a   = *((PetscScalar**)x->data);
1637:       return(0);
1638:     }
1639: #endif
1640: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1641:     if (x->offloadmask == PETSC_OFFLOAD_GPU) {
1642:   #if defined(PETSC_HAVE_VIENNACL)
1643:       PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1644:       if (is_viennacltype) {
1645:         VecViennaCLCopyFromGPU(x);
1646:       } else
1647:   #endif
1648:       {
1649:   #if defined(PETSC_HAVE_CUDA)
1650:         VecCUDACopyFromGPU(x);
1651:   #endif
1652:       }
1653:     } else if (x->offloadmask == PETSC_OFFLOAD_UNALLOCATED) {
1654:   #if defined(PETSC_HAVE_VIENNACL)
1655:       PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1656:       if (is_viennacltype) {
1657:         VecViennaCLAllocateCheckHost(x);
1658:       } else
1659:   #endif
1660:       {
1661:   #if defined(PETSC_HAVE_CUDA)
1662:         VecCUDAAllocateCheckHost(x);
1663:   #endif
1664:       }
1665:     }
1666: #endif
1667:     *a = *((PetscScalar**)x->data);
1668:   } else {
1669:     if (x->ops->getarray) {
1670:       (*x->ops->getarray)(x,a);
1671:     } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array for vector type \"%s\"",((PetscObject)x)->type_name);
1672:   }
1673:   return(0);
1674: }

1676: /*@C
1677:    VecGetArrayInPlace - Like VecGetArray(), but if this is a CUDA vector and it is currently offloaded to GPU,
1678:    the returned pointer will be a GPU pointer to the GPU memory that contains this processor's portion of the
1679:    vector data. Otherwise, it functions as VecGetArray().

1681:    Logically Collective on Vec

1683:    Input Parameter:
1684: .  x - the vector

1686:    Output Parameter:
1687: .  a - location to put pointer to the array

1689:    Level: beginner

1691: .seealso: VecRestoreArrayInPlace(), VecRestoreArrayInPlace(), VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(),
1692:           VecPlaceArray(), VecGetArray2d(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1693: @*/
1694: PetscErrorCode VecGetArrayInPlace(Vec x,PetscScalar **a)
1695: {

1700:   VecSetErrorIfLocked(x,1);

1702: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
1703:   if (x->offloadmask == PETSC_OFFLOAD_VECKOKKOS) {
1704:     VecKokkosGetArrayInPlace(x,a);
1705:     return(0);
1706:   }
1707: #endif

1709: #if defined(PETSC_HAVE_CUDA)
1710:   if (x->petscnative && (x->offloadmask & PETSC_OFFLOAD_GPU)) { /* Prefer working on GPU when offloadmask is PETSC_OFFLOAD_BOTH */
1711:     PetscBool is_cudatype = PETSC_FALSE;
1712:     PetscObjectTypeCompareAny((PetscObject)x,&is_cudatype,VECSEQCUDA,VECMPICUDA,VECCUDA,"");
1713:     if (is_cudatype) {
1714:       VecCUDAGetArray(x,a);
1715:       x->offloadmask = PETSC_OFFLOAD_GPU; /* Change the mask once GPU gets write access, don't wait until restore array */
1716:       return(0);
1717:     }
1718:   }
1719: #endif
1720:   VecGetArray(x,a);
1721:   return(0);
1722: }

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

1729:    Logically Collective on Vec

1731:    Input Parameter:
1732: .  x - the vector

1734:    Output Parameter:
1735: .  a - location to put pointer to the array

1737:    Level: intermediate

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

1742:    Concepts: vector^accessing local values

1744: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1745:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArray(), VecRestoreArrayWrite()
1746: @*/
1747: PetscErrorCode VecGetArrayWrite(Vec x,PetscScalar **a)
1748: {

1753:   VecSetErrorIfLocked(x,1);
1754:   if (!x->ops->getarraywrite) {
1755:     VecGetArray(x,a);
1756:   } else {
1757:     (*x->ops->getarraywrite)(x,a);
1758:   }
1759:   return(0);
1760: }

1762: /*@C
1763:    VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.

1765:    Not Collective

1767:    Input Parameters:
1768: .  x - the vector

1770:    Output Parameter:
1771: .  a - the array

1773:    Level: beginner

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

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

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

1784: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1785: @*/
1786: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1787: {
1789: #if defined(PETSC_HAVE_VIENNACL)
1790:   PetscBool      is_viennacltype = PETSC_FALSE;
1791: #endif

1795:   if (x->petscnative) {
1796: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
1797:     if (x->offloadmask == PETSC_OFFLOAD_VECKOKKOS) {
1798:       VecKokkosSyncHost(x);
1799:       *a   = *((PetscScalar **)x->data);
1800:       return(0);
1801:     }
1802: #endif
1803: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1804:     if (x->offloadmask == PETSC_OFFLOAD_GPU) {
1805: #if defined(PETSC_HAVE_VIENNACL)
1806:       PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1807:       if (is_viennacltype) {
1808:         VecViennaCLCopyFromGPU(x);
1809:       } else
1810: #endif
1811:       {
1812: #if defined(PETSC_HAVE_CUDA)
1813:         VecCUDACopyFromGPU(x);
1814: #endif
1815:       }
1816:     }
1817: #endif
1818:     *a = *((PetscScalar **)x->data);
1819:   } else if (x->ops->getarrayread) {
1820:     (*x->ops->getarrayread)(x,a);
1821:   } else {
1822:     (*x->ops->getarray)(x,(PetscScalar**)a);
1823:   }
1824:   return(0);
1825: }

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

1832:    Not Collective

1834:    Input Parameters:
1835: .  x - the vector

1837:    Output Parameter:
1838: .  a - the array

1840:    Level: beginner

1842:    Notes:
1843:    The array must be returned using a matching call to VecRestoreArrayReadInPlace().


1846: .seealso: VecRestoreArrayReadInPlace(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayInPlace()
1847: @*/
1848: PetscErrorCode VecGetArrayReadInPlace(Vec x,const PetscScalar **a)
1849: {

1854: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
1855:   if (x->offloadmask == PETSC_OFFLOAD_VECKOKKOS) {
1856:     VecKokkosGetArrayReadInPlace(x,a);
1857:     return(0);
1858:   }
1859: #endif

1861: #if defined(PETSC_HAVE_CUDA)
1862:   if (x->petscnative && x->offloadmask & PETSC_OFFLOAD_GPU) {
1863:     PetscBool is_cudatype = PETSC_FALSE;
1864:     PetscObjectTypeCompareAny((PetscObject)x,&is_cudatype,VECSEQCUDA,VECMPICUDA,VECCUDA,"");
1865:     if (is_cudatype) {
1866:       VecCUDAGetArrayRead(x,a);
1867:       return(0);
1868:     }
1869:   }
1870: #endif
1871:   VecGetArrayRead(x,a);
1872:   return(0);
1873: }

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

1880:    Logically Collective on Vec

1882:    Input Parameter:
1883: +  x - the vectors
1884: -  n - the number of vectors

1886:    Output Parameter:
1887: .  a - location to put pointer to the array

1889:    Fortran Note:
1890:    This routine is not supported in Fortran.

1892:    Level: intermediate

1894: .seealso: VecGetArray(), VecRestoreArrays()
1895: @*/
1896: PetscErrorCode  VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1897: {
1899:   PetscInt       i;
1900:   PetscScalar    **q;

1906:   if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1907:   PetscMalloc1(n,&q);
1908:   for (i=0; i<n; ++i) {
1909:     VecGetArray(x[i],&q[i]);
1910:   }
1911:   *a = q;
1912:   return(0);
1913: }

1915: /*@C
1916:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1917:    has been called.

1919:    Logically Collective on Vec

1921:    Input Parameters:
1922: +  x - the vector
1923: .  n - the number of vectors
1924: -  a - location of pointer to arrays obtained from VecGetArrays()

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

1932:    Fortran Note:
1933:    This routine is not supported in Fortran.

1935:    Level: intermediate

1937: .seealso: VecGetArrays(), VecRestoreArray()
1938: @*/
1939: PetscErrorCode  VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1940: {
1942:   PetscInt       i;
1943:   PetscScalar    **q = *a;


1950:   for (i=0; i<n; ++i) {
1951:     VecRestoreArray(x[i],&q[i]);
1952:   }
1953:   PetscFree(q);
1954:   return(0);
1955: }

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

1960:    Logically Collective on Vec

1962:    Input Parameters:
1963: +  x - the vector
1964: -  a - location of pointer to array obtained from VecGetArray()

1966:    Level: beginner

1968: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1969:           VecGetArrayPair(), VecRestoreArrayPair()
1970: @*/
1971: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1972: {

1977:   if (x->petscnative) {
1978:    #if defined(PETSC_HAVE_KOKKOS_KERNELS)
1979:     if (x->offloadmask == PETSC_OFFLOAD_VECKOKKOS) {VecKokkosModifyHost(x);}
1980:     else
1981:    #endif
1982:     {
1983:      #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1984:       x->offloadmask = PETSC_OFFLOAD_CPU;
1985:      #endif
1986:     }
1987:   } else {
1988:     (*x->ops->restorearray)(x,a);
1989:   }
1990:   if (a) *a = NULL;
1991:   PetscObjectStateIncrease((PetscObject)x);
1992:   return(0);
1993: }

1995: /*@C
1996:    VecRestoreArrayInPlace - Restores a vector after VecGetArrayInPlace() has been called.

1998:    Logically Collective on Vec

2000:    Input Parameters:
2001: +  x - the vector
2002: -  a - location of pointer to array obtained from VecGetArrayInPlace()

2004:    Level: beginner

2006: .seealso: VecGetArrayInPlace(), VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(),
2007:           VecPlaceArray(), VecRestoreArray2d(), VecGetArrayPair(), VecRestoreArrayPair()
2008: @*/
2009: PetscErrorCode VecRestoreArrayInPlace(Vec x,PetscScalar **a)
2010: {

2015: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
2016:   if (x->offloadmask == PETSC_OFFLOAD_VECKOKKOS) {
2017:     VecKokkosRestoreArrayInPlace(x,a);
2018:     return(0);
2019:   }
2020: #endif

2022: #if defined(PETSC_HAVE_CUDA)
2023:   if (x->petscnative && x->offloadmask == PETSC_OFFLOAD_GPU) {
2024:     PetscBool is_cudatype = PETSC_FALSE;
2025:     PetscObjectTypeCompareAny((PetscObject)x,&is_cudatype,VECSEQCUDA,VECMPICUDA,VECCUDA,"");
2026:     if (is_cudatype) {
2027:       VecCUDARestoreArray(x,a);
2028:       return(0);
2029:     }
2030:   }
2031: #endif
2032:   VecRestoreArray(x,a);
2033:   return(0);
2034: }


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

2040:    Logically Collective on Vec

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

2046:    Level: beginner

2048: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
2049:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite()
2050: @*/
2051: PetscErrorCode VecRestoreArrayWrite(Vec x,PetscScalar **a)
2052: {


2058: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
2059:   if (x->offloadmask == PETSC_OFFLOAD_VECKOKKOS) {
2060:     VecKokkosModifyHost(x);
2061:   } else
2062: #endif
2063:   if (x->petscnative) {
2064: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
2065:     x->offloadmask = PETSC_OFFLOAD_CPU;
2066: #endif
2067:   } else {
2068:     if (x->ops->restorearraywrite) {
2069:       (*x->ops->restorearraywrite)(x,a);
2070:     } else {
2071:       (*x->ops->restorearray)(x,a);
2072:     }
2073:   }

2075:   if (a) *a = NULL;
2076:   PetscObjectStateIncrease((PetscObject)x);
2077:   return(0);
2078: }

2080: /*@C
2081:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

2083:    Not Collective

2085:    Input Parameters:
2086: +  vec - the vector
2087: -  array - the array

2089:    Level: beginner

2091: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
2092: @*/
2093: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
2094: {

2099:   if (x->petscnative) {
2100:     /* nothing */
2101:   } else if (x->ops->restorearrayread) {
2102:     (*x->ops->restorearrayread)(x,a);
2103:   } else {
2104:     (*x->ops->restorearray)(x,(PetscScalar**)a);
2105:   }
2106:   if (a) *a = NULL;
2107:   return(0);
2108: }

2110: /*@C
2111:    VecRestoreArrayReadInPlace - Restore array obtained with VecGetArrayReadInPlace()

2113:    Not Collective

2115:    Input Parameters:
2116: +  vec - the vector
2117: -  array - the array

2119:    Level: beginner

2121: .seealso: VecGetArrayReadInPlace(), VecRestoreArrayInPlace(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
2122: @*/
2123: PetscErrorCode VecRestoreArrayReadInPlace(Vec x,const PetscScalar **a)
2124: {

2128:   VecRestoreArrayRead(x,a);
2129:   return(0);
2130: }

2132: /*@
2133:    VecPlaceArray - Allows one to replace the array in a vector with an
2134:    array provided by the user. This is useful to avoid copying an array
2135:    into a vector.

2137:    Not Collective

2139:    Input Parameters:
2140: +  vec - the vector
2141: -  array - the array

2143:    Notes:
2144:    You can return to the original array with a call to VecResetArray()

2146:    Level: developer

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

2150: @*/
2151: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
2152: {

2159:   if (vec->ops->placearray) {
2160:     (*vec->ops->placearray)(vec,array);
2161:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
2162:   PetscObjectStateIncrease((PetscObject)vec);
2163:   return(0);
2164: }

2166: /*@C
2167:    VecReplaceArray - Allows one to replace the array in a vector with an
2168:    array provided by the user. This is useful to avoid copying an array
2169:    into a vector.

2171:    Not Collective

2173:    Input Parameters:
2174: +  vec - the vector
2175: -  array - the array

2177:    Notes:
2178:    This permanently replaces the array and frees the memory associated
2179:    with the old array.

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

2184:    Not supported from Fortran

2186:    Level: developer

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

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

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


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

2209:    This function has semantics similar to VecGetArray():  the pointer
2210:    returned by this function points to a consistent view of the vector
2211:    data.  This may involve a copy operation of data from the host to the
2212:    device if the data on the device is out of date.  If the device
2213:    memory hasn't been allocated previously it will be allocated as part
2214:    of this function call.  VecCUDAGetArray() assumes that
2215:    the user will modify the vector data.  This is similar to
2216:    intent(inout) in fortran.

2218:    The CUDA device pointer has to be released by calling
2219:    VecCUDARestoreArray().  Upon restoring the vector data
2220:    the data on the host will be marked as out of date.  A subsequent
2221:    access of the host data will thus incur a data transfer from the
2222:    device to the host.


2225:    Input Parameter:
2226: .  v - the vector

2228:    Output Parameter:
2229: .  a - the CUDA device pointer

2231:    Fortran note:
2232:    This function is not currently available from Fortran.

2234:    Level: intermediate

2236: .seealso: VecCUDARestoreArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2237: @*/
2238: PETSC_EXTERN PetscErrorCode VecCUDAGetArray(Vec v, PetscScalar **a)
2239: {
2240: #if defined(PETSC_HAVE_CUDA)
2242: #endif

2246: #if defined(PETSC_HAVE_CUDA)
2247:   *a   = 0;
2248:   VecCUDACopyToGPU(v);
2249:   *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2250: #endif
2251:   return(0);
2252: }

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

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

2261:    Input Parameter:
2262: +  v - the vector
2263: -  a - the CUDA device pointer.  This pointer is invalid after
2264:        VecCUDARestoreArray() returns.

2266:    Fortran note:
2267:    This function is not currently available from Fortran.

2269:    Level: intermediate

2271: .seealso: VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2272: @*/
2273: PETSC_EXTERN PetscErrorCode VecCUDARestoreArray(Vec v, PetscScalar **a)
2274: {

2279: #if defined(PETSC_HAVE_CUDA)
2280:   v->offloadmask = PETSC_OFFLOAD_GPU;
2281: #endif

2283:   PetscObjectStateIncrease((PetscObject)v);
2284:   return(0);
2285: }

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

2290:    This function is analogous to VecGetArrayRead():  The pointer
2291:    returned by this function points to a consistent view of the vector
2292:    data.  This may involve a copy operation of data from the host to the
2293:    device if the data on the device is out of date.  If the device
2294:    memory hasn't been allocated previously it will be allocated as part
2295:    of this function call.  VecCUDAGetArrayRead() assumes that the
2296:    user will not modify the vector data.  This is analgogous to
2297:    intent(in) in Fortran.

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

2305:    Input Parameter:
2306: .  v - the vector

2308:    Output Parameter:
2309: .  a - the CUDA pointer.

2311:    Fortran note:
2312:    This function is not currently available from Fortran.

2314:    Level: intermediate

2316: .seealso: VecCUDARestoreArrayRead(), VecCUDAGetArray(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2317: @*/
2318: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayRead(Vec v, const PetscScalar **a)
2319: {
2320: #if defined(PETSC_HAVE_CUDA)
2322: #endif

2326: #if defined(PETSC_HAVE_CUDA)
2327:   *a   = 0;
2328:   VecCUDACopyToGPU(v);
2329:   *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2330: #endif
2331:   return(0);
2332: }

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

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

2342:    Input Parameter:
2343: +  v - the vector
2344: -  a - the CUDA device pointer.  This pointer is invalid after
2345:        VecCUDARestoreArrayRead() returns.

2347:    Fortran note:
2348:    This function is not currently available from Fortran.

2350:    Level: intermediate

2352: .seealso: VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2353: @*/
2354: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayRead(Vec v, const PetscScalar **a)
2355: {
2358:   *a = NULL;
2359:   return(0);
2360: }

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

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

2371:    The device pointer needs to be released with
2372:    VecCUDARestoreArrayWrite().  When the pointer is released the
2373:    host data of the vector is marked as out of data.  Subsequent access
2374:    of the host data with e.g. VecGetArray() incurs a device to host data
2375:    transfer.


2378:    Input Parameter:
2379: .  v - the vector

2381:    Output Parameter:
2382: .  a - the CUDA pointer

2384:    Fortran note:
2385:    This function is not currently available from Fortran.

2387:    Level: advanced

2389: .seealso: VecCUDARestoreArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2390: @*/
2391: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayWrite(Vec v, PetscScalar **a)
2392: {
2393: #if defined(PETSC_HAVE_CUDA)
2395: #endif

2399: #if defined(PETSC_HAVE_CUDA)
2400:   *a   = 0;
2401:   VecCUDAAllocateCheck(v);
2402:   *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2403: #endif
2404:   return(0);
2405: }

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

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

2414:    Input Parameter:
2415: +  v - the vector
2416: -  a - the CUDA device pointer.  This pointer is invalid after
2417:        VecCUDARestoreArrayWrite() returns.

2419:    Fortran note:
2420:    This function is not currently available from Fortran.

2422:    Level: intermediate

2424: .seealso: VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2425: @*/
2426: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayWrite(Vec v, PetscScalar **a)
2427: {

2432: #if defined(PETSC_HAVE_CUDA)
2433:   v->offloadmask = PETSC_OFFLOAD_GPU;
2434: #endif

2436:   PetscObjectStateIncrease((PetscObject)v);
2437:   return(0);
2438: }

2440: /*@C
2441:    VecCUDAPlaceArray - Allows one to replace the GPU array in a vector with a
2442:    GPU array provided by the user. This is useful to avoid copying an
2443:    array into a vector.

2445:    Not Collective

2447:    Input Parameters:
2448: +  vec - the vector
2449: -  array - the GPU array

2451:    Notes:
2452:    You can return to the original GPU array with a call to VecCUDAResetArray()
2453:    It is not possible to use VecCUDAPlaceArray() and VecPlaceArray() at the
2454:    same time on the same vector.

2456:    Level: developer

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

2460: @*/
2461: PetscErrorCode VecCUDAPlaceArray(Vec vin,const PetscScalar a[])
2462: {

2467: #if defined(PETSC_HAVE_CUDA)
2468:   VecCUDACopyToGPU(vin);
2469:   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()");
2470:   ((Vec_Seq*)vin->data)->unplacedarray  = (PetscScalar *) ((Vec_CUDA*)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2471:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2472:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2473: #endif
2474:   PetscObjectStateIncrease((PetscObject)vin);
2475:   return(0);
2476: }

2478: /*@C
2479:    VecCUDAReplaceArray - Allows one to replace the GPU array in a vector
2480:    with a GPU array provided by the user. This is useful to avoid copying
2481:    a GPU array into a vector.

2483:    Not Collective

2485:    Input Parameters:
2486: +  vec - the vector
2487: -  array - the GPU array

2489:    Notes:
2490:    This permanently replaces the GPU array and frees the memory associated
2491:    with the old GPU array.

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

2496:    Not supported from Fortran

2498:    Level: developer

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

2502: @*/
2503: PetscErrorCode VecCUDAReplaceArray(Vec vin,const PetscScalar a[])
2504: {
2505: #if defined(PETSC_HAVE_CUDA)
2506:   cudaError_t err;
2507: #endif

2512: #if defined(PETSC_HAVE_CUDA)
2513:   err = cudaFree(((Vec_CUDA*)vin->spptr)->GPUarray);CHKERRCUDA(err);
2514:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2515:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2516: #endif
2517:   PetscObjectStateIncrease((PetscObject)vin);
2518:   return(0);
2519: }

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

2525:    Not Collective

2527:    Input Parameters:
2528: .  vec - the vector

2530:    Level: developer

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

2534: @*/
2535: PetscErrorCode VecCUDAResetArray(Vec vin)
2536: {

2541: #if defined(PETSC_HAVE_CUDA)
2542:   VecCUDACopyToGPU(vin);
2543:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar *) ((Vec_Seq*)vin->data)->unplacedarray;
2544:   ((Vec_Seq*)vin->data)->unplacedarray = 0;
2545:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2546: #endif
2547:   PetscObjectStateIncrease((PetscObject)vin);
2548:   return(0);
2549: }

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

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

2558:     Collective on Vec

2560:     Input Parameters:
2561: +   x - a vector to mimic
2562: -   n - the number of vectors to obtain

2564:     Output Parameters:
2565: +   y - Fortran90 pointer to the array of vectors
2566: -   ierr - error code

2568:     Example of Usage:
2569: .vb
2570: #include <petsc/finclude/petscvec.h>
2571:     use petscvec

2573:     Vec x
2574:     Vec, pointer :: y(:)
2575:     ....
2576:     call VecDuplicateVecsF90(x,2,y,ierr)
2577:     call VecSet(y(2),alpha,ierr)
2578:     call VecSet(y(2),alpha,ierr)
2579:     ....
2580:     call VecDestroyVecsF90(2,y,ierr)
2581: .ve

2583:     Notes:
2584:     Not yet supported for all F90 compilers

2586:     Use VecDestroyVecsF90() to free the space.

2588:     Level: beginner

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

2592: M*/

2594: /*MC
2595:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2596:     VecGetArrayF90().

2598:     Synopsis:
2599:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2601:     Logically Collective on Vec

2603:     Input Parameters:
2604: +   x - vector
2605: -   xx_v - the Fortran90 pointer to the array

2607:     Output Parameter:
2608: .   ierr - error code

2610:     Example of Usage:
2611: .vb
2612: #include <petsc/finclude/petscvec.h>
2613:     use petscvec

2615:     PetscScalar, pointer :: xx_v(:)
2616:     ....
2617:     call VecGetArrayF90(x,xx_v,ierr)
2618:     xx_v(3) = a
2619:     call VecRestoreArrayF90(x,xx_v,ierr)
2620: .ve

2622:     Level: beginner

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

2626: M*/

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

2631:     Synopsis:
2632:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

2634:     Collective on Vec

2636:     Input Parameters:
2637: +   n - the number of vectors previously obtained
2638: -   x - pointer to array of vector pointers

2640:     Output Parameter:
2641: .   ierr - error code

2643:     Notes:
2644:     Not yet supported for all F90 compilers

2646:     Level: beginner

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

2650: M*/

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

2658:     Synopsis:
2659:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2661:     Logically Collective on Vec

2663:     Input Parameter:
2664: .   x - vector

2666:     Output Parameters:
2667: +   xx_v - the Fortran90 pointer to the array
2668: -   ierr - error code

2670:     Example of Usage:
2671: .vb
2672: #include <petsc/finclude/petscvec.h>
2673:     use petscvec

2675:     PetscScalar, pointer :: xx_v(:)
2676:     ....
2677:     call VecGetArrayF90(x,xx_v,ierr)
2678:     xx_v(3) = a
2679:     call VecRestoreArrayF90(x,xx_v,ierr)
2680: .ve

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

2684:     Level: beginner

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

2688: M*/

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

2696:     Synopsis:
2697:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2699:     Logically Collective on Vec

2701:     Input Parameter:
2702: .   x - vector

2704:     Output Parameters:
2705: +   xx_v - the Fortran90 pointer to the array
2706: -   ierr - error code

2708:     Example of Usage:
2709: .vb
2710: #include <petsc/finclude/petscvec.h>
2711:     use petscvec

2713:     PetscScalar, pointer :: xx_v(:)
2714:     ....
2715:     call VecGetArrayReadF90(x,xx_v,ierr)
2716:     a = xx_v(3)
2717:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2718: .ve

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

2722:     Level: beginner

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

2726: M*/

2728: /*MC
2729:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2730:     VecGetArrayReadF90().

2732:     Synopsis:
2733:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2735:     Logically Collective on Vec

2737:     Input Parameters:
2738: +   x - vector
2739: -   xx_v - the Fortran90 pointer to the array

2741:     Output Parameter:
2742: .   ierr - error code

2744:     Example of Usage:
2745: .vb
2746: #include <petsc/finclude/petscvec.h>
2747:     use petscvec

2749:     PetscScalar, pointer :: xx_v(:)
2750:     ....
2751:     call VecGetArrayReadF90(x,xx_v,ierr)
2752:     a = xx_v(3)
2753:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2754: .ve

2756:     Level: beginner

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

2760: M*/

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

2767:    Logically Collective

2769:    Input Parameter:
2770: +  x - the vector
2771: .  m - first dimension of two dimensional array
2772: .  n - second dimension of two dimensional array
2773: .  mstart - first index you will use in first coordinate direction (often 0)
2774: -  nstart - first index in the second coordinate direction (often 0)

2776:    Output Parameter:
2777: .  a - location to put pointer to the array

2779:    Level: developer

2781:   Notes:
2782:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2783:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2784:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2785:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

2789: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2790:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2791:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2792: @*/
2793: PetscErrorCode  VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2794: {
2796:   PetscInt       i,N;
2797:   PetscScalar    *aa;

2803:   VecGetLocalSize(x,&N);
2804:   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);
2805:   VecGetArray(x,&aa);

2807:   PetscMalloc1(m,a);
2808:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2809:   *a -= mstart;
2810:   return(0);
2811: }

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

2818:    Logically Collective

2820:    Input Parameter:
2821: +  x - the vector
2822: .  m - first dimension of two dimensional array
2823: .  n - second dimension of two dimensional array
2824: .  mstart - first index you will use in first coordinate direction (often 0)
2825: -  nstart - first index in the second coordinate direction (often 0)

2827:    Output Parameter:
2828: .  a - location to put pointer to the array

2830:    Level: developer

2832:   Notes:
2833:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2834:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2835:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2836:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

2842: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2843:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2844:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2845: @*/
2846: PetscErrorCode  VecGetArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2847: {
2849:   PetscInt       i,N;
2850:   PetscScalar    *aa;

2856:   VecGetLocalSize(x,&N);
2857:   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);
2858:   VecGetArrayWrite(x,&aa);

2860:   PetscMalloc1(m,a);
2861:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2862:   *a -= mstart;
2863:   return(0);
2864: }

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

2869:    Logically Collective

2871:    Input Parameters:
2872: +  x - the vector
2873: .  m - first dimension of two dimensional array
2874: .  n - second dimension of the two dimensional array
2875: .  mstart - first index you will use in first coordinate direction (often 0)
2876: .  nstart - first index in the second coordinate direction (often 0)
2877: -  a - location of pointer to array obtained from VecGetArray2d()

2879:    Level: developer

2881:    Notes:
2882:    For regular PETSc vectors this routine does not involve any copies. For
2883:    any special vectors that do not store local vector data in a contiguous
2884:    array, this routine will copy the data back into the underlying
2885:    vector data structure from the array obtained with VecGetArray().

2887:    This routine actually zeros out the a pointer.

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

2902:   dummy = (void*)(*a + mstart);
2903:   PetscFree(dummy);
2904:   VecRestoreArray(x,NULL);
2905:   return(0);
2906: }

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

2911:    Logically Collective

2913:    Input Parameters:
2914: +  x - the vector
2915: .  m - first dimension of two dimensional array
2916: .  n - second dimension of the two dimensional array
2917: .  mstart - first index you will use in first coordinate direction (often 0)
2918: .  nstart - first index in the second coordinate direction (often 0)
2919: -  a - location of pointer to array obtained from VecGetArray2d()

2921:    Level: developer

2923:    Notes:
2924:    For regular PETSc vectors this routine does not involve any copies. For
2925:    any special vectors that do not store local vector data in a contiguous
2926:    array, this routine will copy the data back into the underlying
2927:    vector data structure from the array obtained with VecGetArray().

2929:    This routine actually zeros out the a pointer.

2931: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2932:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2933:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2934: @*/
2935: PetscErrorCode  VecRestoreArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2936: {
2938:   void           *dummy;

2944:   dummy = (void*)(*a + mstart);
2945:   PetscFree(dummy);
2946:   VecRestoreArrayWrite(x,NULL);
2947:   return(0);
2948: }

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

2955:    Logically Collective

2957:    Input Parameter:
2958: +  x - the vector
2959: .  m - first dimension of two dimensional array
2960: -  mstart - first index you will use in first coordinate direction (often 0)

2962:    Output Parameter:
2963: .  a - location to put pointer to the array

2965:    Level: developer

2967:   Notes:
2968:    For a vector obtained from DMCreateLocalVector() mstart are likely
2969:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2970:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

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

2987:   VecGetLocalSize(x,&N);
2988:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2989:   VecGetArray(x,a);
2990:   *a  -= mstart;
2991:   return(0);
2992: }

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

2999:    Logically Collective

3001:    Input Parameter:
3002: +  x - the vector
3003: .  m - first dimension of two dimensional array
3004: -  mstart - first index you will use in first coordinate direction (often 0)

3006:    Output Parameter:
3007: .  a - location to put pointer to the array

3009:    Level: developer

3011:   Notes:
3012:    For a vector obtained from DMCreateLocalVector() mstart are likely
3013:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3014:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3018: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3019:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3020:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3021: @*/
3022: PetscErrorCode  VecGetArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3023: {
3025:   PetscInt       N;

3031:   VecGetLocalSize(x,&N);
3032:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3033:   VecGetArrayWrite(x,a);
3034:   *a  -= mstart;
3035:   return(0);
3036: }

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

3041:    Logically Collective

3043:    Input Parameters:
3044: +  x - the vector
3045: .  m - first dimension of two dimensional array
3046: .  mstart - first index you will use in first coordinate direction (often 0)
3047: -  a - location of pointer to array obtained from VecGetArray21()

3049:    Level: developer

3051:    Notes:
3052:    For regular PETSc vectors this routine does not involve any copies. For
3053:    any special vectors that do not store local vector data in a contiguous
3054:    array, this routine will copy the data back into the underlying
3055:    vector data structure from the array obtained with VecGetArray1d().

3057:    This routine actually zeros out the a pointer.

3059: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3060:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3061:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3062: @*/
3063: PetscErrorCode  VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3064: {

3070:   VecRestoreArray(x,NULL);
3071:   return(0);
3072: }

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

3077:    Logically Collective

3079:    Input Parameters:
3080: +  x - the vector
3081: .  m - first dimension of two dimensional array
3082: .  mstart - first index you will use in first coordinate direction (often 0)
3083: -  a - location of pointer to array obtained from VecGetArray21()

3085:    Level: developer

3087:    Notes:
3088:    For regular PETSc vectors this routine does not involve any copies. For
3089:    any special vectors that do not store local vector data in a contiguous
3090:    array, this routine will copy the data back into the underlying
3091:    vector data structure from the array obtained with VecGetArray1d().

3093:    This routine actually zeros out the a pointer.

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

3097: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3098:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3099:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3100: @*/
3101: PetscErrorCode  VecRestoreArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3102: {

3108:   VecRestoreArrayWrite(x,NULL);
3109:   return(0);
3110: }

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

3117:    Logically Collective

3119:    Input Parameter:
3120: +  x - the vector
3121: .  m - first dimension of three dimensional array
3122: .  n - second dimension of three dimensional array
3123: .  p - third dimension of three dimensional array
3124: .  mstart - first index you will use in first coordinate direction (often 0)
3125: .  nstart - first index in the second coordinate direction (often 0)
3126: -  pstart - first index in the third coordinate direction (often 0)

3128:    Output Parameter:
3129: .  a - location to put pointer to the array

3131:    Level: developer

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

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

3141: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3142:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3143:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3144: @*/
3145: PetscErrorCode  VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3146: {
3148:   PetscInt       i,N,j;
3149:   PetscScalar    *aa,**b;

3155:   VecGetLocalSize(x,&N);
3156:   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);
3157:   VecGetArray(x,&aa);

3159:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3160:   b    = (PetscScalar**)((*a) + m);
3161:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3162:   for (i=0; i<m; i++)
3163:     for (j=0; j<n; j++)
3164:       b[i*n+j] = aa + i*n*p + j*p - pstart;

3166:   *a -= mstart;
3167:   return(0);
3168: }

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

3175:    Logically Collective

3177:    Input Parameter:
3178: +  x - the vector
3179: .  m - first dimension of three dimensional array
3180: .  n - second dimension of three dimensional array
3181: .  p - third dimension of three dimensional array
3182: .  mstart - first index you will use in first coordinate direction (often 0)
3183: .  nstart - first index in the second coordinate direction (often 0)
3184: -  pstart - first index in the third coordinate direction (often 0)

3186:    Output Parameter:
3187: .  a - location to put pointer to the array

3189:    Level: developer

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

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

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

3201: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3202:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3203:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3204: @*/
3205: PetscErrorCode  VecGetArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3206: {
3208:   PetscInt       i,N,j;
3209:   PetscScalar    *aa,**b;

3215:   VecGetLocalSize(x,&N);
3216:   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);
3217:   VecGetArrayWrite(x,&aa);

3219:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3220:   b    = (PetscScalar**)((*a) + m);
3221:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3222:   for (i=0; i<m; i++)
3223:     for (j=0; j<n; j++)
3224:       b[i*n+j] = aa + i*n*p + j*p - pstart;

3226:   *a -= mstart;
3227:   return(0);
3228: }

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

3233:    Logically Collective

3235:    Input Parameters:
3236: +  x - the vector
3237: .  m - first dimension of three dimensional array
3238: .  n - second dimension of the three dimensional array
3239: .  p - third dimension of the three dimensional array
3240: .  mstart - first index you will use in first coordinate direction (often 0)
3241: .  nstart - first index in the second coordinate direction (often 0)
3242: .  pstart - first index in the third coordinate direction (often 0)
3243: -  a - location of pointer to array obtained from VecGetArray3d()

3245:    Level: developer

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

3253:    This routine actually zeros out the a pointer.

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

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

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

3277:    Logically Collective

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

3289:    Level: developer

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

3297:    This routine actually zeros out the a pointer.

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

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

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

3323:    Logically Collective

3325:    Input Parameter:
3326: +  x - the vector
3327: .  m - first dimension of four dimensional array
3328: .  n - second dimension of four dimensional array
3329: .  p - third dimension of four dimensional array
3330: .  q - fourth dimension of four dimensional array
3331: .  mstart - first index you will use in first coordinate direction (often 0)
3332: .  nstart - first index in the second coordinate direction (often 0)
3333: .  pstart - first index in the third coordinate direction (often 0)
3334: -  qstart - first index in the fourth coordinate direction (often 0)

3336:    Output Parameter:
3337: .  a - location to put pointer to the array

3339:    Level: beginner

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

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

3349: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3350:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3351:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3352: @*/
3353: PetscErrorCode  VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3354: {
3356:   PetscInt       i,N,j,k;
3357:   PetscScalar    *aa,***b,**c;

3363:   VecGetLocalSize(x,&N);
3364:   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);
3365:   VecGetArray(x,&aa);

3367:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3368:   b    = (PetscScalar***)((*a) + m);
3369:   c    = (PetscScalar**)(b + m*n);
3370:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3371:   for (i=0; i<m; i++)
3372:     for (j=0; j<n; j++)
3373:       b[i*n+j] = c + i*n*p + j*p - pstart;
3374:   for (i=0; i<m; i++)
3375:     for (j=0; j<n; j++)
3376:       for (k=0; k<p; k++)
3377:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3378:   *a -= mstart;
3379:   return(0);
3380: }

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

3387:    Logically Collective

3389:    Input Parameter:
3390: +  x - the vector
3391: .  m - first dimension of four dimensional array
3392: .  n - second dimension of four dimensional array
3393: .  p - third dimension of four dimensional array
3394: .  q - fourth dimension of four dimensional array
3395: .  mstart - first index you will use in first coordinate direction (often 0)
3396: .  nstart - first index in the second coordinate direction (often 0)
3397: .  pstart - first index in the third coordinate direction (often 0)
3398: -  qstart - first index in the fourth coordinate direction (often 0)

3400:    Output Parameter:
3401: .  a - location to put pointer to the array

3403:    Level: beginner

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

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

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

3415: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3416:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3417:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3418: @*/
3419: PetscErrorCode  VecGetArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3420: {
3422:   PetscInt       i,N,j,k;
3423:   PetscScalar    *aa,***b,**c;

3429:   VecGetLocalSize(x,&N);
3430:   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);
3431:   VecGetArrayWrite(x,&aa);

3433:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3434:   b    = (PetscScalar***)((*a) + m);
3435:   c    = (PetscScalar**)(b + m*n);
3436:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3437:   for (i=0; i<m; i++)
3438:     for (j=0; j<n; j++)
3439:       b[i*n+j] = c + i*n*p + j*p - pstart;
3440:   for (i=0; i<m; i++)
3441:     for (j=0; j<n; j++)
3442:       for (k=0; k<p; k++)
3443:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3444:   *a -= mstart;
3445:   return(0);
3446: }

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

3451:    Logically Collective

3453:    Input Parameters:
3454: +  x - the vector
3455: .  m - first dimension of four dimensional array
3456: .  n - second dimension of the four dimensional array
3457: .  p - third dimension of the four dimensional array
3458: .  q - fourth dimension of the four dimensional array
3459: .  mstart - first index you will use in first coordinate direction (often 0)
3460: .  nstart - first index in the second coordinate direction (often 0)
3461: .  pstart - first index in the third coordinate direction (often 0)
3462: .  qstart - first index in the fourth coordinate direction (often 0)
3463: -  a - location of pointer to array obtained from VecGetArray4d()

3465:    Level: beginner

3467:    Notes:
3468:    For regular PETSc vectors this routine does not involve any copies. For
3469:    any special vectors that do not store local vector data in a contiguous
3470:    array, this routine will copy the data back into the underlying
3471:    vector data structure from the array obtained with VecGetArray().

3473:    This routine actually zeros out the a pointer.

3475: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3476:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3477:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3478: @*/
3479: PetscErrorCode  VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3480: {
3482:   void           *dummy;

3488:   dummy = (void*)(*a + mstart);
3489:   PetscFree(dummy);
3490:   VecRestoreArray(x,NULL);
3491:   return(0);
3492: }

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

3497:    Logically Collective

3499:    Input Parameters:
3500: +  x - the vector
3501: .  m - first dimension of four dimensional array
3502: .  n - second dimension of the four dimensional array
3503: .  p - third dimension of the four dimensional array
3504: .  q - fourth dimension of the four dimensional array
3505: .  mstart - first index you will use in first coordinate direction (often 0)
3506: .  nstart - first index in the second coordinate direction (often 0)
3507: .  pstart - first index in the third coordinate direction (often 0)
3508: .  qstart - first index in the fourth coordinate direction (often 0)
3509: -  a - location of pointer to array obtained from VecGetArray4d()

3511:    Level: beginner

3513:    Notes:
3514:    For regular PETSc vectors this routine does not involve any copies. For
3515:    any special vectors that do not store local vector data in a contiguous
3516:    array, this routine will copy the data back into the underlying
3517:    vector data structure from the array obtained with VecGetArray().

3519:    This routine actually zeros out the a pointer.

3521: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3522:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3523:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3524: @*/
3525: PetscErrorCode  VecRestoreArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3526: {
3528:   void           *dummy;

3534:   dummy = (void*)(*a + mstart);
3535:   PetscFree(dummy);
3536:   VecRestoreArrayWrite(x,NULL);
3537:   return(0);
3538: }

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

3545:    Logically Collective

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

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

3557:    Level: developer

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

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

3567: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3568:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3569:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3570: @*/
3571: PetscErrorCode  VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3572: {
3573:   PetscErrorCode    ierr;
3574:   PetscInt          i,N;
3575:   const PetscScalar *aa;

3581:   VecGetLocalSize(x,&N);
3582:   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);
3583:   VecGetArrayRead(x,&aa);

3585:   PetscMalloc1(m,a);
3586:   for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
3587:   *a -= mstart;
3588:   return(0);
3589: }

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

3594:    Logically Collective

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

3604:    Level: developer

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

3612:    This routine actually zeros out the a pointer.

3614: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3615:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3616:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3617: @*/
3618: PetscErrorCode  VecRestoreArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3619: {
3621:   void           *dummy;

3627:   dummy = (void*)(*a + mstart);
3628:   PetscFree(dummy);
3629:   VecRestoreArrayRead(x,NULL);
3630:   return(0);
3631: }

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

3638:    Logically Collective

3640:    Input Parameter:
3641: +  x - the vector
3642: .  m - first dimension of two dimensional array
3643: -  mstart - first index you will use in first 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 are likely
3652:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3653:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3657: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3658:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3659:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3660: @*/
3661: PetscErrorCode  VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3662: {
3664:   PetscInt       N;

3670:   VecGetLocalSize(x,&N);
3671:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3672:   VecGetArrayRead(x,(const PetscScalar**)a);
3673:   *a  -= mstart;
3674:   return(0);
3675: }

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

3680:    Logically Collective

3682:    Input Parameters:
3683: +  x - the vector
3684: .  m - first dimension of two dimensional array
3685: .  mstart - first index you will use in first coordinate direction (often 0)
3686: -  a - location of pointer to array obtained from VecGetArray21()

3688:    Level: developer

3690:    Notes:
3691:    For regular PETSc vectors this routine does not involve any copies. For
3692:    any special vectors that do not store local vector data in a contiguous
3693:    array, this routine will copy the data back into the underlying
3694:    vector data structure from the array obtained with VecGetArray1dRead().

3696:    This routine actually zeros out the a pointer.

3698: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3699:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3700:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3701: @*/
3702: PetscErrorCode  VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3703: {

3709:   VecRestoreArrayRead(x,NULL);
3710:   return(0);
3711: }


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

3719:    Logically Collective

3721:    Input Parameter:
3722: +  x - the vector
3723: .  m - first dimension of three dimensional array
3724: .  n - second dimension of three dimensional array
3725: .  p - third dimension of three dimensional array
3726: .  mstart - first index you will use in first coordinate direction (often 0)
3727: .  nstart - first index in the second coordinate direction (often 0)
3728: -  pstart - first index in the third coordinate direction (often 0)

3730:    Output Parameter:
3731: .  a - location to put pointer to the array

3733:    Level: developer

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

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

3743: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3744:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3745:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3746: @*/
3747: PetscErrorCode  VecGetArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3748: {
3749:   PetscErrorCode    ierr;
3750:   PetscInt          i,N,j;
3751:   const PetscScalar *aa;
3752:   PetscScalar       **b;

3758:   VecGetLocalSize(x,&N);
3759:   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);
3760:   VecGetArrayRead(x,&aa);

3762:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3763:   b    = (PetscScalar**)((*a) + m);
3764:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3765:   for (i=0; i<m; i++)
3766:     for (j=0; j<n; j++)
3767:       b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;

3769:   *a -= mstart;
3770:   return(0);
3771: }

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

3776:    Logically Collective

3778:    Input Parameters:
3779: +  x - the vector
3780: .  m - first dimension of three dimensional array
3781: .  n - second dimension of the three dimensional array
3782: .  p - third dimension of the three dimensional array
3783: .  mstart - first index you will use in first coordinate direction (often 0)
3784: .  nstart - first index in the second coordinate direction (often 0)
3785: .  pstart - first index in the third coordinate direction (often 0)
3786: -  a - location of pointer to array obtained from VecGetArray3dRead()

3788:    Level: developer

3790:    Notes:
3791:    For regular PETSc vectors this routine does not involve any copies. For
3792:    any special vectors that do not store local vector data in a contiguous
3793:    array, this routine will copy the data back into the underlying
3794:    vector data structure from the array obtained with VecGetArray().

3796:    This routine actually zeros out the a pointer.

3798: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3799:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3800:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3801: @*/
3802: PetscErrorCode  VecRestoreArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3803: {
3805:   void           *dummy;

3811:   dummy = (void*)(*a + mstart);
3812:   PetscFree(dummy);
3813:   VecRestoreArrayRead(x,NULL);
3814:   return(0);
3815: }

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

3822:    Logically Collective

3824:    Input Parameter:
3825: +  x - the vector
3826: .  m - first dimension of four dimensional array
3827: .  n - second dimension of four dimensional array
3828: .  p - third dimension of four dimensional array
3829: .  q - fourth dimension of four dimensional array
3830: .  mstart - first index you will use in first coordinate direction (often 0)
3831: .  nstart - first index in the second coordinate direction (often 0)
3832: .  pstart - first index in the third coordinate direction (often 0)
3833: -  qstart - first index in the fourth coordinate direction (often 0)

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

3838:    Level: beginner

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

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

3848: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3849:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3850:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3851: @*/
3852: PetscErrorCode  VecGetArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3853: {
3854:   PetscErrorCode    ierr;
3855:   PetscInt          i,N,j,k;
3856:   const PetscScalar *aa;
3857:   PetscScalar       ***b,**c;

3863:   VecGetLocalSize(x,&N);
3864:   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);
3865:   VecGetArrayRead(x,&aa);

3867:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3868:   b    = (PetscScalar***)((*a) + m);
3869:   c    = (PetscScalar**)(b + m*n);
3870:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3871:   for (i=0; i<m; i++)
3872:     for (j=0; j<n; j++)
3873:       b[i*n+j] = c + i*n*p + j*p - pstart;
3874:   for (i=0; i<m; i++)
3875:     for (j=0; j<n; j++)
3876:       for (k=0; k<p; k++)
3877:         c[i*n*p+j*p+k] = (PetscScalar*) aa + i*n*p*q + j*p*q + k*q - qstart;
3878:   *a -= mstart;
3879:   return(0);
3880: }

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

3885:    Logically Collective

3887:    Input Parameters:
3888: +  x - the vector
3889: .  m - first dimension of four dimensional array
3890: .  n - second dimension of the four dimensional array
3891: .  p - third dimension of the four dimensional array
3892: .  q - fourth dimension of the four dimensional array
3893: .  mstart - first index you will use in first coordinate direction (often 0)
3894: .  nstart - first index in the second coordinate direction (often 0)
3895: .  pstart - first index in the third coordinate direction (often 0)
3896: .  qstart - first index in the fourth coordinate direction (often 0)
3897: -  a - location of pointer to array obtained from VecGetArray4dRead()

3899:    Level: beginner

3901:    Notes:
3902:    For regular PETSc vectors this routine does not involve any copies. For
3903:    any special vectors that do not store local vector data in a contiguous
3904:    array, this routine will copy the data back into the underlying
3905:    vector data structure from the array obtained with VecGetArray().

3907:    This routine actually zeros out the a pointer.

3909: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3910:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3911:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3912: @*/
3913: PetscErrorCode  VecRestoreArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3914: {
3916:   void           *dummy;

3922:   dummy = (void*)(*a + mstart);
3923:   PetscFree(dummy);
3924:   VecRestoreArrayRead(x,NULL);
3925:   return(0);
3926: }

3928: #if defined(PETSC_USE_DEBUG)

3930: /*@
3931:    VecLockGet  - Gets the current lock status of a vector

3933:    Logically Collective on Vec

3935:    Input Parameter:
3936: .  x - the vector

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

3942:    Level: beginner

3944: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop()
3945: @*/
3946: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
3947: {
3950:   *state = x->lock;
3951:   return(0);
3952: }

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

3957:    Logically Collective on Vec

3959:    Input Parameter:
3960: .  x - the vector

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

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

3968:    Level: beginner

3970: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPop(), VecLockGet()
3971: @*/
3972: PetscErrorCode VecLockReadPush(Vec x)
3973: {
3976:   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");
3977:   x->lock++;
3978:   return(0);
3979: }

3981: /*@
3982:    VecLockReadPop  - Pops a read-only lock from a vector

3984:    Logically Collective on Vec

3986:    Input Parameter:
3987: .  x - the vector

3989:    Level: beginner

3991: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockGet()
3992: @*/
3993: PetscErrorCode VecLockReadPop(Vec x)
3994: {
3997:   x->lock--;
3998:   if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector has been unlocked from read-only access too many times");
3999:   return(0);
4000: }

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

4005:    Logically Collective on Vec

4007:    Input Parameter:
4008: +  x   - the vector
4009: -  flg - PETSC_TRUE to lock the vector for writing; PETSC_FALSE to unlock it.

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


4018:        VecGetArray(x,&xdata); // begin phase
4019:        VecLockWriteSet_Private(v,PETSC_TRUE);

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

4023:        VecRestoreArray(x,&vdata); // end phase
4024:        VecLockWriteSet_Private(v,PETSC_FALSE);

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

4029:    Level: beginner

4031: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop(), VecLockGet()
4032: @*/
4033: PetscErrorCode VecLockWriteSet_Private(Vec x,PetscBool flg)
4034: {
4037:   if (flg) {
4038:     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");
4039:     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");
4040:     else x->lock = -1;
4041:   } else {
4042:     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");
4043:     x->lock = 0;
4044:   }
4045:   return(0);
4046: }

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

4051:    Level: deprecated

4053: .seealso: VecLockReadPush()
4054: @*/
4055: PetscErrorCode VecLockPush(Vec x)
4056: {
4059:   VecLockReadPush(x);
4060:   return(0);
4061: }

4063: /*@
4064:    VecLockPop  - Pops a read-only lock from a vector

4066:    Level: deprecated

4068: .seealso: VecLockReadPop()
4069: @*/
4070: PetscErrorCode VecLockPop(Vec x)
4071: {
4074:   VecLockReadPop(x);
4075:   return(0);
4076: }

4078: #endif