Actual source code: rvector.c

petsc-3.14.3 2021-01-09
Report Typos and Errors
  1: /*
  2:      Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
  3:    These are the vector functions the user calls.
  4: */
  5: #include "petsc/private/sfimpl.h"
  6: #include <petsc/private/vecimpl.h>
  7: #if defined(PETSC_HAVE_CUDA)
  8: #include <../src/vec/vec/impls/dvecimpl.h>
  9: #include <petsc/private/cudavecimpl.h>
 10: #endif
 11: static PetscInt VecGetSubVectorSavedStateId = -1;

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

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

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

 46:    Logically Collective on Vec

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

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

 54:    Level: advanced

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

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

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

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

 81:    Collective on Vec

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

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

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

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

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

105:    Level: intermediate


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

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

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

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

132:    Collective on Vec

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

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

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

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

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

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

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

155:    Level: intermediate


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

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

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

174:    Collective on Vec

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

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

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

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

195:    Level: intermediate

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


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

206: @*/

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


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

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

237:    Not Collective

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

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

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

254:    Level: intermediate

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

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


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

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


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

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

290:    Collective on Vec

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

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

299:    Level: intermediate


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

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

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

327:    Collective on Vec

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

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

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

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


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

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

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

362:    Collective on Vec

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

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

371:    Level: intermediate

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

376:    This returns the smallest index with the minumum value


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

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

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

399:    Collective on Vec

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

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

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

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

416:    Level: intermediate

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

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

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

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

442:    Not collective on Vec

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

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

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

455:    Level: intermediate


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

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

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

493:    Logically Collective on Vec

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

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

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

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

512:    Level: beginner

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

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

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

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

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


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

557:    Logically Collective on Vec

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

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

566:    Level: intermediate

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

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


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

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

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

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

610:    Logically Collective on Vec

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

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

619:    Level: intermediate

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


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

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

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

654:    Logically Collective on Vec

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

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

663:    Level: intermediate

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


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

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

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

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

705:    Logically Collective on Vec

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

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

714:    Level: intermediate

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


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

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

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


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

749:    Logically Collective on Vec

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

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

758:    Level: intermediate

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


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

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

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


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

798:    Not Collective

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

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

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

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

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

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

827:    Level: beginner

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

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

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

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

854:     Not Collective

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

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

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

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

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

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

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

877:    Level: beginner

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

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

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

898:    Not Collective

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

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

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

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

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

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

927:    Level: intermediate

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

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

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


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

955:    Not Collective

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

966:    Level: intermediate

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

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

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

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

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

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

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

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

1018:    Not Collective

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

1029:    Level: intermediate

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

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

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

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


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

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

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

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

1078:    Collective on Vec

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

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

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

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

1097:    Level: intermediate


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

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

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

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

1127:    Collective on Vec

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

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

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

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

1146:    Level: intermediate


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

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

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

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

1177:    Logically Collective on Vec

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

1185:    Level: intermediate

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

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

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

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

1224:    Collective on X and IS

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

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

1233:    Level: advanced

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

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

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

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

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

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

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

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

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

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

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

1361:    Collective on IS

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

1368:    Level: advanced

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

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

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

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

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

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

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

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

1454:    Not collective.

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

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

1462:    Level: beginner

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

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

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

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

1498: /*@
1499:    VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1500:    previously mapped into a vector using VecGetLocalVectorRead().

1502:    Not collective.

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

1508:    Level: beginner

1510: .seealso: VecGetLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1511: @*/
1512: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec w)
1513: {
1515:   PetscScalar    *a;

1520:   if (v->ops->restorelocalvectorread) {
1521:     (*v->ops->restorelocalvectorread)(v,w);
1522:   } else {
1523:     VecGetArrayRead(w,(const PetscScalar**)&a);
1524:     VecRestoreArrayRead(v,(const PetscScalar**)&a);
1525:     VecResetArray(w);
1526:   }
1527:   return(0);
1528: }

1530: /*@
1531:    VecGetLocalVector - Maps the local portion of a vector into a
1532:    vector.

1534:    Collective on v, not collective on w.

1536:    Input parameter:
1537: .  v - The vector for which the local vector is desired.

1539:    Output parameter:
1540: .  w - Upon exit this contains the local vector.

1542:    Level: beginner

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

1555: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1556: @*/
1557: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1558: {
1560:   PetscScalar    *a;

1565:   VecCheckSameLocalSize(v,1,w,2);
1566:   if (v->ops->getlocalvector) {
1567:     (*v->ops->getlocalvector)(v,w);
1568:   } else {
1569:     VecGetArray(v,&a);
1570:     VecPlaceArray(w,a);
1571:   }
1572:   return(0);
1573: }

1575: /*@
1576:    VecRestoreLocalVector - Unmaps the local portion of a vector
1577:    previously mapped into a vector using VecGetLocalVector().

1579:    Logically collective.

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

1585:    Level: beginner

1587: .seealso: VecGetLocalVector(), VecGetLocalVectorRead(), VecRestoreLocalVectorRead(), LocalVectorRead(), VecGetArrayRead(), VecGetArray()
1588: @*/
1589: PetscErrorCode VecRestoreLocalVector(Vec v,Vec w)
1590: {
1592:   PetscScalar    *a;

1597:   if (v->ops->restorelocalvector) {
1598:     (*v->ops->restorelocalvector)(v,w);
1599:   } else {
1600:     VecGetArray(w,&a);
1601:     VecRestoreArray(v,&a);
1602:     VecResetArray(w);
1603:   }
1604:   return(0);
1605: }

