Actual source code: rvector.c

  1: /*
  2:      Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
  3:    These are the vector functions the user calls.
  4: */
  5: #include "petsc/private/sfimpl.h"
  6: #include "petscsystypes.h"
  7: #include <petsc/private/vecimpl.h>
  8: #if defined(PETSC_HAVE_CUDA)
  9: #include <../src/vec/vec/impls/dvecimpl.h>
 10: #include <petsc/private/cudavecimpl.h>
 11: #endif
 12: #if defined(PETSC_HAVE_HIP)
 13: #include <../src/vec/vec/impls/dvecimpl.h>
 14: #include <petsc/private/hipvecimpl.h>
 15: #endif
 16: static PetscInt VecGetSubVectorSavedStateId = -1;

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

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

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

 51:    Logically Collective on Vec

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

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

 59:    Level: advanced

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

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

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

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

 86:    Collective on Vec

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

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

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

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

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

110:    Level: intermediate

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

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

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

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

136:    Collective on Vec

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

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

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

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

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

154:      For complex numbers in C^n (that is a vector of n components with a complex number for each component) this is equal to the usual real dot product on the
155:      the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)

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

159:    Level: intermediate

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

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

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

177:    Collective on Vec

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

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

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

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

198:    Level: intermediate

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

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

208: @*/

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


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

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

239:    Not Collective

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

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

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

256:    Level: intermediate

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

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

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

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


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

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

291:    Collective on Vec

293:    Input Parameter:
294: .  x - the vector

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

299:    Level: intermediate

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

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

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

326:    Collective on Vec

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

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

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

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

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

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

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

360:    Collective on Vec

362:    Input Parameter:
363: .  x - the vector

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

369:    Level: intermediate

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

374:    This returns the smallest index with the minumum value

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

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

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

396:    Collective on Vec

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

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

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

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

413:    Level: intermediate

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

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

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

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

439:    Not collective on Vec

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

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

449:    Level: intermediate