1607: /*@C
1608:    VecGetArray - Returns a pointer to a contiguous array that contains this
1609:    processor's portion of the vector data. For the standard PETSc
1610:    vectors, VecGetArray() returns a pointer to the local data array and
1611:    does not use any copies. If the underlying vector data is not stored
1612:    in a contiguous array this routine will copy the data to a contiguous
1613:    array and return a pointer to that. You MUST call VecRestoreArray()
1614:    when you no longer need access to the array.

1616:    Logically Collective on Vec

1618:    Input Parameter:
1619: .  x - the vector

1621:    Output Parameter:
1622: .  a - location to put pointer to the array

1624:    Fortran Note:
1625:    This routine is used differently from Fortran 77
1626: $    Vec         x
1627: $    PetscScalar x_array(1)
1628: $    PetscOffset i_x
1629: $    PetscErrorCode ierr
1630: $       call VecGetArray(x,x_array,i_x,ierr)
1631: $
1632: $   Access first local entry in vector with
1633: $      value = x_array(i_x + 1)
1634: $
1635: $      ...... other code
1636: $       call VecRestoreArray(x,x_array,i_x,ierr)
1637:    For Fortran 90 see VecGetArrayF90()

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

1642:    Level: beginner

1644: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1645:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1646: @*/
1647: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1648: {
1650: #if defined(PETSC_HAVE_VIENNACL)
1651:   PetscBool      is_viennacltype = PETSC_FALSE;
1652: #endif

1656:   VecSetErrorIfLocked(x,1);
1657:   if (x->petscnative) {
1658: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
1659:     if (x->offloadmask == PETSC_OFFLOAD_VECKOKKOS) { /* offloadmask here works as a tag quickly saying this is a VecKokkos */
1660:       VecKokkosSyncHost(x);
1661:       *a   = *((PetscScalar**)x->data);
1662:       return(0);
1663:     }
1664: #endif
1665: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1666:     if (x->offloadmask == PETSC_OFFLOAD_GPU) {
1667:   #if defined(PETSC_HAVE_VIENNACL)
1668:       PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1669:       if (is_viennacltype) {
1670:         VecViennaCLCopyFromGPU(x);
1671:       } else
1672:   #endif
1673:       {
1674:   #if defined(PETSC_HAVE_CUDA)
1675:         VecCUDACopyFromGPU(x);
1676:   #endif
1677:       }
1678:     } else if (x->offloadmask == PETSC_OFFLOAD_UNALLOCATED) {
1679:   #if defined(PETSC_HAVE_VIENNACL)
1680:       PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1681:       if (is_viennacltype) {
1682:         VecViennaCLAllocateCheckHost(x);
1683:       } else
1684:   #endif
1685:       {
1686:   #if defined(PETSC_HAVE_CUDA)
1687:         VecCUDAAllocateCheckHost(x);
1688:   #endif
1689:       }
1690:     }
1691: #endif
1692:     *a = *((PetscScalar**)x->data);
1693:   } else {
1694:     if (x->ops->getarray) {
1695:       (*x->ops->getarray)(x,a);
1696:     } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array for vector type \"%s\"",((PetscObject)x)->type_name);
1697:   }
1698:   return(0);
1699: }

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

1706:    Logically Collective on Vec

1708:    Input Parameter:
1709: .  x - the vector

1711:    Output Parameter:
1712: .  a - location to put pointer to the array

1714:    Level: beginner