451: @*/
452: PetscErrorCode  VecScale(Vec x, PetscScalar alpha)
453: {
454:   PetscReal      norms[4] = {0.0,0.0,0.0, 0.0};
455:   PetscBool      flgs[4];
457:   PetscInt       i;

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

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

486:    Logically Collective on Vec

488:    Input Parameters:
489: +  x  - the vector
490: -  alpha - the scalar

492:    Output Parameter:
493: .  x  - the vector

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

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

505:    Level: beginner

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

509: @*/
510: PetscErrorCode  VecSet(Vec x,PetscScalar alpha)
511: {
512:   PetscReal      val;

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

522:   PetscLogEventBegin(VEC_Set,x,0,0,0);
523:   (*x->ops->set)(x,alpha);
524:   PetscLogEventEnd(VEC_Set,x,0,0,0);
525:   PetscObjectStateIncrease((PetscObject)x);

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

546: /*@
547:    VecAXPY - Computes y = alpha x + y.

549:    Logically Collective on Vec

551:    Input Parameters:
552: +  alpha - the scalar
553: -  x, y  - the vectors

555:    Output Parameter:
556: .  y - output vector

558:    Level: intermediate

560:    Notes:
561:     x and y MUST be different vectors
562:     This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine

564: $    VecAXPY(y,alpha,x)                   y = alpha x           +      y
565: $    VecAYPX(y,beta,x)                    y =       x           + beta y
566: $    VecAXPBY(y,alpha,beta,x)             y = alpha x           + beta y
567: $    VecWAXPY(w,alpha,x,y)                w = alpha x           +      y
568: $    VecAXPBYPCZ(w,alpha,beta,gamma,x,y)  z = alpha x           + beta y + gamma z
569: $    VecMAXPY(y,nv,alpha[],x[])           y = sum alpha[i] x[i] +      y

571: .seealso:  VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPBYPCZ(), VecAXPBY()
572: @*/
573: PetscErrorCode  VecAXPY(Vec y,PetscScalar alpha,Vec x)
574: {

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

589:   VecLockReadPush(x);
590:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
591:   (*y->ops->axpy)(y,alpha,x);
592:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
593:   VecLockReadPop(x);
594:   PetscObjectStateIncrease((PetscObject)y);
595:   return(0);
596: }

598: /*@
599:    VecAXPBY - Computes y = alpha x + beta y.

601:    Logically Collective on Vec

603:    Input Parameters:
604: +  alpha,beta - the scalars
605: -  x, y  - the vectors

607:    Output Parameter:
608: .  y - output vector

610:    Level: intermediate

612:    Notes:
613:     x and y MUST be different vectors
614:     The implementation is optimized for alpha and/or beta values of 0.0 and 1.0

616: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ()
617: @*/
618: PetscErrorCode  VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
619: {

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

641: /*@
642:    VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z

644:    Logically Collective on Vec

646:    Input Parameters:
647: +  alpha,beta, gamma - the scalars
648: -  x, y, z  - the vectors

650:    Output Parameter:
651: .  z - output vector

653:    Level: intermediate

655:    Notes:
656:     x, y and z must be different vectors
657:     The implementation is optimized for alpha of 1.0 and gamma of 1.0 or 0.0

659: .seealso:  VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBY()
660: @*/
661: PetscErrorCode  VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
662: {

674:   VecCheckSameSize(x,1,y,5);
675:   VecCheckSameSize(x,1,z,6);
676:   if (x == y || x == z) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
677:   if (y == z) SETERRQ(PetscObjectComm((PetscObject)y),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
681:   if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) return(0);
682:   VecSetErrorIfLocked(z,1);

684:   PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
685:   (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
686:   PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
687:   PetscObjectStateIncrease((PetscObject)z);
688:   return(0);
689: }

691: /*@
692:    VecAYPX - Computes y = x + beta y.

694:    Logically Collective on Vec

696:    Input Parameters:
697: +  beta - the scalar
698: -  x, y  - the vectors

700:    Output Parameter:
701: .  y - output vector

703:    Level: intermediate

705:    Notes:
706:     x and y MUST be different vectors
707:     The implementation is optimized for beta of -1.0, 0.0, and 1.0

709: .seealso:  VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
710: @*/
711: PetscErrorCode  VecAYPX(Vec y,PetscScalar beta,Vec x)
712: {

721:   VecCheckSameSize(x,1,y,3);
722:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
724:   VecSetErrorIfLocked(y,1);

726:   PetscLogEventBegin(VEC_AYPX,x,y,0,0);
727:    (*y->ops->aypx)(y,beta,x);
728:   PetscLogEventEnd(VEC_AYPX,x,y,0,0);
729:   PetscObjectStateIncrease((PetscObject)y);
730:   return(0);
731: }

733: /*@
734:    VecWAXPY - Computes w = alpha x + y.

736:    Logically Collective on Vec

738:    Input Parameters:
739: +  alpha - the scalar
740: -  x, y  - the vectors

742:    Output Parameter:
743: .  w - the result

745:    Level: intermediate

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

751: .seealso: VecAXPY(), VecAYPX(), VecAXPBY(), VecMAXPY(), VecAXPBYPCZ()
752: @*/
753: PetscErrorCode  VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
754: {

766:   VecCheckSameSize(x,3,y,4);
767:   VecCheckSameSize(x,3,w,1);
768:   if (w == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector y, suggest VecAXPY()");
769:   if (w == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector x, suggest VecAYPX()");
771:   VecSetErrorIfLocked(w,1);

773:   PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
774:    (*w->ops->waxpy)(w,alpha,x,y);
775:   PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
776:   PetscObjectStateIncrease((PetscObject)w);
777:   return(0);
778: }

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

783:    Not Collective

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

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

797:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
798:    options cannot be mixed without intervening calls to the assembly
799:    routines.

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

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

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

812:    Level: beginner

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

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

828:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
829:   (*x->ops->setvalues)(x,ni,ix,y,iora);
830:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
831:   PetscObjectStateIncrease((PetscObject)x);
832:   return(0);
833: }

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

839:     Not Collective

841:    Input Parameters:
842: +  x - vector to get values from
843: .  ni - number of elements to get
844: -  ix - indices where to get them from (in global 1d numbering)

846:    Output Parameter:
847: .   y - array of values

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

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

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

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

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

862:    Level: beginner

864: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues()
865: @*/
866: PetscErrorCode  VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
867: {

872:   if (!ni) return(0);
876:   (*x->ops->getvalues)(x,ni,ix,y);
877:   return(0);
878: }

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

883:    Not Collective

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

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

898:    Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
899:    options cannot be mixed without intervening calls to the assembly
900:    routines.

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

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

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

912:    Level: intermediate

914: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
915:            VecSetValues()
916: @*/
917: PetscErrorCode  VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
918: {

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

928:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
929:   (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
930:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
931:   PetscObjectStateIncrease((PetscObject)x);
932:   return(0);
933: }

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

939:    Not Collective

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

950:    Level: intermediate

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

955:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
956:    options cannot be mixed without intervening calls to the assembly
957:    routines.

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

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

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

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

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

998: /*@
999:    VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1000:    using a local ordering of the nodes.

1002:    Not Collective

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

1013:    Level: intermediate

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

1019:    Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1020:    options cannot be mixed without intervening calls to the assembly
1021:    routines.

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

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

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

1038:   if (!ni) return(0);
1042:   if (ni > 128) {
1043:     PetscMalloc1(ni,&lix);
1044:   }

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

1057: /*@
1058:    VecMTDot - Computes indefinite vector multiple dot products.
1059:    That is, it does NOT use the complex conjugate.

1061:    Collective on Vec

1063:    Input Parameters:
1064: +  x - one vector
1065: .  nv - number of vectors
1066: -  y - array of vectors.  Note that vectors are pointers

1068:    Output Parameter:
1069: .  val - array of the dot products

1071:    Notes for Users of Complex Numbers:
1072:    For complex vectors, VecMTDot() computes the indefinite form
1073: $      val = (x,y) = y^T x,
1074:    where y^T denotes the transpose of y.

1076:    Use VecMDot() for the inner product
1077: $      val = (x,y) = y^H x,
1078:    where y^H denotes the conjugate transpose of y.

1080:    Level: intermediate

1082: .seealso: VecMDot(), VecTDot()
1083: @*/
1084: PetscErrorCode  VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1085: {

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

1100:   PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1101:   (*x->ops->mtdot)(x,nv,y,val);
1102:   PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1103:   return(0);
1104: }

1106: /*@
1107:    VecMDot - Computes vector multiple dot products.

1109:    Collective on Vec

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

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

1119:    Notes for Users of Complex Numbers:
1120:    For complex vectors, VecMDot() computes
1121: $     val = (x,y) = y^H x,
1122:    where y^H denotes the conjugate transpose of y.

1124:    Use VecMTDot() for the indefinite form
1125: $     val = (x,y) = y^T x,
1126:    where y^T denotes the transpose of y.

1128:    Level: intermediate

1130: .seealso: VecMTDot(), VecDot()
1131: @*/
1132: PetscErrorCode  VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1133: {

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

1149:   PetscLogEventBegin(VEC_MDot,x,*y,0,0);
1150:   (*x->ops->mdot)(x,nv,y,val);
1151:   PetscLogEventEnd(VEC_MDot,x,*y,0,0);
1152:   return(0);
1153: }

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

1158:    Logically Collective on Vec

1160:    Input Parameters:
1161: +  nv - number of scalars and x-vectors
1162: .  alpha - array of scalars
1163: .  y - one vector
1164: -  x - array of vectors

1166:    Level: intermediate

1168:    Notes:
1169:     y cannot be any of the x vectors

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

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

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

1207:    Collective on X

1209:    Input Parameters:
1210: +  nx   - number of vectors to be concatenated
1211: -  X    - array containing the vectors to be concatenated in the order of concatenation

1213:    Output Parameters:
1214: +  Y    - concatenated vector
1215: -  x_is - array of index sets corresponding to the concatenated components of Y (NULL if not needed)

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

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

1227:    Level: advanced

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


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

1285: /*@
1286:    VecGetSubVector - Gets a vector representing part of another vector

1288:    Collective on X and IS

1290:    Input Parameters:
1291: + X - vector from which to extract a subvector
1292: - is - index set representing portion of X to extract

1294:    Output Parameter:
1295: . Y - subvector corresponding to is

1297:    Level: advanced

1299:    Notes:
1300:    The subvector Y should be returned with VecRestoreSubVector().
1301:    X and is must be defined on the same communicator

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

1307: .seealso: MatCreateSubMatrix()
1308: @*/
1309: PetscErrorCode  VecGetSubVector(Vec X,IS is,Vec *Y)
1310: {
1311:   PetscErrorCode   ierr;
1312:   Vec              Z;

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

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

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

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

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

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

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

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

1437: /*@
1438:    VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()

1440:    Collective on IS

1442:    Input Parameters:
1443: + X - vector from which subvector was obtained
1444: . is - index set representing the subset of X
1445: - Y - subvector being restored

1447:    Level: advanced

1449: .seealso: VecGetSubVector()
1450: @*/
1451: PetscErrorCode  VecRestoreSubVector(Vec X,IS is,Vec *Y)
1452: {

1461:   if (X->ops->restoresubvector) {
1462:     (*X->ops->restoresubvector)(X,is,Y);
1463:   } else {
1464:     PETSC_UNUSED PetscObjectState dummystate = 0;
1465:     PetscBool valid;

1467:     PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,dummystate,valid);
1468:     if (!valid) {
1469:       VecScatter scatter;
1470:       PetscInt   state;

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

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

1484:         if (iscuda) {
1485: #if defined(PETSC_HAVE_CUDA)
1486:           PetscOffloadMask ymask = (*Y)->offloadmask;

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

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

1556: /*@
1557:    VecGetLocalVectorRead - Maps the local portion of a vector into a
1558:    vector.  You must call VecRestoreLocalVectorRead() when the local
1559:    vector is no longer needed.

1561:    Not collective.

1563:    Input parameter:
1564: .  v - The vector for which the local vector is desired.

1566:    Output parameter:
1567: .  w - Upon exit this contains the local vector.

1569:    Level: beginner

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

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

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

1595:   VecCheckSameLocalSize(v,1,w,2);
1596:   if (v->ops->getlocalvectorread) {
1597:     (*v->ops->getlocalvectorread)(v,w);
1598:   } else {
1599:     VecGetArrayRead(v,(const PetscScalar**)&a);
1600:     VecPlaceArray(w,a);
1601:   }
1602:   PetscObjectStateIncrease((PetscObject)w);
1603:   VecLockReadPush(v);
1604:   VecLockReadPush(w);
1605:   return(0);
1606: }

1608: /*@
1609:    VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1610:    previously mapped into a vector using VecGetLocalVectorRead().

1612:    Not collective.

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

1618:    Level: beginner

1620: .seealso: VecGetLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1621: @*/
1622: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec w)
1623: {
1625:   PetscScalar    *a;

1630:   if (v->ops->restorelocalvectorread) {
1631:     (*v->ops->restorelocalvectorread)(v,w);
1632:   } else {
1633:     VecGetArrayRead(w,(const PetscScalar**)&a);
1634:     VecRestoreArrayRead(v,(const PetscScalar**)&a);
1635:     VecResetArray(w);
1636:   }
1637:   VecLockReadPop(v);
1638:   VecLockReadPop(w);
1639:   PetscObjectStateIncrease((PetscObject)w);
1640:   return(0);
1641: }

1643: /*@
1644:    VecGetLocalVector - Maps the local portion of a vector into a
1645:    vector.

1647:    Collective on v, not collective on w.

1649:    Input parameter:
1650: .  v - The vector for which the local vector is desired.

1652:    Output parameter:
1653: .  w - Upon exit this contains the local vector.

1655:    Level: beginner

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

1668: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1669: @*/
1670: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1671: {
1673:   PetscScalar    *a;

1678:   VecCheckSameLocalSize(v,1,w,2);
1679:   if (v->ops->getlocalvector) {
1680:     (*v->ops->getlocalvector)(v,w);
1681:   } else {
1682:     VecGetArray(v,&a);
1683:     VecPlaceArray(w,a);
1684:   }
1685:   PetscObjectStateIncrease((PetscObject)w);
1686:   return(0);
1687: }

1689: /*@
1690:    VecRestoreLocalVector - Unmaps the local portion of a vector
1691:    previously mapped into a vector using VecGetLocalVector().

1693:    Logically collective.

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

1699:    Level: beginner

1701: .seealso: VecGetLocalVector(), VecGetLocalVectorRead(), VecRestoreLocalVectorRead(), LocalVectorRead(), VecGetArrayRead(), VecGetArray()
1702: @*/
1703: PetscErrorCode VecRestoreLocalVector(Vec v,Vec w)
1704: {
1706:   PetscScalar    *a;

1711:   if (v->ops->restorelocalvector) {
1712:     (*v->ops->restorelocalvector)(v,w);
1713:   } else {
1714:     VecGetArray(w,&a);
1715:     VecRestoreArray(v,&a);
1716:     VecResetArray(w);
1717:   }
1718:   PetscObjectStateIncrease((PetscObject)w);
1719:   PetscObjectStateIncrease((PetscObject)v);
1720:   return(0);
1721: }

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

1732:    Logically Collective on Vec

1734:    Input Parameter:
1735: .  x - the vector

1737:    Output Parameter:
1738: .  a - location to put pointer to the array

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

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

1758:    Level: beginner

1760: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1761:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1762: @*/
1763: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1764: {

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

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

1781:    Logically Collective on Vec

1783:    Input Parameters:
1784: +  x - the vector
1785: -  a - location of pointer to array obtained from VecGetArray()

1787:    Level: beginner

1789: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1790:           VecGetArrayPair(), VecRestoreArrayPair()
1791: @*/
1792: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1793: {

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

1810:    Not Collective

1812:    Input Parameter:
1813: .  x - the vector

1815:    Output Parameter:
1816: .  a - the array

1818:    Level: beginner

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

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

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

1829: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1830: @*/
1831: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1832: {

1837:   if (x->ops->getarray) { /* VECNEST, VECCUDA, VECKOKKOS etc */
1838:     (*x->ops->getarray)(x,(PetscScalar**)a);
1839:   } else if (x->petscnative) { /* VECSTANDARD */
1840:     *a = *((PetscScalar**)x->data);
1841:   } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array read for vector type \"%s\"",((PetscObject)x)->type_name);
1842:   return(0);
1843: }

1845: /*@C
1846:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

1848:    Not Collective

1850:    Input Parameters:
1851: +  vec - the vector
1852: -  array - the array

1854:    Level: beginner

1856: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1857: @*/
1858: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1859: {

1864:   if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
1865:     /* nothing */
1866:   } else if (x->ops->restorearrayread) { /* VECNEST */
1867:     (*x->ops->restorearrayread)(x,a);
1868:   } else { /* No one? */
1869:     (*x->ops->restorearray)(x,(PetscScalar**)a);
1870:   }
1871:   if (a) *a = NULL;
1872:   return(0);
1873: }

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

1880:    Logically Collective on Vec

1882:    Input Parameter:
1883: .  x - the vector

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

1888:    Level: intermediate

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

1893:    Concepts: vector^accessing local values

1895: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1896:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArray(), VecRestoreArrayWrite()
1897: @*/
1898: PetscErrorCode VecGetArrayWrite(Vec x,PetscScalar **a)
1899: {

1904:   VecSetErrorIfLocked(x,1);
1905:   if (x->ops->getarraywrite) {
1906:     (*x->ops->getarraywrite)(x,a);
1907:   } else {
1908:     VecGetArray(x,a);
1909:   }
1910:   return(0);
1911: }

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

1916:    Logically Collective on Vec

1918:    Input Parameters:
1919: +  x - the vector
1920: -  a - location of pointer to array obtained from VecGetArray()

1922:    Level: beginner

1924: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1925:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite()
1926: @*/
1927: PetscErrorCode VecRestoreArrayWrite(Vec x,PetscScalar **a)
1928: {

1933:   if (x->ops->restorearraywrite) {
1934:     (*x->ops->restorearraywrite)(x,a);
1935:   } else if (x->ops->restorearray) {
1936:     (*x->ops->restorearray)(x,a);
1937:   }
1938:   if (a) *a = NULL;
1939:   PetscObjectStateIncrease((PetscObject)x);
1940:   return(0);
1941: }

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

1948:    Logically Collective on Vec

1950:    Input Parameters:
1951: +  x - the vectors
1952: -  n - the number of vectors

1954:    Output Parameter:
1955: .  a - location to put pointer to the array

1957:    Fortran Note:
1958:    This routine is not supported in Fortran.

1960:    Level: intermediate

1962: .seealso: VecGetArray(), VecRestoreArrays()
1963: @*/
1964: PetscErrorCode  VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1965: {
1967:   PetscInt       i;
1968:   PetscScalar    **q;

1974:   if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1975:   PetscMalloc1(n,&q);
1976:   for (i=0; i<n; ++i) {
1977:     VecGetArray(x[i],&q[i]);
1978:   }
1979:   *a = q;
1980:   return(0);
1981: }

1983: /*@C
1984:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1985:    has been called.

1987:    Logically Collective on Vec

1989:    Input Parameters:
1990: +  x - the vector
1991: .  n - the number of vectors
1992: -  a - location of pointer to arrays obtained from VecGetArrays()

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

2000:    Fortran Note:
2001:    This routine is not supported in Fortran.

2003:    Level: intermediate

2005: .seealso: VecGetArrays(), VecRestoreArray()
2006: @*/
2007: PetscErrorCode  VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
2008: {
2010:   PetscInt       i;
2011:   PetscScalar    **q = *a;


2018:   for (i=0; i<n; ++i) {
2019:     VecRestoreArray(x[i],&q[i]);
2020:   }
2021:   PetscFree(q);
2022:   return(0);
2023: }

2025: /*@C
2026:    VecGetArrayAndMemType - Like VecGetArray(), but if this is a device vector (e.g., VECCUDA) and the device has up-to-date data,
2027:    the returned pointer will be a device pointer to the device memory that contains this processor's portion of the vector data.
2028:    Otherwise, when this is a host vector (e.g., VECMPI), or a device vector with the host having newer data than device,
2029:    it functions as VecGetArray() and returns a host pointer.

2031:    Logically Collective on Vec

2033:    Input Parameter:
2034: .  x - the vector

2036:    Output Parameters:
2037: +  a - location to put pointer to the array
2038: -  mtype - memory type of the array

2040:    Level: beginner

2042: .seealso: VecRestoreArrayAndMemType(), VecRestoreArrayAndMemType(), VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(),
2043:           VecPlaceArray(), VecGetArray2d(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
2044: @*/
2045: PetscErrorCode VecGetArrayAndMemType(Vec x,PetscScalar **a,PetscMemType *mtype)
2046: {

2052:   VecSetErrorIfLocked(x,1);
2053:   if (x->ops->getarrayandmemtype) { /* VECCUDA, VECKOKKOS etc */
2054:     (*x->ops->getarrayandmemtype)(x,a,mtype);
2055:   } else { /* VECSTANDARD, VECNEST, VECVIENNACL */
2056:     VecGetArray(x,a);
2057:     if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2058:   }
2059:   return(0);
2060: }

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

2065:    Logically Collective on Vec

2067:    Input Parameters:
2068: +  x - the vector
2069: -  a - location of pointer to array obtained from VecGetArrayAndMemType()

2071:    Level: beginner

2073: .seealso: VecGetArrayAndMemType(), VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(),
2074:           VecPlaceArray(), VecRestoreArray2d(), VecGetArrayPair(), VecRestoreArrayPair()
2075: @*/
2076: PetscErrorCode VecRestoreArrayAndMemType(Vec x,PetscScalar **a)
2077: {

2083:   if (x->ops->restorearrayandmemtype) { /* VECCUDA, VECKOKKOS etc */
2084:     (*x->ops->restorearrayandmemtype)(x,a);
2085:   } else if (x->ops->restorearray) { /* VECNEST, VECVIENNACL */
2086:     (*x->ops->restorearray)(x,a);
2087:   } /* VECSTANDARD does nothing */
2088:   if (a) *a = NULL;
2089:   PetscObjectStateIncrease((PetscObject)x);
2090:   return(0);
2091: }

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

2098:    Not Collective

2100:    Input Parameter:
2101: .  x - the vector

2103:    Output Parameters:
2104: +  a - the array
2105: -  mtype - memory type of the array

2107:    Level: beginner

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

2112: .seealso: VecRestoreArrayReadAndMemType(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayAndMemType()
2113: @*/
2114: PetscErrorCode VecGetArrayReadAndMemType(Vec x,const PetscScalar **a,PetscMemType *mtype)
2115: {

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

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

2137: /*@C
2138:    VecRestoreArrayReadAndMemType - Restore array obtained with VecGetArrayReadAndMemType()

2140:    Not Collective

2142:    Input Parameters:
2143: +  vec - the vector
2144: -  array - the array

2146:    Level: beginner

2148: .seealso: VecGetArrayReadAndMemType(), VecRestoreArrayAndMemType(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
2149: @*/
2150: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x,const PetscScalar **a)
2151: {

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

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

2171:    Not Collective

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

2177:    Notes:
2178:    You can return to the original array with a call to VecResetArray()

2180:    Level: developer

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

2184: @*/
2185: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
2186: {

2193:   if (vec->ops->placearray) {
2194:     (*vec->ops->placearray)(vec,array);
2195:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
2196:   PetscObjectStateIncrease((PetscObject)vec);
2197:   return(0);
2198: }

2200: /*@C
2201:    VecReplaceArray - Allows one to replace the array in a vector with an
2202:    array provided by the user. This is useful to avoid copying an array
2203:    into a vector.

2205:    Not Collective

2207:    Input Parameters:
2208: +  vec - the vector
2209: -  array - the array

2211:    Notes:
2212:    This permanently replaces the array and frees the memory associated
2213:    with the old array.

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

2218:    Not supported from Fortran

2220:    Level: developer

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

2224: @*/
2225: PetscErrorCode  VecReplaceArray(Vec vec,const PetscScalar array[])
2226: {

2232:   if (vec->ops->replacearray) {
2233:     (*vec->ops->replacearray)(vec,array);
2234:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
2235:   PetscObjectStateIncrease((PetscObject)vec);
2236:   return(0);
2237: }

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

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

2251:    The CUDA device pointer has to be released by calling
2252:    VecCUDARestoreArray().  Upon restoring the vector data
2253:    the data on the host will be marked as out of date.  A subsequent
2254:    access of the host data will thus incur a data transfer from the
2255:    device to the host.

2257:    Input Parameter:
2258: .  v - the vector

2260:    Output Parameter:
2261: .  a - the CUDA device pointer

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

2266:    Level: intermediate

2268: .seealso: VecCUDARestoreArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2269: @*/
2270: PETSC_EXTERN PetscErrorCode VecCUDAGetArray(Vec v, PetscScalar **a)
2271: {
2274:  #if defined(PETSC_HAVE_CUDA)
2275:   {
2277:     VecCUDACopyToGPU(v);
2278:     *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2279:   }
2280:  #endif
2281:   return(0);
2282: }

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

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

2291:    Input Parameters:
2292: +  v - the vector
2293: -  a - the CUDA device pointer.  This pointer is invalid after
2294:        VecCUDARestoreArray() returns.

2296:    Fortran note:
2297:    This function is not currently available from Fortran.

2299:    Level: intermediate

2301: .seealso: VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2302: @*/
2303: PETSC_EXTERN PetscErrorCode VecCUDARestoreArray(Vec v, PetscScalar **a)
2304: {

2309: #if defined(PETSC_HAVE_CUDA)
2310:   v->offloadmask = PETSC_OFFLOAD_GPU;
2311: #endif
2312:   PetscObjectStateIncrease((PetscObject)v);
2313:   return(0);
2314: }

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

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

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

2334:    Input Parameter:
2335: .  v - the vector

2337:    Output Parameter:
2338: .  a - the CUDA pointer.

2340:    Fortran note:
2341:    This function is not currently available from Fortran.

2343:    Level: intermediate

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

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

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

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

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

2371:    Level: intermediate

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

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

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

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

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

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

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

2407:    Level: advanced

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

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

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

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

2437:    Fortran note:
2438:    This function is not currently available from Fortran.

2440:    Level: intermediate

2442: .seealso: VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2443: @*/
2444: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayWrite(Vec v, PetscScalar **a)
2445: {

2450:  #if defined(PETSC_HAVE_CUDA)
2451:   v->offloadmask = PETSC_OFFLOAD_GPU;
2452:   if (a) *a = NULL;
2453:  #endif
2454:   PetscObjectStateIncrease((PetscObject)v);
2455:   return(0);
2456: }

2458: /*@C
2459:    VecCUDAPlaceArray - Allows one to replace the GPU array in a vector with a
2460:    GPU array provided by the user. This is useful to avoid copying an
2461:    array into a vector.

2463:    Not Collective

2465:    Input Parameters:
2466: +  vec - the vector
2467: -  array - the GPU array

2469:    Notes:
2470:    You can return to the original GPU array with a call to VecCUDAResetArray()
2471:    It is not possible to use VecCUDAPlaceArray() and VecPlaceArray() at the
2472:    same time on the same vector.

2474:    Level: developer

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

2478: @*/
2479: PetscErrorCode VecCUDAPlaceArray(Vec vin,const PetscScalar a[])
2480: {

2485: #if defined(PETSC_HAVE_CUDA)
2486:   VecCUDACopyToGPU(vin);
2487:   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()");
2488:   ((Vec_Seq*)vin->data)->unplacedarray  = (PetscScalar *) ((Vec_CUDA*)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2489:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2490:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2491: #endif
2492:   PetscObjectStateIncrease((PetscObject)vin);
2493:   return(0);
2494: }

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

2501:    Not Collective

2503:    Input Parameters:
2504: +  vec - the vector
2505: -  array - the GPU array

2507:    Notes:
2508:    This permanently replaces the GPU array and frees the memory associated
2509:    with the old GPU array.

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

2514:    Not supported from Fortran

2516:    Level: developer

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

2520: @*/
2521: PetscErrorCode VecCUDAReplaceArray(Vec vin,const PetscScalar a[])
2522: {
2523: #if defined(PETSC_HAVE_CUDA)
2524:   cudaError_t err;
2525: #endif

2530: #if defined(PETSC_HAVE_CUDA)
2531:   if (((Vec_CUDA*)vin->spptr)->GPUarray_allocated) {
2532:     err = cudaFree(((Vec_CUDA*)vin->spptr)->GPUarray_allocated);CHKERRCUDA(err);
2533:   }
2534:   ((Vec_CUDA*)vin->spptr)->GPUarray_allocated = ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2535:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2536: #endif
2537:   PetscObjectStateIncrease((PetscObject)vin);
2538:   return(0);
2539: }

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

2545:    Not Collective

2547:    Input Parameters:
2548: .  vec - the vector

2550:    Level: developer

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

2554: @*/
2555: PetscErrorCode VecCUDAResetArray(Vec vin)
2556: {

2561: #if defined(PETSC_HAVE_CUDA)
2562:   VecCUDACopyToGPU(vin);
2563:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar *) ((Vec_Seq*)vin->data)->unplacedarray;
2564:   ((Vec_Seq*)vin->data)->unplacedarray = 0;
2565:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2566: #endif
2567:   PetscObjectStateIncrease((PetscObject)vin);
2568:   return(0);
2569: }

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

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

2583:    The HIP device pointer has to be released by calling
2584:    VecHIPRestoreArray().  Upon restoring the vector data
2585:    the data on the host will be marked as out of date.  A subsequent
2586:    access of the host data will thus incur a data transfer from the
2587:    device to the host.

2589:    Input Parameter:
2590: .  v - the vector

2592:    Output Parameter:
2593: .  a - the HIP device pointer

2595:    Fortran note:
2596:    This function is not currently available from Fortran.

2598:    Level: intermediate

2600: .seealso: VecHIPRestoreArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2601: @*/
2602: PETSC_EXTERN PetscErrorCode VecHIPGetArray(Vec v, PetscScalar **a)
2603: {
2604: #if defined(PETSC_HAVE_HIP)
2606: #endif

2610: #if defined(PETSC_HAVE_HIP)
2611:   *a   = 0;
2612:   VecHIPCopyToGPU(v);
2613:   *a   = ((Vec_HIP*)v->spptr)->GPUarray;
2614: #endif
2615:   return(0);
2616: }

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

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

2625:    Input Parameters:
2626: +  v - the vector
2627: -  a - the HIP device pointer.  This pointer is invalid after
2628:        VecHIPRestoreArray() returns.

2630:    Fortran note:
2631:    This function is not currently available from Fortran.

2633:    Level: intermediate

2635: .seealso: VecHIPGetArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2636: @*/
2637: PETSC_EXTERN PetscErrorCode VecHIPRestoreArray(Vec v, PetscScalar **a)
2638: {

2643: #if defined(PETSC_HAVE_HIP)
2644:   v->offloadmask = PETSC_OFFLOAD_GPU;
2645: #endif

2647:   PetscObjectStateIncrease((PetscObject)v);
2648:   return(0);
2649: }

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

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

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

2669:    Input Parameter:
2670: .  v - the vector

2672:    Output Parameter:
2673: .  a - the HIP pointer.

2675:    Fortran note:
2676:    This function is not currently available from Fortran.

2678:    Level: intermediate

2680: .seealso: VecHIPRestoreArrayRead(), VecHIPGetArray(), VecHIPGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2681: @*/
2682: PETSC_EXTERN PetscErrorCode VecHIPGetArrayRead(Vec v, const PetscScalar **a)
2683: {
2684: #if defined(PETSC_HAVE_HIP)
2686: #endif

2690: #if defined(PETSC_HAVE_HIP)
2691:   *a   = 0;
2692:   VecHIPCopyToGPU(v);
2693:   *a   = ((Vec_HIP*)v->spptr)->GPUarray;
2694: #endif
2695:   return(0);
2696: }

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

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

2706:    Input Parameters:
2707: +  v - the vector
2708: -  a - the HIP device pointer.  This pointer is invalid after
2709:        VecHIPRestoreArrayRead() returns.

2711:    Fortran note:
2712:    This function is not currently available from Fortran.

2714:    Level: intermediate

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

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

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

2735:    The device pointer needs to be released with
2736:    VecHIPRestoreArrayWrite().  When the pointer is released the
2737:    host data of the vector is marked as out of data.  Subsequent access
2738:    of the host data with e.g. VecGetArray() incurs a device to host data
2739:    transfer.

2741:    Input Parameter:
2742: .  v - the vector

2744:    Output Parameter:
2745: .  a - the HIP pointer

2747:    Fortran note:
2748:    This function is not currently available from Fortran.

2750:    Level: advanced

2752: .seealso: VecHIPRestoreArrayWrite(), VecHIPGetArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2753: @*/
2754: PETSC_EXTERN PetscErrorCode VecHIPGetArrayWrite(Vec v, PetscScalar **a)
2755: {
2756: #if defined(PETSC_HAVE_HIP)
2758: #endif

2762: #if defined(PETSC_HAVE_HIP)
2763:   *a   = 0;
2764:   VecHIPAllocateCheck(v);
2765:   *a   = ((Vec_HIP*)v->spptr)->GPUarray;
2766: #endif
2767:   return(0);
2768: }

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

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

2777:    Input Parameters:
2778: +  v - the vector
2779: -  a - the HIP device pointer.  This pointer is invalid after
2780:        VecHIPRestoreArrayWrite() returns.

2782:    Fortran note:
2783:    This function is not currently available from Fortran.

2785:    Level: intermediate

2787: .seealso: VecHIPGetArrayWrite(), VecHIPGetArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2788: @*/
2789: PETSC_EXTERN PetscErrorCode VecHIPRestoreArrayWrite(Vec v, PetscScalar **a)
2790: {

2795: #if defined(PETSC_HAVE_HIP)
2796:   v->offloadmask = PETSC_OFFLOAD_GPU;
2797: #endif

2799:   PetscObjectStateIncrease((PetscObject)v);
2800:   return(0);
2801: }

2803: /*@C
2804:    VecHIPPlaceArray - Allows one to replace the GPU array in a vector with a
2805:    GPU array provided by the user. This is useful to avoid copying an
2806:    array into a vector.

2808:    Not Collective

2810:    Input Parameters:
2811: +  vec - the vector
2812: -  array - the GPU array

2814:    Notes:
2815:    You can return to the original GPU array with a call to VecHIPResetArray()
2816:    It is not possible to use VecHIPPlaceArray() and VecPlaceArray() at the
2817:    same time on the same vector.

2819:    Level: developer

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

2823: @*/
2824: PetscErrorCode VecHIPPlaceArray(Vec vin,const PetscScalar a[])
2825: {

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

2841: /*@C
2842:    VecHIPReplaceArray - Allows one to replace the GPU array in a vector
2843:    with a GPU array provided by the user. This is useful to avoid copying
2844:    a GPU array into a vector.

2846:    Not Collective

2848:    Input Parameters:
2849: +  vec - the vector
2850: -  array - the GPU array

2852:    Notes:
2853:    This permanently replaces the GPU array and frees the memory associated
2854:    with the old GPU array.

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

2859:    Not supported from Fortran

2861:    Level: developer

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

2865: @*/
2866: PetscErrorCode VecHIPReplaceArray(Vec vin,const PetscScalar a[])
2867: {
2868: #if defined(PETSC_HAVE_HIP)
2869:   hipError_t err;
2870: #endif

2875: #if defined(PETSC_HAVE_HIP)
2876:   err = hipFree(((Vec_HIP*)vin->spptr)->GPUarray);CHKERRHIP(err);
2877:   ((Vec_HIP*)vin->spptr)->GPUarray = (PetscScalar*)a;
2878:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2879: #endif
2880:   PetscObjectStateIncrease((PetscObject)vin);
2881:   return(0);
2882: }

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

2888:    Not Collective

2890:    Input Parameters:
2891: .  vec - the vector

2893:    Level: developer

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

2897: @*/
2898: PetscErrorCode VecHIPResetArray(Vec vin)
2899: {

2904: #if defined(PETSC_HAVE_HIP)
2905:   VecHIPCopyToGPU(vin);
2906:   ((Vec_HIP*)vin->spptr)->GPUarray = (PetscScalar *) ((Vec_Seq*)vin->data)->unplacedarray;
2907:   ((Vec_Seq*)vin->data)->unplacedarray = 0;
2908:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2909: #endif
2910:   PetscObjectStateIncrease((PetscObject)vin);
2911:   return(0);
2912: }

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

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

2921:     Collective on Vec

2923:     Input Parameters:
2924: +   x - a vector to mimic
2925: -   n - the number of vectors to obtain

2927:     Output Parameters:
2928: +   y - Fortran90 pointer to the array of vectors
2929: -   ierr - error code

2931:     Example of Usage:
2932: .vb
2933: #include <petsc/finclude/petscvec.h>
2934:     use petscvec

2936:     Vec x
2937:     Vec, pointer :: y(:)
2938:     ....
2939:     call VecDuplicateVecsF90(x,2,y,ierr)
2940:     call VecSet(y(2),alpha,ierr)
2941:     call VecSet(y(2),alpha,ierr)
2942:     ....
2943:     call VecDestroyVecsF90(2,y,ierr)
2944: .ve

2946:     Notes:
2947:     Not yet supported for all F90 compilers

2949:     Use VecDestroyVecsF90() to free the space.

2951:     Level: beginner

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

2955: M*/

2957: /*MC
2958:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2959:     VecGetArrayF90().

2961:     Synopsis:
2962:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2964:     Logically Collective on Vec

2966:     Input Parameters:
2967: +   x - vector
2968: -   xx_v - the Fortran90 pointer to the array

2970:     Output Parameter:
2971: .   ierr - error code

2973:     Example of Usage:
2974: .vb
2975: #include <petsc/finclude/petscvec.h>
2976:     use petscvec

2978:     PetscScalar, pointer :: xx_v(:)
2979:     ....
2980:     call VecGetArrayF90(x,xx_v,ierr)
2981:     xx_v(3) = a
2982:     call VecRestoreArrayF90(x,xx_v,ierr)
2983: .ve

2985:     Level: beginner

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

2989: M*/

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

2994:     Synopsis:
2995:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

2997:     Collective on Vec

2999:     Input Parameters:
3000: +   n - the number of vectors previously obtained
3001: -   x - pointer to array of vector pointers

3003:     Output Parameter:
3004: .   ierr - error code

3006:     Notes:
3007:     Not yet supported for all F90 compilers

3009:     Level: beginner

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

3013: M*/

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

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

3024:     Logically Collective on Vec

3026:     Input Parameter:
3027: .   x - vector

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

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

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

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

3047:     Level: beginner

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

3051: M*/

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

3059:     Synopsis:
3060:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

3062:     Logically Collective on Vec

3064:     Input Parameter:
3065: .   x - vector

3067:     Output Parameters:
3068: +   xx_v - the Fortran90 pointer to the array
3069: -   ierr - error code

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

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

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

3085:     Level: beginner

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

3089: M*/

3091: /*MC
3092:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
3093:     VecGetArrayReadF90().

3095:     Synopsis:
3096:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

3098:     Logically Collective on Vec

3100:     Input Parameters:
3101: +   x - vector
3102: -   xx_v - the Fortran90 pointer to the array

3104:     Output Parameter:
3105: .   ierr - error code

3107:     Example of Usage:
3108: .vb
3109: #include <petsc/finclude/petscvec.h>
3110:     use petscvec

3112:     PetscScalar, pointer :: xx_v(:)
3113:     ....
3114:     call VecGetArrayReadF90(x,xx_v,ierr)
3115:     a = xx_v(3)
3116:     call VecRestoreArrayReadF90(x,xx_v,ierr)
3117: .ve

3119:     Level: beginner

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

3123: M*/

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

3130:    Logically Collective

3132:    Input Parameters:
3133: +  x - the vector
3134: .  m - first dimension of two dimensional array
3135: .  n - second dimension of two dimensional array
3136: .  mstart - first index you will use in first coordinate direction (often 0)
3137: -  nstart - first index in the second coordinate direction (often 0)

3139:    Output Parameter:
3140: .  a - location to put pointer to the array

3142:    Level: developer

3144:   Notes:
3145:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
3146:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3147:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3148:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

3166:   VecGetLocalSize(x,&N);
3167:   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);
3168:   VecGetArray(x,&aa);

3170:   PetscMalloc1(m,a);
3171:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
3172:   *a -= mstart;
3173:   return(0);
3174: }

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

3181:    Logically Collective

3183:    Input Parameters:
3184: +  x - the vector
3185: .  m - first dimension of two dimensional array
3186: .  n - second dimension of two dimensional array
3187: .  mstart - first index you will use in first coordinate direction (often 0)
3188: -  nstart - first index in the second coordinate direction (often 0)

3190:    Output Parameter:
3191: .  a - location to put pointer to the array

3193:    Level: developer

3195:   Notes:
3196:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
3197:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3198:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3199:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

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

3219:   VecGetLocalSize(x,&N);
3220:   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);
3221:   VecGetArrayWrite(x,&aa);

3223:   PetscMalloc1(m,a);
3224:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
3225:   *a -= mstart;
3226:   return(0);
3227: }

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

3232:    Logically Collective

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

3242:    Level: developer

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

3250:    This routine actually zeros out the a pointer.

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

3265:   dummy = (void*)(*a + mstart);
3266:   PetscFree(dummy);
3267:   VecRestoreArray(x,NULL);
3268:   return(0);
3269: }

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

3274:    Logically Collective

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

3284:    Level: developer

3286:    Notes:
3287:    For regular PETSc vectors this routine does not involve any copies. For
3288:    any special vectors that do not store local vector data in a contiguous
3289:    array, this routine will copy the data back into the underlying
3290:    vector data structure from the array obtained with VecGetArray().

3292:    This routine actually zeros out the a pointer.

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

3307:   dummy = (void*)(*a + mstart);
3308:   PetscFree(dummy);
3309:   VecRestoreArrayWrite(x,NULL);
3310:   return(0);
3311: }

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

3318:    Logically Collective

3320:    Input Parameters:
3321: +  x - the vector
3322: .  m - first dimension of two dimensional array
3323: -  mstart - first index you will use in first coordinate direction (often 0)

3325:    Output Parameter:
3326: .  a - location to put pointer to the array

3328:    Level: developer

3330:   Notes:
3331:    For a vector obtained from DMCreateLocalVector() mstart are likely
3332:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3333:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3337: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3338:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3339:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3340: @*/
3341: PetscErrorCode  VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3342: {
3344:   PetscInt       N;

3350:   VecGetLocalSize(x,&N);
3351:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3352:   VecGetArray(x,a);
3353:   *a  -= mstart;
3354:   return(0);
3355: }

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

3362:    Logically Collective

3364:    Input Parameters:
3365: +  x - the vector
3366: .  m - first dimension of two dimensional array
3367: -  mstart - first index you will use in first coordinate direction (often 0)

3369:    Output Parameter:
3370: .  a - location to put pointer to the array

3372:    Level: developer

3374:   Notes:
3375:    For a vector obtained from DMCreateLocalVector() mstart are likely
3376:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3377:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3381: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3382:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3383:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3384: @*/
3385: PetscErrorCode  VecGetArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3386: {
3388:   PetscInt       N;

3394:   VecGetLocalSize(x,&N);
3395:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3396:   VecGetArrayWrite(x,a);
3397:   *a  -= mstart;
3398:   return(0);
3399: }

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

3404:    Logically Collective

3406:    Input Parameters:
3407: +  x - the vector
3408: .  m - first dimension of two dimensional array
3409: .  mstart - first index you will use in first coordinate direction (often 0)
3410: -  a - location of pointer to array obtained from VecGetArray21()

3412:    Level: developer

3414:    Notes:
3415:    For regular PETSc vectors this routine does not involve any copies. For
3416:    any special vectors that do not store local vector data in a contiguous
3417:    array, this routine will copy the data back into the underlying
3418:    vector data structure from the array obtained with VecGetArray1d().

3420:    This routine actually zeros out the a pointer.

3422: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3423:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3424:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3425: @*/
3426: PetscErrorCode  VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3427: {

3433:   VecRestoreArray(x,NULL);
3434:   return(0);
3435: }

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

3440:    Logically Collective

3442:    Input Parameters:
3443: +  x - the vector
3444: .  m - first dimension of two dimensional array
3445: .  mstart - first index you will use in first coordinate direction (often 0)
3446: -  a - location of pointer to array obtained from VecGetArray21()

3448:    Level: developer

3450:    Notes:
3451:    For regular PETSc vectors this routine does not involve any copies. For
3452:    any special vectors that do not store local vector data in a contiguous
3453:    array, this routine will copy the data back into the underlying
3454:    vector data structure from the array obtained with VecGetArray1d().

3456:    This routine actually zeros out the a pointer.

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

3460: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3461:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3462:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3463: @*/
3464: PetscErrorCode  VecRestoreArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3465: {

3471:   VecRestoreArrayWrite(x,NULL);
3472:   return(0);
3473: }

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

3480:    Logically Collective

3482:    Input Parameters:
3483: +  x - the vector
3484: .  m - first dimension of three dimensional array
3485: .  n - second dimension of three dimensional array
3486: .  p - third dimension of three dimensional array
3487: .  mstart - first index you will use in first coordinate direction (often 0)
3488: .  nstart - first index in the second coordinate direction (often 0)
3489: -  pstart - first index in the third coordinate direction (often 0)

3491:    Output Parameter:
3492: .  a - location to put pointer to the array

3494:    Level: developer

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

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

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

3518:   VecGetLocalSize(x,&N);
3519:   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);
3520:   VecGetArray(x,&aa);

3522:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3523:   b    = (PetscScalar**)((*a) + m);
3524:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3525:   for (i=0; i<m; i++)
3526:     for (j=0; j<n; j++)
3527:       b[i*n+j] = aa + i*n*p + j*p - pstart;

3529:   *a -= mstart;
3530:   return(0);
3531: }

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

3538:    Logically Collective

3540:    Input Parameters:
3541: +  x - the vector
3542: .  m - first dimension of three dimensional array
3543: .  n - second dimension of three dimensional array
3544: .  p - third dimension of three dimensional array
3545: .  mstart - first index you will use in first coordinate direction (often 0)
3546: .  nstart - first index in the second coordinate direction (often 0)
3547: -  pstart - first index in the third coordinate direction (often 0)

3549:    Output Parameter:
3550: .  a - location to put pointer to the array

3552:    Level: developer

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

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

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

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

3578:   VecGetLocalSize(x,&N);
3579:   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);
3580:   VecGetArrayWrite(x,&aa);

3582:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3583:   b    = (PetscScalar**)((*a) + m);
3584:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3585:   for (i=0; i<m; i++)
3586:     for (j=0; j<n; j++)
3587:       b[i*n+j] = aa + i*n*p + j*p - pstart;

3589:   *a -= mstart;
3590:   return(0);
3591: }

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

3596:    Logically Collective

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

3608:    Level: developer

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

3616:    This routine actually zeros out the a pointer.

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

3631:   dummy = (void*)(*a + mstart);
3632:   PetscFree(dummy);
3633:   VecRestoreArray(x,NULL);
3634:   return(0);
3635: }

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

3640:    Logically Collective

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

3652:    Level: developer

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

3660:    This routine actually zeros out the a pointer.

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

3675:   dummy = (void*)(*a + mstart);
3676:   PetscFree(dummy);
3677:   VecRestoreArrayWrite(x,NULL);
3678:   return(0);
3679: }

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

3686:    Logically Collective

3688:    Input Parameters:
3689: +  x - the vector
3690: .  m - first dimension of four dimensional array
3691: .  n - second dimension of four dimensional array
3692: .  p - third dimension of four dimensional array
3693: .  q - fourth dimension of four dimensional array
3694: .  mstart - first index you will use in first coordinate direction (often 0)
3695: .  nstart - first index in the second coordinate direction (often 0)
3696: .  pstart - first index in the third coordinate direction (often 0)
3697: -  qstart - first index in the fourth coordinate direction (often 0)

3699:    Output Parameter:
3700: .  a - location to put pointer to the array

3702:    Level: beginner

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

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

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

3726:   VecGetLocalSize(x,&N);
3727:   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);
3728:   VecGetArray(x,&aa);

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

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

3750:    Logically Collective

3752:    Input Parameters:
3753: +  x - the vector
3754: .  m - first dimension of four dimensional array
3755: .  n - second dimension of four dimensional array
3756: .  p - third dimension of four dimensional array
3757: .  q - fourth dimension of four dimensional array
3758: .  mstart - first index you will use in first coordinate direction (often 0)
3759: .  nstart - first index in the second coordinate direction (often 0)
3760: .  pstart - first index in the third coordinate direction (often 0)
3761: -  qstart - first index in the fourth coordinate direction (often 0)

3763:    Output Parameter:
3764: .  a - location to put pointer to the array

3766:    Level: beginner

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

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

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

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

3792:   VecGetLocalSize(x,&N);
3793:   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);
3794:   VecGetArrayWrite(x,&aa);

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

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

3814:    Logically Collective

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

3828:    Level: beginner

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

3836:    This routine actually zeros out the a pointer.

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

3851:   dummy = (void*)(*a + mstart);
3852:   PetscFree(dummy);
3853:   VecRestoreArray(x,NULL);
3854:   return(0);
3855: }

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

3860:    Logically Collective

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

3874:    Level: beginner

3876:    Notes:
3877:    For regular PETSc vectors this routine does not involve any copies. For
3878:    any special vectors that do not store local vector data in a contiguous
3879:    array, this routine will copy the data back into the underlying
3880:    vector data structure from the array obtained with VecGetArray().

3882:    This routine actually zeros out the a pointer.

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

3897:   dummy = (void*)(*a + mstart);
3898:   PetscFree(dummy);
3899:   VecRestoreArrayWrite(x,NULL);
3900:   return(0);
3901: }

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

3908:    Logically Collective

3910:    Input Parameters:
3911: +  x - the vector
3912: .  m - first dimension of two dimensional array
3913: .  n - second dimension of two dimensional array
3914: .  mstart - first index you will use in first coordinate direction (often 0)
3915: -  nstart - first index in the second coordinate direction (often 0)

3917:    Output Parameter:
3918: .  a - location to put pointer to the array

3920:    Level: developer