1716: .seealso: VecRestoreArrayInPlace(), VecRestoreArrayInPlace(), VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(),
1717:           VecPlaceArray(), VecGetArray2d(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1718: @*/
1719: PetscErrorCode VecGetArrayInPlace(Vec x,PetscScalar **a)
1720: {
1723:   VecGetArrayInPlace_Internal(x,a,NULL);
1724:   return(0);
1725: }

1727: PetscErrorCode VecGetArrayInPlace_Internal(Vec x,PetscScalar **a,PetscMemType *mtype)
1728: {

1733:   VecSetErrorIfLocked(x,1);
1734:   if (mtype) *mtype = PETSC_MEMTYPE_HOST;
1735: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
1736:   if (x->offloadmask == PETSC_OFFLOAD_VECKOKKOS) {
1737:     VecKokkosGetArrayInPlace_Internal(x,a,mtype);
1738:     return(0);
1739:   }
1740: #endif

1742: #if defined(PETSC_HAVE_CUDA)
1743:   if (x->petscnative && (x->offloadmask & PETSC_OFFLOAD_GPU)) { /* Prefer working on GPU when offloadmask is PETSC_OFFLOAD_BOTH */
1744:     PetscBool is_cudatype = PETSC_FALSE;
1745:     PetscObjectTypeCompareAny((PetscObject)x,&is_cudatype,VECSEQCUDA,VECMPICUDA,VECCUDA,"");
1746:     if (is_cudatype) {
1747:       VecCUDAGetArray(x,a);
1748:       x->offloadmask = PETSC_OFFLOAD_GPU; /* Change the mask once GPU gets write access, don't wait until restore array */
1749:       if (mtype) *mtype = PETSC_MEMTYPE_DEVICE;
1750:       return(0);
1751:     }
1752:   }
1753: #endif
1754:   VecGetArray(x,a);
1755:   return(0);
1756: }

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

1763:    Logically Collective on Vec

1765:    Input Parameter:
1766: .  x - the vector

1768:    Output Parameter:
1769: .  a - location to put pointer to the array

1771:    Level: intermediate

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

1776:    Concepts: vector^accessing local values

1778: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1779:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArray(), VecRestoreArrayWrite()
1780: @*/
1781: PetscErrorCode VecGetArrayWrite(Vec x,PetscScalar **a)
1782: {

1787:   VecSetErrorIfLocked(x,1);
1788:   if (!x->ops->getarraywrite) {
1789:     VecGetArray(x,a);
1790:   } else {
1791:     (*x->ops->getarraywrite)(x,a);
1792:   }
1793:   return(0);
1794: }

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

1799:    Not Collective

1801:    Input Parameters:
1802: .  x - the vector

1804:    Output Parameter:
1805: .  a - the array

1807:    Level: beginner

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

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

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

1818: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1819: @*/
1820: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1821: {
1823: #if defined(PETSC_HAVE_VIENNACL)
1824:   PetscBool      is_viennacltype = PETSC_FALSE;
1825: #endif

1829:   if (x->petscnative) {
1830: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
1831:     if (x->offloadmask == PETSC_OFFLOAD_VECKOKKOS) {
1832:       VecKokkosSyncHost(x);
1833:       *a   = *((PetscScalar **)x->data);
1834:       return(0);
1835:     }
1836: #endif
1837: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1838:     if (x->offloadmask == PETSC_OFFLOAD_GPU) {
1839: #if defined(PETSC_HAVE_VIENNACL)
1840:       PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1841:       if (is_viennacltype) {
1842:         VecViennaCLCopyFromGPU(x);
1843:       } else
1844: #endif
1845:       {
1846: #if defined(PETSC_HAVE_CUDA)
1847:         VecCUDACopyFromGPU(x);
1848: #endif
1849:       }
1850:     }
1851: #endif
1852:     *a = *((PetscScalar **)x->data);
1853:   } else if (x->ops->getarrayread) {
1854:     (*x->ops->getarrayread)(x,a);
1855:   } else {
1856:     (*x->ops->getarray)(x,(PetscScalar**)a);
1857:   }
1858:   return(0);
1859: }

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

1866:    Not Collective

1868:    Input Parameters:
1869: .  x - the vector

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

1874:    Level: beginner

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


1880: .seealso: VecRestoreArrayReadInPlace(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayInPlace()
1881: @*/
1882: PetscErrorCode VecGetArrayReadInPlace(Vec x,const PetscScalar **a)
1883: {
1886:   VecGetArrayReadInPlace_Internal(x,a,NULL);
1887:   return(0);
1888: }

1890: PetscErrorCode VecGetArrayReadInPlace_Internal(Vec x,const PetscScalar **a,PetscMemType *mtype)
1891: {

1896:   if (mtype) *mtype = PETSC_MEMTYPE_HOST;
1897: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
1898:   if (x->offloadmask == PETSC_OFFLOAD_VECKOKKOS) {
1899:     VecKokkosGetArrayReadInPlace_Internal(x,a,mtype);
1900:     return(0);
1901:   }
1902: #endif

1904: #if defined(PETSC_HAVE_CUDA)
1905:   if (x->petscnative && x->offloadmask & PETSC_OFFLOAD_GPU) {
1906:     PetscBool is_cudatype = PETSC_FALSE;
1907:     PetscObjectTypeCompareAny((PetscObject)x,&is_cudatype,VECSEQCUDA,VECMPICUDA,VECCUDA,"");
1908:     if (is_cudatype) {
1909:       VecCUDAGetArrayRead(x,a);
1910:       if (mtype) *mtype = PETSC_MEMTYPE_DEVICE;
1911:       return(0);
1912:     }
1913:   }
1914: #endif
1915:   VecGetArrayRead(x,a);
1916:   return(0);
1917: }

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

1924:    Logically Collective on Vec

1926:    Input Parameter:
1927: +  x - the vectors
1928: -  n - the number of vectors

1930:    Output Parameter:
1931: .  a - location to put pointer to the array

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

1936:    Level: intermediate

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

1950:   if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1951:   PetscMalloc1(n,&q);
1952:   for (i=0; i<n; ++i) {
1953:     VecGetArray(x[i],&q[i]);
1954:   }
1955:   *a = q;
1956:   return(0);
1957: }

1959: /*@C
1960:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1961:    has been called.

1963:    Logically Collective on Vec

1965:    Input Parameters:
1966: +  x - the vector
1967: .  n - the number of vectors
1968: -  a - location of pointer to arrays obtained from VecGetArrays()

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

1976:    Fortran Note:
1977:    This routine is not supported in Fortran.

1979:    Level: intermediate

1981: .seealso: VecGetArrays(), VecRestoreArray()
1982: @*/
1983: PetscErrorCode  VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1984: {
1986:   PetscInt       i;
1987:   PetscScalar    **q = *a;


1994:   for (i=0; i<n; ++i) {
1995:     VecRestoreArray(x[i],&q[i]);
1996:   }
1997:   PetscFree(q);
1998:   return(0);
1999: }

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

2004:    Logically Collective on Vec

2006:    Input Parameters:
2007: +  x - the vector
2008: -  a - location of pointer to array obtained from VecGetArray()

2010:    Level: beginner

2012: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
2013:           VecGetArrayPair(), VecRestoreArrayPair()
2014: @*/
2015: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
2016: {

2021:   if (x->petscnative) {
2022:    #if defined(PETSC_HAVE_KOKKOS_KERNELS)
2023:     if (x->offloadmask == PETSC_OFFLOAD_VECKOKKOS) {VecKokkosModifyHost(x);}
2024:     else
2025:    #endif
2026:     {
2027:      #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
2028:       x->offloadmask = PETSC_OFFLOAD_CPU;
2029:      #endif
2030:     }
2031:   } else {
2032:     (*x->ops->restorearray)(x,a);
2033:   }
2034:   if (a) *a = NULL;
2035:   PetscObjectStateIncrease((PetscObject)x);
2036:   return(0);
2037: }

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

2042:    Logically Collective on Vec

2044:    Input Parameters:
2045: +  x - the vector
2046: -  a - location of pointer to array obtained from VecGetArrayInPlace()

2048:    Level: beginner

2050: .seealso: VecGetArrayInPlace(), VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(),
2051:           VecPlaceArray(), VecRestoreArray2d(), VecGetArrayPair(), VecRestoreArrayPair()
2052: @*/
2053: PetscErrorCode VecRestoreArrayInPlace(Vec x,PetscScalar **a)
2054: {

2059: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
2060:   if (x->offloadmask == PETSC_OFFLOAD_VECKOKKOS) {
2061:     VecKokkosRestoreArrayInPlace(x,a);
2062:     return(0);
2063:   }
2064: #endif

2066: #if defined(PETSC_HAVE_CUDA)
2067:   if (x->petscnative && x->offloadmask == PETSC_OFFLOAD_GPU) {
2068:     PetscBool is_cudatype = PETSC_FALSE;
2069:     PetscObjectTypeCompareAny((PetscObject)x,&is_cudatype,VECSEQCUDA,VECMPICUDA,VECCUDA,"");
2070:     if (is_cudatype) {
2071:       VecCUDARestoreArray(x,a);
2072:       return(0);
2073:     }
2074:   }
2075: #endif
2076:   VecRestoreArray(x,a);
2077:   return(0);
2078: }


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

2084:    Logically Collective on Vec

2086:    Input Parameters:
2087: +  x - the vector
2088: -  a - location of pointer to array obtained from VecGetArray()

2090:    Level: beginner

2092: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
2093:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite()
2094: @*/
2095: PetscErrorCode VecRestoreArrayWrite(Vec x,PetscScalar **a)
2096: {


2102: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
2103:   if (x->offloadmask == PETSC_OFFLOAD_VECKOKKOS) {
2104:     VecKokkosModifyHost(x);
2105:   } else
2106: #endif
2107:   if (x->petscnative) {
2108: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
2109:     x->offloadmask = PETSC_OFFLOAD_CPU;
2110: #endif
2111:   } else {
2112:     if (x->ops->restorearraywrite) {
2113:       (*x->ops->restorearraywrite)(x,a);
2114:     } else {
2115:       (*x->ops->restorearray)(x,a);
2116:     }
2117:   }

2119:   if (a) *a = NULL;
2120:   PetscObjectStateIncrease((PetscObject)x);
2121:   return(0);
2122: }

2124: /*@C
2125:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

2127:    Not Collective

2129:    Input Parameters:
2130: +  vec - the vector
2131: -  array - the array

2133:    Level: beginner

2135: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
2136: @*/
2137: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
2138: {

2143:   if (x->petscnative) {
2144:     /* nothing */
2145:   } else if (x->ops->restorearrayread) {
2146:     (*x->ops->restorearrayread)(x,a);
2147:   } else {
2148:     (*x->ops->restorearray)(x,(PetscScalar**)a);
2149:   }
2150:   if (a) *a = NULL;
2151:   return(0);
2152: }

2154: /*@C
2155:    VecRestoreArrayReadInPlace - Restore array obtained with VecGetArrayReadInPlace()

2157:    Not Collective

2159:    Input Parameters:
2160: +  vec - the vector
2161: -  array - the array

2163:    Level: beginner

2165: .seealso: VecGetArrayReadInPlace(), VecRestoreArrayInPlace(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
2166: @*/
2167: PetscErrorCode VecRestoreArrayReadInPlace(Vec x,const PetscScalar **a)
2168: {

2172:   VecRestoreArrayRead(x,a);
2173:   return(0);
2174: }

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

2181:    Not Collective

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

2187:    Notes:
2188:    You can return to the original array with a call to VecResetArray()

2190:    Level: developer

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

2194: @*/
2195: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
2196: {

2203:   if (vec->ops->placearray) {
2204:     (*vec->ops->placearray)(vec,array);
2205:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
2206:   PetscObjectStateIncrease((PetscObject)vec);
2207:   return(0);
2208: }

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

2215:    Not Collective

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

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

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

2228:    Not supported from Fortran

2230:    Level: developer

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

2234: @*/
2235: PetscErrorCode  VecReplaceArray(Vec vec,const PetscScalar array[])
2236: {

2242:   if (vec->ops->replacearray) {
2243:     (*vec->ops->replacearray)(vec,array);
2244:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
2245:   PetscObjectStateIncrease((PetscObject)vec);
2246:   return(0);
2247: }


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

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

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


2269:    Input Parameter:
2270: .  v - the vector

2272:    Output Parameter:
2273: .  a - the CUDA device pointer

2275:    Fortran note:
2276:    This function is not currently available from Fortran.

2278:    Level: intermediate

2280: .seealso: VecCUDARestoreArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2281: @*/
2282: PETSC_EXTERN PetscErrorCode VecCUDAGetArray(Vec v, PetscScalar **a)
2283: {
2284: #if defined(PETSC_HAVE_CUDA)
2286: #endif

2290: #if defined(PETSC_HAVE_CUDA)
2291:   *a   = 0;
2292:   VecCUDACopyToGPU(v);
2293:   *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2294: #endif
2295:   return(0);
2296: }

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

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

2305:    Input Parameter:
2306: +  v - the vector
2307: -  a - the CUDA device pointer.  This pointer is invalid after
2308:        VecCUDARestoreArray() returns.

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

2313:    Level: intermediate

2315: .seealso: VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2316: @*/
2317: PETSC_EXTERN PetscErrorCode VecCUDARestoreArray(Vec v, PetscScalar **a)
2318: {

2323: #if defined(PETSC_HAVE_CUDA)
2324:   v->offloadmask = PETSC_OFFLOAD_GPU;
2325: #endif

2327:   PetscObjectStateIncrease((PetscObject)v);
2328:   return(0);
2329: }

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

2334:    This function is analogous to VecGetArrayRead():  The pointer
2335:    returned by this function points to a consistent view of the vector
2336:    data.  This may involve a copy operation of data from the host to the
2337:    device if the data on the device is out of date.  If the device
2338:    memory hasn't been allocated previously it will be allocated as part
2339:    of this function call.  VecCUDAGetArrayRead() assumes that the
2340:    user will not modify the vector data.  This is analgogous to
2341:    intent(in) in Fortran.

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

2349:    Input Parameter:
2350: .  v - the vector

2352:    Output Parameter:
2353: .  a - the CUDA pointer.

2355:    Fortran note:
2356:    This function is not currently available from Fortran.

2358:    Level: intermediate

2360: .seealso: VecCUDARestoreArrayRead(), VecCUDAGetArray(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2361: @*/
2362: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayRead(Vec v, const PetscScalar **a)
2363: {
2364: #if defined(PETSC_HAVE_CUDA)
2366: #endif

2370: #if defined(PETSC_HAVE_CUDA)
2371:   *a   = 0;
2372:   VecCUDACopyToGPU(v);
2373:   *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2374: #endif
2375:   return(0);
2376: }

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

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

2386:    Input Parameter:
2387: +  v - the vector
2388: -  a - the CUDA device pointer.  This pointer is invalid after
2389:        VecCUDARestoreArrayRead() returns.

2391:    Fortran note:
2392:    This function is not currently available from Fortran.

2394:    Level: intermediate

2396: .seealso: VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2397: @*/
2398: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayRead(Vec v, const PetscScalar **a)
2399: {
2402:   *a = NULL;
2403:   return(0);
2404: }

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

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

2415:    The device pointer needs to be released with
2416:    VecCUDARestoreArrayWrite().  When the pointer is released the
2417:    host data of the vector is marked as out of data.  Subsequent access
2418:    of the host data with e.g. VecGetArray() incurs a device to host data
2419:    transfer.


2422:    Input Parameter:
2423: .  v - the vector

2425:    Output Parameter:
2426: .  a - the CUDA pointer

2428:    Fortran note:
2429:    This function is not currently available from Fortran.

2431:    Level: advanced

2433: .seealso: VecCUDARestoreArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2434: @*/
2435: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayWrite(Vec v, PetscScalar **a)
2436: {
2437: #if defined(PETSC_HAVE_CUDA)
2439: #endif

2443: #if defined(PETSC_HAVE_CUDA)
2444:   *a   = 0;
2445:   VecCUDAAllocateCheck(v);
2446:   *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2447: #endif
2448:   return(0);
2449: }

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

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

2458:    Input Parameter:
2459: +  v - the vector
2460: -  a - the CUDA device pointer.  This pointer is invalid after
2461:        VecCUDARestoreArrayWrite() returns.

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

2466:    Level: intermediate

2468: .seealso: VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2469: @*/
2470: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayWrite(Vec v, PetscScalar **a)
2471: {

2476: #if defined(PETSC_HAVE_CUDA)
2477:   v->offloadmask = PETSC_OFFLOAD_GPU;
2478: #endif

2480:   PetscObjectStateIncrease((PetscObject)v);
2481:   return(0);
2482: }

2484: /*@C
2485:    VecCUDAPlaceArray - Allows one to replace the GPU array in a vector with a
2486:    GPU array provided by the user. This is useful to avoid copying an
2487:    array into a vector.

2489:    Not Collective

2491:    Input Parameters:
2492: +  vec - the vector
2493: -  array - the GPU array

2495:    Notes:
2496:    You can return to the original GPU array with a call to VecCUDAResetArray()
2497:    It is not possible to use VecCUDAPlaceArray() and VecPlaceArray() at the
2498:    same time on the same vector.

2500:    Level: developer

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

2504: @*/
2505: PetscErrorCode VecCUDAPlaceArray(Vec vin,const PetscScalar a[])
2506: {

2511: #if defined(PETSC_HAVE_CUDA)
2512:   VecCUDACopyToGPU(vin);
2513:   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()");
2514:   ((Vec_Seq*)vin->data)->unplacedarray  = (PetscScalar *) ((Vec_CUDA*)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2515:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2516:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2517: #endif
2518:   PetscObjectStateIncrease((PetscObject)vin);
2519:   return(0);
2520: }

2522: /*@C
2523:    VecCUDAReplaceArray - Allows one to replace the GPU array in a vector
2524:    with a GPU array provided by the user. This is useful to avoid copying
2525:    a GPU array into a vector.

2527:    Not Collective

2529:    Input Parameters:
2530: +  vec - the vector
2531: -  array - the GPU array

2533:    Notes:
2534:    This permanently replaces the GPU array and frees the memory associated
2535:    with the old GPU array.

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

2540:    Not supported from Fortran

2542:    Level: developer

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

2546: @*/
2547: PetscErrorCode VecCUDAReplaceArray(Vec vin,const PetscScalar a[])
2548: {
2549: #if defined(PETSC_HAVE_CUDA)
2550:   cudaError_t err;
2551: #endif

2556: #if defined(PETSC_HAVE_CUDA)
2557:   err = cudaFree(((Vec_CUDA*)vin->spptr)->GPUarray);CHKERRCUDA(err);
2558:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2559:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2560: #endif
2561:   PetscObjectStateIncrease((PetscObject)vin);
2562:   return(0);
2563: }

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

2569:    Not Collective

2571:    Input Parameters:
2572: .  vec - the vector

2574:    Level: developer

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

2578: @*/
2579: PetscErrorCode VecCUDAResetArray(Vec vin)
2580: {

2585: #if defined(PETSC_HAVE_CUDA)
2586:   VecCUDACopyToGPU(vin);
2587:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar *) ((Vec_Seq*)vin->data)->unplacedarray;
2588:   ((Vec_Seq*)vin->data)->unplacedarray = 0;
2589:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2590: #endif
2591:   PetscObjectStateIncrease((PetscObject)vin);
2592:   return(0);
2593: }

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

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

2602:     Collective on Vec

2604:     Input Parameters:
2605: +   x - a vector to mimic
2606: -   n - the number of vectors to obtain

2608:     Output Parameters:
2609: +   y - Fortran90 pointer to the array of vectors
2610: -   ierr - error code

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

2617:     Vec x
2618:     Vec, pointer :: y(:)
2619:     ....
2620:     call VecDuplicateVecsF90(x,2,y,ierr)
2621:     call VecSet(y(2),alpha,ierr)
2622:     call VecSet(y(2),alpha,ierr)
2623:     ....
2624:     call VecDestroyVecsF90(2,y,ierr)
2625: .ve

2627:     Notes:
2628:     Not yet supported for all F90 compilers

2630:     Use VecDestroyVecsF90() to free the space.

2632:     Level: beginner

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

2636: M*/

2638: /*MC
2639:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2640:     VecGetArrayF90().

2642:     Synopsis:
2643:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2645:     Logically Collective on Vec

2647:     Input Parameters:
2648: +   x - vector
2649: -   xx_v - the Fortran90 pointer to the array

2651:     Output Parameter:
2652: .   ierr - error code

2654:     Example of Usage:
2655: .vb
2656: #include <petsc/finclude/petscvec.h>
2657:     use petscvec

2659:     PetscScalar, pointer :: xx_v(:)
2660:     ....
2661:     call VecGetArrayF90(x,xx_v,ierr)
2662:     xx_v(3) = a
2663:     call VecRestoreArrayF90(x,xx_v,ierr)
2664: .ve

2666:     Level: beginner

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

2670: M*/

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

2675:     Synopsis:
2676:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

2678:     Collective on Vec

2680:     Input Parameters:
2681: +   n - the number of vectors previously obtained
2682: -   x - pointer to array of vector pointers

2684:     Output Parameter:
2685: .   ierr - error code

2687:     Notes:
2688:     Not yet supported for all F90 compilers

2690:     Level: beginner

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

2694: M*/

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

2702:     Synopsis:
2703:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2705:     Logically Collective on Vec

2707:     Input Parameter:
2708: .   x - vector

2710:     Output Parameters:
2711: +   xx_v - the Fortran90 pointer to the array
2712: -   ierr - error code

2714:     Example of Usage:
2715: .vb
2716: #include <petsc/finclude/petscvec.h>
2717:     use petscvec

2719:     PetscScalar, pointer :: xx_v(:)
2720:     ....
2721:     call VecGetArrayF90(x,xx_v,ierr)
2722:     xx_v(3) = a
2723:     call VecRestoreArrayF90(x,xx_v,ierr)
2724: .ve

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

2728:     Level: beginner

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

2732: M*/

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

2740:     Synopsis:
2741:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2743:     Logically Collective on Vec

2745:     Input Parameter:
2746: .   x - vector

2748:     Output Parameters:
2749: +   xx_v - the Fortran90 pointer to the array
2750: -   ierr - error code

2752:     Example of Usage:
2753: .vb
2754: #include <petsc/finclude/petscvec.h>
2755:     use petscvec

2757:     PetscScalar, pointer :: xx_v(:)
2758:     ....
2759:     call VecGetArrayReadF90(x,xx_v,ierr)
2760:     a = xx_v(3)
2761:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2762: .ve

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

2766:     Level: beginner

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

2770: M*/

2772: /*MC
2773:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2774:     VecGetArrayReadF90().

2776:     Synopsis:
2777:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2779:     Logically Collective on Vec

2781:     Input Parameters:
2782: +   x - vector
2783: -   xx_v - the Fortran90 pointer to the array

2785:     Output Parameter:
2786: .   ierr - error code

2788:     Example of Usage:
2789: .vb
2790: #include <petsc/finclude/petscvec.h>
2791:     use petscvec

2793:     PetscScalar, pointer :: xx_v(:)
2794:     ....
2795:     call VecGetArrayReadF90(x,xx_v,ierr)
2796:     a = xx_v(3)
2797:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2798: .ve

2800:     Level: beginner

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

2804: M*/

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

2811:    Logically Collective

2813:    Input Parameter:
2814: +  x - the vector
2815: .  m - first dimension of two dimensional array
2816: .  n - second dimension of two dimensional array
2817: .  mstart - first index you will use in first coordinate direction (often 0)
2818: -  nstart - first index in the second coordinate direction (often 0)

2820:    Output Parameter:
2821: .  a - location to put pointer to the array

2823:    Level: developer

2825:   Notes:
2826:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2827:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2828:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2829:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

2833: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2834:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2835:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2836: @*/
2837: PetscErrorCode  VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2838: {
2840:   PetscInt       i,N;
2841:   PetscScalar    *aa;

2847:   VecGetLocalSize(x,&N);
2848:   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);
2849:   VecGetArray(x,&aa);

2851:   PetscMalloc1(m,a);
2852:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2853:   *a -= mstart;
2854:   return(0);
2855: }

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

2862:    Logically Collective

2864:    Input Parameter:
2865: +  x - the vector
2866: .  m - first dimension of two dimensional array
2867: .  n - second dimension of two dimensional array
2868: .  mstart - first index you will use in first coordinate direction (often 0)
2869: -  nstart - first index in the second coordinate direction (often 0)

2871:    Output Parameter:
2872: .  a - location to put pointer to the array

2874:    Level: developer

2876:   Notes:
2877:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2878:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2879:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2880:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

2886: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2887:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2888:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2889: @*/
2890: PetscErrorCode  VecGetArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2891: {
2893:   PetscInt       i,N;
2894:   PetscScalar    *aa;

2900:   VecGetLocalSize(x,&N);
2901:   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);
2902:   VecGetArrayWrite(x,&aa);

2904:   PetscMalloc1(m,a);
2905:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2906:   *a -= mstart;
2907:   return(0);
2908: }

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

2913:    Logically Collective

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

2923:    Level: developer

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

2931:    This routine actually zeros out the a pointer.

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

2946:   dummy = (void*)(*a + mstart);
2947:   PetscFree(dummy);
2948:   VecRestoreArray(x,NULL);
2949:   return(0);
2950: }

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

2955:    Logically Collective

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

2965:    Level: developer

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

2973:    This routine actually zeros out the a pointer.

2975: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2976:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2977:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2978: @*/
2979: PetscErrorCode  VecRestoreArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2980: {
2982:   void           *dummy;

2988:   dummy = (void*)(*a + mstart);
2989:   PetscFree(dummy);
2990:   VecRestoreArrayWrite(x,NULL);
2991:   return(0);
2992: }

2994: /*@C
2995:    VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2996:    processor's portion of the vector data.  You MUST call VecRestoreArray1d()
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  VecGetArray1d(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:   VecGetArray(x,a);
3034:   *a  -= mstart;
3035:   return(0);
3036: }

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

3043:    Logically Collective

3045:    Input Parameter:
3046: +  x - the vector
3047: .  m - first dimension of two dimensional array
3048: -  mstart - first index you will use in first coordinate direction (often 0)

3050:    Output Parameter:
3051: .  a - location to put pointer to the array

3053:    Level: developer

3055:   Notes:
3056:    For a vector obtained from DMCreateLocalVector() mstart are likely
3057:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3058:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3062: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3063:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3064:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3065: @*/
3066: PetscErrorCode  VecGetArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3067: {
3069:   PetscInt       N;

3075:   VecGetLocalSize(x,&N);
3076:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3077:   VecGetArrayWrite(x,a);
3078:   *a  -= mstart;
3079:   return(0);
3080: }

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

3085:    Logically Collective

3087:    Input Parameters:
3088: +  x - the vector
3089: .  m - first dimension of two dimensional array
3090: .  mstart - first index you will use in first coordinate direction (often 0)
3091: -  a - location of pointer to array obtained from VecGetArray21()

3093:    Level: developer

3095:    Notes:
3096:    For regular PETSc vectors this routine does not involve any copies. For
3097:    any special vectors that do not store local vector data in a contiguous
3098:    array, this routine will copy the data back into the underlying
3099:    vector data structure from the array obtained with VecGetArray1d().

3101:    This routine actually zeros out the a pointer.

3103: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3104:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3105:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3106: @*/
3107: PetscErrorCode  VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3108: {

3114:   VecRestoreArray(x,NULL);
3115:   return(0);
3116: }

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

3121:    Logically Collective

3123:    Input Parameters:
3124: +  x - the vector
3125: .  m - first dimension of two dimensional array
3126: .  mstart - first index you will use in first coordinate direction (often 0)
3127: -  a - location of pointer to array obtained from VecGetArray21()

3129:    Level: developer

3131:    Notes:
3132:    For regular PETSc vectors this routine does not involve any copies. For
3133:    any special vectors that do not store local vector data in a contiguous
3134:    array, this routine will copy the data back into the underlying
3135:    vector data structure from the array obtained with VecGetArray1d().

3137:    This routine actually zeros out the a pointer.

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

3141: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3142:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3143:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3144: @*/
3145: PetscErrorCode  VecRestoreArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3146: {

3152:   VecRestoreArrayWrite(x,NULL);
3153:   return(0);
3154: }

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

3161:    Logically Collective

3163:    Input Parameter:
3164: +  x - the vector
3165: .  m - first dimension of three dimensional array
3166: .  n - second dimension of three dimensional array
3167: .  p - third dimension of three dimensional array
3168: .  mstart - first index you will use in first coordinate direction (often 0)
3169: .  nstart - first index in the second coordinate direction (often 0)
3170: -  pstart - first index in the third coordinate direction (often 0)

3172:    Output Parameter:
3173: .  a - location to put pointer to the array

3175:    Level: developer

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

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

3185: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3186:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3187:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3188: @*/
3189: PetscErrorCode  VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3190: {
3192:   PetscInt       i,N,j;
3193:   PetscScalar    *aa,**b;

3199:   VecGetLocalSize(x,&N);
3200:   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);
3201:   VecGetArray(x,&aa);

3203:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3204:   b    = (PetscScalar**)((*a) + m);
3205:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3206:   for (i=0; i<m; i++)
3207:     for (j=0; j<n; j++)
3208:       b[i*n+j] = aa + i*n*p + j*p - pstart;

3210:   *a -= mstart;
3211:   return(0);
3212: }

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