3922:   Notes:
3923:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
3924:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3925:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3926:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

3930: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3931:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3932:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3933: @*/
3934: PetscErrorCode  VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3935: {
3936:   PetscErrorCode    ierr;
3937:   PetscInt          i,N;
3938:   const PetscScalar *aa;

3944:   VecGetLocalSize(x,&N);
3945:   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);
3946:   VecGetArrayRead(x,&aa);

3948:   PetscMalloc1(m,a);
3949:   for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
3950:   *a -= mstart;
3951:   return(0);
3952: }

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

3957:    Logically Collective

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

3967:    Level: developer

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

3975:    This routine actually zeros out the a pointer.

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

3990:   dummy = (void*)(*a + mstart);
3991:   PetscFree(dummy);
3992:   VecRestoreArrayRead(x,NULL);
3993:   return(0);
3994: }

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

4001:    Logically Collective

4003:    Input Parameters:
4004: +  x - the vector
4005: .  m - first dimension of two dimensional array
4006: -  mstart - first index you will use in first coordinate direction (often 0)

4008:    Output Parameter:
4009: .  a - location to put pointer to the array

4011:    Level: developer

4013:   Notes:
4014:    For a vector obtained from DMCreateLocalVector() mstart are likely
4015:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
4016:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

4020: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
4021:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
4022:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
4023: @*/
4024: PetscErrorCode  VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
4025: {
4027:   PetscInt       N;

4033:   VecGetLocalSize(x,&N);
4034:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
4035:   VecGetArrayRead(x,(const PetscScalar**)a);
4036:   *a  -= mstart;
4037:   return(0);
4038: }

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