3219:    Logically Collective

3221:    Input Parameter:
3222: +  x - the vector
3223: .  m - first dimension of three dimensional array
3224: .  n - second dimension of three dimensional array
3225: .  p - third dimension of three dimensional array
3226: .  mstart - first index you will use in first coordinate direction (often 0)
3227: .  nstart - first index in the second coordinate direction (often 0)
3228: -  pstart - first index in the third coordinate direction (often 0)

3230:    Output Parameter:
3231: .  a - location to put pointer to the array

3233:    Level: developer

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

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

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

3245: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3246:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3247:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3248: @*/
3249: PetscErrorCode  VecGetArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3250: {
3252:   PetscInt       i,N,j;
3253:   PetscScalar    *aa,**b;

3259:   VecGetLocalSize(x,&N);
3260:   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);
3261:   VecGetArrayWrite(x,&aa);

3263:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3264:   b    = (PetscScalar**)((*a) + m);
3265:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3266:   for (i=0; i<m; i++)
3267:     for (j=0; j<n; j++)
3268:       b[i*n+j] = aa + i*n*p + j*p - pstart;

3270:   *a -= mstart;
3271:   return(0);
3272: }

3274: /*@C
3275:    VecRestoreArray3d - Restores a vector after VecGetArray3d() 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  VecRestoreArray3d(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:   VecRestoreArray(x,NULL);
3315:   return(0);
3316: }

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

3321:    Logically Collective

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

3333:    Level: developer

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

3341:    This routine actually zeros out the a pointer.

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

3356:   dummy = (void*)(*a + mstart);
3357:   PetscFree(dummy);
3358:   VecRestoreArrayWrite(x,NULL);
3359:   return(0);
3360: }

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

3367:    Logically Collective

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

3380:    Output Parameter:
3381: .  a - location to put pointer to the array

3383:    Level: beginner

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

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

3393: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3394:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3395:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3396: @*/
3397: PetscErrorCode  VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3398: {
3400:   PetscInt       i,N,j,k;
3401:   PetscScalar    *aa,***b,**c;

3407:   VecGetLocalSize(x,&N);
3408:   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);
3409:   VecGetArray(x,&aa);

3411:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3412:   b    = (PetscScalar***)((*a) + m);
3413:   c    = (PetscScalar**)(b + m*n);
3414:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3415:   for (i=0; i<m; i++)
3416:     for (j=0; j<n; j++)
3417:       b[i*n+j] = c + i*n*p + j*p - pstart;
3418:   for (i=0; i<m; i++)
3419:     for (j=0; j<n; j++)
3420:       for (k=0; k<p; k++)
3421:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3422:   *a -= mstart;
3423:   return(0);
3424: }

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

3431:    Logically Collective

3433:    Input Parameter:
3434: +  x - the vector
3435: .  m - first dimension of four dimensional array
3436: .  n - second dimension of four dimensional array
3437: .  p - third dimension of four dimensional array
3438: .  q - fourth dimension of four dimensional array
3439: .  mstart - first index you will use in first coordinate direction (often 0)
3440: .  nstart - first index in the second coordinate direction (often 0)
3441: .  pstart - first index in the third coordinate direction (often 0)
3442: -  qstart - first index in the fourth coordinate direction (often 0)

3444:    Output Parameter:
3445: .  a - location to put pointer to the array

3447:    Level: beginner

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

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

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

3459: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3460:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3461:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3462: @*/
3463: PetscErrorCode  VecGetArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3464: {
3466:   PetscInt       i,N,j,k;
3467:   PetscScalar    *aa,***b,**c;

3473:   VecGetLocalSize(x,&N);
3474:   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);
3475:   VecGetArrayWrite(x,&aa);