4043:    Logically Collective

4045:    Input Parameters:
4046: +  x - the vector
4047: .  m - first dimension of two dimensional array
4048: .  mstart - first index you will use in first coordinate direction (often 0)
4049: -  a - location of pointer to array obtained from VecGetArray21()

4051:    Level: developer

4053:    Notes:
4054:    For regular PETSc vectors this routine does not involve any copies. For
4055:    any special vectors that do not store local vector data in a contiguous
4056:    array, this routine will copy the data back into the underlying
4057:    vector data structure from the array obtained with VecGetArray1dRead().

4059:    This routine actually zeros out the a pointer.

4061: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
4062:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
4063:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
4064: @*/
4065: PetscErrorCode  VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
4066: {

4072:   VecRestoreArrayRead(x,NULL);
4073:   return(0);
4074: }

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

4081:    Logically Collective

4083:    Input Parameters:
4084: +  x - the vector
4085: .  m - first dimension of three dimensional array
4086: .  n - second dimension of three dimensional array
4087: .  p - third dimension of three dimensional array
4088: .  mstart - first index you will use in first coordinate direction (often 0)
4089: .  nstart - first index in the second coordinate direction (often 0)
4090: -  pstart - first index in the third coordinate direction (often 0)

4092:    Output Parameter:
4093: .  a - location to put pointer to the array