3477:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3478:   b    = (PetscScalar***)((*a) + m);
3479:   c    = (PetscScalar**)(b + m*n);
3480:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3481:   for (i=0; i<m; i++)
3482:     for (j=0; j<n; j++)
3483:       b[i*n+j] = c + i*n*p + j*p - pstart;
3484:   for (i=0; i<m; i++)
3485:     for (j=0; j<n; j++)
3486:       for (k=0; k<p; k++)
3487:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3488:   *a -= mstart;
3489:   return(0);
3490: }

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

3495:    Logically Collective

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

3509:    Level: beginner

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

3517:    This routine actually zeros out the a pointer.

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

3532:   dummy = (void*)(*a + mstart);
3533:   PetscFree(dummy);
3534:   VecRestoreArray(x,NULL);
3535:   return(0);
3536: }

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

3541:    Logically Collective

3543:    Input Parameters:
3544: +  x - the vector
3545: .  m - first dimension of four dimensional array
3546: .  n - second dimension of the four dimensional array
3547: .  p - third dimension of the four dimensional array
3548: .  q - fourth dimension of the four dimensional array
3549: .  mstart - first index you will use in first coordinate direction (often 0)
3550: .  nstart - first index in the second coordinate direction (often 0)
3551: .  pstart - first index in the third coordinate direction (often 0)
3552: .  qstart - first index in the fourth coordinate direction (often 0)
3553: -  a - location of pointer to array obtained from VecGetArray4d()

3555:    Level: beginner

3557:    Notes:
3558:    For regular PETSc vectors this routine does not involve any copies. For
3559:    any special vectors that do not store local vector data in a contiguous
3560:    array, this routine will copy the data back into the underlying
3561:    vector data structure from the array obtained with VecGetArray().

3563:    This routine actually zeros out the a pointer.

3565: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3566:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3567:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3568: @*/
3569: PetscErrorCode  VecRestoreArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3570: {
3572:   void           *dummy;

3578:   dummy = (void*)(*a + mstart);
3579:   PetscFree(dummy);
3580:   VecRestoreArrayWrite(x,NULL);
3581:   return(0);
3582: }

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

3589:    Logically Collective

3591:    Input Parameter:
3592: +  x - the vector
3593: .  m - first dimension of two dimensional array
3594: .  n - second dimension of two dimensional array
3595: .  mstart - first index you will use in first coordinate direction (often 0)
3596: -  nstart - first index in the second coordinate direction (often 0)

3598:    Output Parameter:
3599: .  a - location to put pointer to the array

3601:    Level: developer

3603:   Notes:
3604:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
3605:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3606:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3607:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

3611: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3612:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3613:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3614: @*/
3615: PetscErrorCode  VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3616: {
3617:   PetscErrorCode    ierr;
3618:   PetscInt          i,N;
3619:   const PetscScalar *aa;

3625:   VecGetLocalSize(x,&N);
3626:   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);
3627:   VecGetArrayRead(x,&aa);

3629:   PetscMalloc1(m,a);
3630:   for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
3631:   *a -= mstart;
3632:   return(0);
3633: }

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

3638:    Logically Collective

3640:    Input Parameters:
3641: +  x - the vector
3642: .  m - first dimension of two dimensional array
3643: .  n - second dimension of the two dimensional array
3644: .  mstart - first index you will use in first coordinate direction (often 0)
3645: .  nstart - first index in the second coordinate direction (often 0)
3646: -  a - location of pointer to array obtained from VecGetArray2d()

3648:    Level: developer

3650:    Notes:
3651:    For regular PETSc vectors this routine does not involve any copies. For
3652:    any special vectors that do not store local vector data in a contiguous
3653:    array, this routine will copy the data back into the underlying
3654:    vector data structure from the array obtained with VecGetArray().

3656:    This routine actually zeros out the a pointer.

3658: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3659:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3660:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3661: @*/
3662: PetscErrorCode  VecRestoreArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3663: {
3665:   void           *dummy;

3671:   dummy = (void*)(*a + mstart);
3672:   PetscFree(dummy);
3673:   VecRestoreArrayRead(x,NULL);
3674:   return(0);
3675: }

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

3682:    Logically Collective

3684:    Input Parameter:
3685: +  x - the vector
3686: .  m - first dimension of two dimensional array
3687: -  mstart - first index you will use in first coordinate direction (often 0)

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

3692:    Level: developer

3694:   Notes:
3695:    For a vector obtained from DMCreateLocalVector() mstart are likely
3696:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3697:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3701: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3702:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3703:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3704: @*/
3705: PetscErrorCode  VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3706: {
3708:   PetscInt       N;

3714:   VecGetLocalSize(x,&N);
3715:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3716:   VecGetArrayRead(x,(const PetscScalar**)a);
3717:   *a  -= mstart;
3718:   return(0);
3719: }

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

3724:    Logically Collective

3726:    Input Parameters:
3727: +  x - the vector
3728: .  m - first dimension of two dimensional array
3729: .  mstart - first index you will use in first coordinate direction (often 0)
3730: -  a - location of pointer to array obtained from VecGetArray21()

3732:    Level: developer

3734:    Notes:
3735:    For regular PETSc vectors this routine does not involve any copies. For
3736:    any special vectors that do not store local vector data in a contiguous
3737:    array, this routine will copy the data back into the underlying
3738:    vector data structure from the array obtained with VecGetArray1dRead().