4095:    Level: developer

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

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

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

4120:   VecGetLocalSize(x,&N);
4121:   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);
4122:   VecGetArrayRead(x,&aa);

4124:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
4125:   b    = (PetscScalar**)((*a) + m);
4126:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
4127:   for (i=0; i<m; i++)
4128:     for (j=0; j<n; j++)
4129:       b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;

4131:   *a -= mstart;
4132:   return(0);
4133: }

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

4138:    Logically Collective

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

4150:    Level: developer

4152:    Notes:
4153:    For regular PETSc vectors this routine does not involve any copies. For
4154:    any special vectors that do not store local vector data in a contiguous
4155:    array, this routine will copy the data back into the underlying
4156:    vector data structure from the array obtained with VecGetArray().

4158:    This routine actually zeros out the a pointer.

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

4173:   dummy = (void*)(*a + mstart);
4174:   PetscFree(dummy);
4175:   VecRestoreArrayRead(x,NULL);
4176:   return(0);
4177: }

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

4184:    Logically Collective

4186:    Input Parameters:
4187: +  x - the vector
4188: .  m - first dimension of four dimensional array
4189: .  n - second dimension of four dimensional array
4190: .  p - third dimension of four dimensional array
4191: .  q - fourth dimension of four dimensional array
4192: .  mstart - first index you will use in first coordinate direction (often 0)
4193: .  nstart - first index in the second coordinate direction (often 0)
4194: .  pstart - first index in the third coordinate direction (often 0)
4195: -  qstart - first index in the fourth coordinate direction (often 0)

4197:    Output Parameter:
4198: .  a - location to put pointer to the array

4200:    Level: beginner

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

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

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

4225:   VecGetLocalSize(x,&N);
4226:   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);
4227:   VecGetArrayRead(x,&aa);

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

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

4247:    Logically Collective

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

4261:    Level: beginner

4263:    Notes:
4264:    For regular PETSc vectors this routine does not involve any copies. For
4265:    any special vectors that do not store local vector data in a contiguous
4266:    array, this routine will copy the data back into the underlying
4267:    vector data structure from the array obtained with VecGetArray().

4269:    This routine actually zeros out the a pointer.

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

4284:   dummy = (void*)(*a + mstart);
4285:   PetscFree(dummy);
4286:   VecRestoreArrayRead(x,NULL);
4287:   return(0);
4288: }

4290: #if defined(PETSC_USE_DEBUG)

4292: /*@
4293:    VecLockGet  - Gets the current lock status of a vector

4295:    Logically Collective on Vec

4297:    Input Parameter:
4298: .  x - the vector

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

4304:    Level: beginner

4306: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop()
4307: @*/
4308: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
4309: {
4312:   *state = x->lock;
4313:   return(0);
4314: }

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

4319:    Logically Collective on Vec

4321:    Input Parameter:
4322: .  x - the vector

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

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