3740:    This routine actually zeros out the a pointer.

3742: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3743:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3744:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3745: @*/
3746: PetscErrorCode  VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3747: {

3753:   VecRestoreArrayRead(x,NULL);
3754:   return(0);
3755: }


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

3763:    Logically Collective

3765:    Input Parameter:
3766: +  x - the vector
3767: .  m - first dimension of three dimensional array
3768: .  n - second dimension of three dimensional array
3769: .  p - third dimension of three dimensional array
3770: .  mstart - first index you will use in first coordinate direction (often 0)
3771: .  nstart - first index in the second coordinate direction (often 0)
3772: -  pstart - first index in the third coordinate direction (often 0)

3774:    Output Parameter:
3775: .  a - location to put pointer to the array

3777:    Level: developer

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

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

3787: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3788:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3789:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3790: @*/
3791: PetscErrorCode  VecGetArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3792: {
3793:   PetscErrorCode    ierr;
3794:   PetscInt          i,N,j;
3795:   const PetscScalar *aa;
3796:   PetscScalar       **b;

3802:   VecGetLocalSize(x,&N);
3803:   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);
3804:   VecGetArrayRead(x,&aa);

3806:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3807:   b    = (PetscScalar**)((*a) + m);
3808:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3809:   for (i=0; i<m; i++)
3810:     for (j=0; j<n; j++)
3811:       b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;

3813:   *a -= mstart;
3814:   return(0);
3815: }

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

3820:    Logically Collective

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

3832:    Level: developer

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

3840:    This routine actually zeros out the a pointer.

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

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

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

3866:    Logically Collective

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

3879:    Output Parameter:
3880: .  a - location to put pointer to the array

3882:    Level: beginner

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

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

3892: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3893:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3894:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3895: @*/
3896: PetscErrorCode  VecGetArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3897: {
3898:   PetscErrorCode    ierr;
3899:   PetscInt          i,N,j,k;
3900:   const PetscScalar *aa;
3901:   PetscScalar       ***b,**c;

3907:   VecGetLocalSize(x,&N);
3908:   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);
3909:   VecGetArrayRead(x,&aa);

3911:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3912:   b    = (PetscScalar***)((*a) + m);
3913:   c    = (PetscScalar**)(b + m*n);
3914:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3915:   for (i=0; i<m; i++)
3916:     for (j=0; j<n; j++)
3917:       b[i*n+j] = c + i*n*p + j*p - pstart;
3918:   for (i=0; i<m; i++)
3919:     for (j=0; j<n; j++)
3920:       for (k=0; k<p; k++)
3921:         c[i*n*p+j*p+k] = (PetscScalar*) aa + i*n*p*q + j*p*q + k*q - qstart;
3922:   *a -= mstart;
3923:   return(0);
3924: }

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

3929:    Logically Collective

3931:    Input Parameters:
3932: +  x - the vector
3933: .  m - first dimension of four dimensional array
3934: .  n - second dimension of the four dimensional array
3935: .  p - third dimension of the four dimensional array
3936: .  q - fourth dimension of the four dimensional array
3937: .  mstart - first index you will use in first coordinate direction (often 0)
3938: .  nstart - first index in the second coordinate direction (often 0)
3939: .  pstart - first index in the third coordinate direction (often 0)
3940: .  qstart - first index in the fourth coordinate direction (often 0)
3941: -  a - location of pointer to array obtained from VecGetArray4dRead()

3943:    Level: beginner

3945:    Notes:
3946:    For regular PETSc vectors this routine does not involve any copies. For
3947:    any special vectors that do not store local vector data in a contiguous
3948:    array, this routine will copy the data back into the underlying
3949:    vector data structure from the array obtained with VecGetArray().

3951:    This routine actually zeros out the a pointer.

3953: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3954:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3955:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3956: @*/
3957: PetscErrorCode  VecRestoreArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3958: {
3960:   void           *dummy;

3966:   dummy = (void*)(*a + mstart);
3967:   PetscFree(dummy);
3968:   VecRestoreArrayRead(x,NULL);
3969:   return(0);
3970: }

3972: #if defined(PETSC_USE_DEBUG)

3974: /*@
3975:    VecLockGet  - Gets the current lock status of a vector

3977:    Logically Collective on Vec

3979:    Input Parameter:
3980: .  x - the vector

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

3986:    Level: beginner

3988: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop()
3989: @*/
3990: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
3991: {
3994:   *state = x->lock;
3995:   return(0);
3996: }

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

4001:    Logically Collective on Vec

4003:    Input Parameter:
4004: .  x - the vector

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

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

4012:    Level: beginner

4014: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPop(), VecLockGet()
4015: @*/
4016: PetscErrorCode VecLockReadPush(Vec x)
4017: {
4020:   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");
4021:   x->lock++;
4022:   return(0);
4023: }

4025: /*@
4026:    VecLockReadPop  - Pops a read-only lock from a vector

4028:    Logically Collective on Vec

4030:    Input Parameter:
4031: .  x - the vector

4033:    Level: beginner

4035: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockGet()
4036: @*/
4037: PetscErrorCode VecLockReadPop(Vec x)
4038: {
4041:   x->lock--;
4042:   if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector has been unlocked from read-only access too many times");
4043:   return(0);
4044: }

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

4049:    Logically Collective on Vec

4051:    Input Parameter:
4052: +  x   - the vector
4053: -  flg - PETSC_TRUE to lock the vector for writing; PETSC_FALSE to unlock it.

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


4062:        VecGetArray(x,&xdata); // begin phase
4063:        VecLockWriteSet_Private(v,PETSC_TRUE);

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

4067:        VecRestoreArray(x,&vdata); // end phase
4068:        VecLockWriteSet_Private(v,PETSC_FALSE);

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

4073:    Level: beginner

4075: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop(), VecLockGet()
4076: @*/
4077: PetscErrorCode VecLockWriteSet_Private(Vec x,PetscBool flg)
4078: {
4081:   if (flg) {
4082:     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");
4083:     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");
4084:     else x->lock = -1;
4085:   } else {
4086:     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");
4087:     x->lock = 0;
4088:   }
4089:   return(0);
4090: }

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

4095:    Level: deprecated

4097: .seealso: VecLockReadPush()
4098: @*/
4099: PetscErrorCode VecLockPush(Vec x)
4100: {
4103:   VecLockReadPush(x);
4104:   return(0);
4105: }

4107: /*@
4108:    VecLockPop  - Pops a read-only lock from a vector

4110:    Level: deprecated

4112: .seealso: VecLockReadPop()
4113: @*/
4114: PetscErrorCode VecLockPop(Vec x)
4115: {
4118:   VecLockReadPop(x);
4119:   return(0);
4120: }

4122: #endif