4330:    Level: beginner

4332: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPop(), VecLockGet()
4333: @*/
4334: PetscErrorCode VecLockReadPush(Vec x)
4335: {
4338:   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");
4339:   x->lock++;
4340:   return(0);
4341: }

4343: /*@
4344:    VecLockReadPop  - Pops a read-only lock from a vector

4346:    Logically Collective on Vec

4348:    Input Parameter:
4349: .  x - the vector

4351:    Level: beginner

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

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

4367:    Logically Collective on Vec

4369:    Input Parameters:
4370: +  x   - the vector
4371: -  flg - PETSC_TRUE to lock the vector for writing; PETSC_FALSE to unlock it.

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

4379:        VecGetArray(x,&xdata); // begin phase
4380:        VecLockWriteSet_Private(v,PETSC_TRUE);

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

4384:        VecRestoreArray(x,&vdata); // end phase
4385:        VecLockWriteSet_Private(v,PETSC_FALSE);

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

4390:    Level: beginner

4392: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop(), VecLockGet()
4393: @*/
4394: PetscErrorCode VecLockWriteSet_Private(Vec x,PetscBool flg)
4395: {
4398:   if (flg) {
4399:     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");
4400:     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");
4401:     else x->lock = -1;
4402:   } else {
4403:     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");
4404:     x->lock = 0;
4405:   }
4406:   return(0);
4407: }

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

4412:    Level: deprecated

4414: .seealso: VecLockReadPush()
4415: @*/
4416: PetscErrorCode VecLockPush(Vec x)
4417: {
4420:   VecLockReadPush(x);
4421:   return(0);
4422: }

4424: /*@
4425:    VecLockPop  - Pops a read-only lock from a vector

4427:    Level: deprecated

4429: .seealso: VecLockReadPop()
4430: @*/
4431: PetscErrorCode VecLockPop(Vec x)
4432: {
4435:   VecLockReadPop(x);
4436:   return(0);
4437: }

4439: #endif