Actual source code: rvector.c

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

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

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

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

 45:    Logically Collective on Vec

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

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

 53:    Level: advanced

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

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

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

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

 80:    Collective on Vec

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

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

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

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

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

104:    Level: intermediate


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

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

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

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

131:    Collective on Vec

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

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

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

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

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

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

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

154:    Level: intermediate


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

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

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

173:    Collective on Vec

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

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

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

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

194:    Level: intermediate

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


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

205: @*/

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


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

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

236:    Not Collective

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

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

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

253:    Level: intermediate

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

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


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

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


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

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

289:    Collective on Vec

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

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

298:    Level: intermediate


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

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

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

326:    Collective on Vec

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

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

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

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


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

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

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

361:    Collective on Vec

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

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

370:    Level: intermediate

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

375:    This returns the smallest index with the minumum value


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

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

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

398:    Collective on Vec

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

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

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

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

415:    Level: intermediate

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

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

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

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

441:    Not collective on Vec

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

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

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

454:    Level: intermediate


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

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

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

492:    Logically Collective on Vec

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

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

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

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

511:    Level: beginner

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

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

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

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

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


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

556:    Logically Collective on Vec

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

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

565:    Level: intermediate

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

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


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

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

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

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

609:    Logically Collective on Vec

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

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

618:    Level: intermediate

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


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

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

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

653:    Logically Collective on Vec

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

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

662:    Level: intermediate

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


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

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

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

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

704:    Logically Collective on Vec

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

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

713:    Level: intermediate

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


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

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

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


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

748:    Logically Collective on Vec

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

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

757:    Level: intermediate

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


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

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

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


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

797:    Not Collective

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

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

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

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

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

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

826:    Level: beginner

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

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

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

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

853:     Not Collective

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

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

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

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

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

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

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

876:    Level: beginner

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

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

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

897:    Not Collective

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

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

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

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

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

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

926:    Level: intermediate

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

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

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


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

954:    Not Collective

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

965:    Level: intermediate

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

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

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

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

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

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

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

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

1017:    Not Collective

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

1028:    Level: intermediate

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

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

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

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


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

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

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

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

1077:    Collective on Vec

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

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

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

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

1096:    Level: intermediate


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

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

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

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

1126:    Collective on Vec

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

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

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

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

1145:    Level: intermediate


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

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

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

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

1176:    Logically Collective on Vec

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

1184:    Level: intermediate

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

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

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

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

1223:    Collective on IS

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

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

1232:    Level: advanced

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

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

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

1251:   if (X->ops->getsubvector) {
1252:     (*X->ops->getsubvector)(X,is,&Z);
1253:   } else { /* Default implementation currently does no caching */
1254:     PetscInt  gstart,gend,start;
1255:     PetscBool contiguous,gcontiguous;
1256:     VecGetOwnershipRange(X,&gstart,&gend);
1257:     ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1258:     MPIU_Allreduce(&contiguous,&gcontiguous,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1259:     if (gcontiguous) { /* We can do a no-copy implementation */
1260:       PetscInt n,N,bs;
1261:       PetscInt state;

1263:       ISGetSize(is,&N);
1264:       ISGetLocalSize(is,&n);
1265:       VecGetBlockSize(X,&bs);
1266:       if (n%bs || bs == 1 || !n) bs = -1; /* Do not decide block size if we do not have to */
1267:       VecLockGet(X,&state);
1268:       if (state) {
1269:         const PetscScalar *x;
1270:         VecGetArrayRead(X,&x);
1271:         VecCreate(PetscObjectComm((PetscObject)X),&Z);
1272:         VecSetType(Z,((PetscObject)X)->type_name);
1273:         VecSetSizes(Z,n,N);
1274:         VecSetBlockSize(Z,bs);
1275:         VecPlaceArray(Z,(PetscScalar*)x+start);
1276:         VecLockReadPush(Z);
1277:         VecRestoreArrayRead(X,&x);
1278:       } else {
1279:         PetscScalar *x;
1280:         VecGetArray(X,&x);
1281:         VecCreate(PetscObjectComm((PetscObject)X),&Z);
1282:         VecSetType(Z,((PetscObject)X)->type_name);
1283:         VecSetSizes(Z,n,N);
1284:         VecSetBlockSize(Z,bs);
1285:         VecPlaceArray(Z,(PetscScalar*)x+start);
1286:         VecRestoreArray(X,&x);
1287:       }
1288:       Z->ops->placearray = NULL;
1289:       Z->ops->replacearray = NULL;
1290:     } else { /* Have to create a scatter and do a copy */
1291:       VecScatter scatter;
1292:       PetscInt   n,N;
1293:       ISGetLocalSize(is,&n);
1294:       ISGetSize(is,&N);
1295:       VecCreate(PetscObjectComm((PetscObject)is),&Z);
1296:       VecSetSizes(Z,n,N);
1297:       VecSetType(Z,((PetscObject)X)->type_name);
1298:       VecScatterCreate(X,is,Z,NULL,&scatter);
1299:       VecScatterBegin(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1300:       VecScatterEnd(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1301:       PetscObjectCompose((PetscObject)Z,"VecGetSubVector_Scatter",(PetscObject)scatter);
1302:       VecScatterDestroy(&scatter);
1303:     }
1304:   }
1305:   /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1306:   if (VecGetSubVectorSavedStateId < 0) {PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);}
1307:   PetscObjectComposedDataSetInt((PetscObject)Z,VecGetSubVectorSavedStateId,1);
1308:   *Y   = Z;
1309:   return(0);
1310: }

1312: /*@
1313:    VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()

1315:    Collective on IS

1317:    Input Arguments:
1318: + X - vector from which subvector was obtained
1319: . is - index set representing the subset of X
1320: - Y - subvector being restored

1322:    Level: advanced

1324: .seealso: VecGetSubVector()
1325: @*/
1326: PetscErrorCode  VecRestoreSubVector(Vec X,IS is,Vec *Y)
1327: {

1335:   if (X->ops->restoresubvector) {
1336:     (*X->ops->restoresubvector)(X,is,Y);
1337:   } else {
1338:     PETSC_UNUSED PetscObjectState dummystate = 0;
1339:     PetscBool valid;
1340:     PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,dummystate,valid);
1341:     if (!valid) {
1342:       VecScatter scatter;

1344:       PetscObjectQuery((PetscObject)*Y,"VecGetSubVector_Scatter",(PetscObject*)&scatter);
1345:       if (scatter) {
1346:         VecScatterBegin(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1347:         VecScatterEnd(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1348:       }
1349:     }
1350:     VecDestroy(Y);
1351:   }
1352:   return(0);
1353: }

1355: /*@
1356:    VecGetLocalVectorRead - Maps the local portion of a vector into a
1357:    vector.  You must call VecRestoreLocalVectorRead() when the local
1358:    vector is no longer needed.

1360:    Not collective.

1362:    Input parameter:
1363: .  v - The vector for which the local vector is desired.

1365:    Output parameter:
1366: .  w - Upon exit this contains the local vector.

1368:    Level: beginner

1370:    Notes:
1371:    This function is similar to VecGetArrayRead() which maps the local
1372:    portion into a raw pointer.  VecGetLocalVectorRead() is usually
1373:    almost as efficient as VecGetArrayRead() but in certain circumstances
1374:    VecGetLocalVectorRead() can be much more efficient than
1375:    VecGetArrayRead().  This is because the construction of a contiguous
1376:    array representing the vector data required by VecGetArrayRead() can
1377:    be an expensive operation for certain vector types.  For example, for
1378:    GPU vectors VecGetArrayRead() requires that the data between device
1379:    and host is synchronized.

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

1384: .seealso: VecRestoreLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1385: @*/
1386: PetscErrorCode VecGetLocalVectorRead(Vec v,Vec w)
1387: {
1389:   PetscScalar    *a;

1394:   VecCheckSameLocalSize(v,1,w,2);
1395:   if (v->ops->getlocalvectorread) {
1396:     (*v->ops->getlocalvectorread)(v,w);
1397:   } else {
1398:     VecGetArrayRead(v,(const PetscScalar**)&a);
1399:     VecPlaceArray(w,a);
1400:   }
1401:   return(0);
1402: }

1404: /*@
1405:    VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1406:    previously mapped into a vector using VecGetLocalVectorRead().

1408:    Not collective.

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

1414:    Level: beginner

1416: .seealso: VecGetLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1417: @*/
1418: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec w)
1419: {
1421:   PetscScalar    *a;

1426:   if (v->ops->restorelocalvectorread) {
1427:     (*v->ops->restorelocalvectorread)(v,w);
1428:   } else {
1429:     VecGetArrayRead(w,(const PetscScalar**)&a);
1430:     VecRestoreArrayRead(v,(const PetscScalar**)&a);
1431:     VecResetArray(w);
1432:   }
1433:   return(0);
1434: }

1436: /*@
1437:    VecGetLocalVector - Maps the local portion of a vector into a
1438:    vector.

1440:    Collective on v, not collective on w.

1442:    Input parameter:
1443: .  v - The vector for which the local vector is desired.

1445:    Output parameter:
1446: .  w - Upon exit this contains the local vector.

1448:    Level: beginner

1450:    Notes:
1451:    This function is similar to VecGetArray() which maps the local
1452:    portion into a raw pointer.  VecGetLocalVector() is usually about as
1453:    efficient as VecGetArray() but in certain circumstances
1454:    VecGetLocalVector() can be much more efficient than VecGetArray().
1455:    This is because the construction of a contiguous array representing
1456:    the vector data required by VecGetArray() can be an expensive
1457:    operation for certain vector types.  For example, for GPU vectors
1458:    VecGetArray() requires that the data between device and host is
1459:    synchronized.

1461: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1462: @*/
1463: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1464: {
1466:   PetscScalar    *a;

1471:   VecCheckSameLocalSize(v,1,w,2);
1472:   if (v->ops->getlocalvector) {
1473:     (*v->ops->getlocalvector)(v,w);
1474:   } else {
1475:     VecGetArray(v,&a);
1476:     VecPlaceArray(w,a);
1477:   }
1478:   return(0);
1479: }

1481: /*@
1482:    VecRestoreLocalVector - Unmaps the local portion of a vector
1483:    previously mapped into a vector using VecGetLocalVector().

1485:    Logically collective.

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

1491:    Level: beginner

1493: .seealso: VecGetLocalVector(), VecGetLocalVectorRead(), VecRestoreLocalVectorRead(), LocalVectorRead(), VecGetArrayRead(), VecGetArray()
1494: @*/
1495: PetscErrorCode VecRestoreLocalVector(Vec v,Vec w)
1496: {
1498:   PetscScalar    *a;

1503:   if (v->ops->restorelocalvector) {
1504:     (*v->ops->restorelocalvector)(v,w);
1505:   } else {
1506:     VecGetArray(w,&a);
1507:     VecRestoreArray(v,&a);
1508:     VecResetArray(w);
1509:   }
1510:   return(0);
1511: }

1513: /*@C
1514:    VecGetArray - Returns a pointer to a contiguous array that contains this
1515:    processor's portion of the vector data. For the standard PETSc
1516:    vectors, VecGetArray() returns a pointer to the local data array and
1517:    does not use any copies. If the underlying vector data is not stored
1518:    in a contiguous array this routine will copy the data to a contiguous
1519:    array and return a pointer to that. You MUST call VecRestoreArray()
1520:    when you no longer need access to the array.

1522:    Logically Collective on Vec

1524:    Input Parameter:
1525: .  x - the vector

1527:    Output Parameter:
1528: .  a - location to put pointer to the array

1530:    Fortran Note:
1531:    This routine is used differently from Fortran 77
1532: $    Vec         x
1533: $    PetscScalar x_array(1)
1534: $    PetscOffset i_x
1535: $    PetscErrorCode ierr
1536: $       call VecGetArray(x,x_array,i_x,ierr)
1537: $
1538: $   Access first local entry in vector with
1539: $      value = x_array(i_x + 1)
1540: $
1541: $      ...... other code
1542: $       call VecRestoreArray(x,x_array,i_x,ierr)
1543:    For Fortran 90 see VecGetArrayF90()

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

1548:    Level: beginner

1550: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1551:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1552: @*/
1553: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1554: {
1556: #if defined(PETSC_HAVE_VIENNACL)
1557:   PetscBool      is_viennacltype = PETSC_FALSE;
1558: #endif

1562:   VecSetErrorIfLocked(x,1);
1563:   if (x->petscnative) {
1564: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1565:     if (x->offloadmask == PETSC_OFFLOAD_GPU) {
1566: #if defined(PETSC_HAVE_VIENNACL)
1567:       PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1568:       if (is_viennacltype) {
1569:         VecViennaCLCopyFromGPU(x);
1570:       } else
1571: #endif
1572:       {
1573: #if defined(PETSC_HAVE_CUDA)
1574:         VecCUDACopyFromGPU(x);
1575: #endif
1576:       }
1577:     } else if (x->offloadmask == PETSC_OFFLOAD_UNALLOCATED) {
1578: #if defined(PETSC_HAVE_VIENNACL)
1579:       PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1580:       if (is_viennacltype) {
1581:         VecViennaCLAllocateCheckHost(x);
1582:       } else
1583: #endif
1584:       {
1585: #if defined(PETSC_HAVE_CUDA)
1586:         VecCUDAAllocateCheckHost(x);
1587: #endif
1588:       }
1589:     }
1590: #endif
1591:     *a = *((PetscScalar**)x->data);
1592:   } else {
1593:     if (x->ops->getarray) {
1594:       (*x->ops->getarray)(x,a);
1595:     } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array for vector type \"%s\"",((PetscObject)x)->type_name);
1596:   }
1597:   return(0);
1598: }

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

1605:    Logically Collective on Vec

1607:    Input Parameter:
1608: .  x - the vector

1610:    Output Parameter:
1611: .  a - location to put pointer to the array

1613:    Level: beginner

1615: .seealso: VecRestoreArrayInPlace(), VecRestoreArrayInPlace(), VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(),
1616:           VecPlaceArray(), VecGetArray2d(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1617: @*/
1618: PetscErrorCode VecGetArrayInPlace(Vec x,PetscScalar **a)
1619: {

1624:   VecSetErrorIfLocked(x,1);

1626: #if defined(PETSC_HAVE_CUDA)
1627:   if (x->petscnative && (x->offloadmask & PETSC_OFFLOAD_GPU)) { /* Prefer working on GPU when offloadmask is PETSC_OFFLOAD_BOTH */
1628:     PetscBool is_cudatype = PETSC_FALSE;
1629:     PetscObjectTypeCompareAny((PetscObject)x,&is_cudatype,VECSEQCUDA,VECMPICUDA,VECCUDA,"");
1630:     if (is_cudatype) {
1631:       VecCUDAGetArray(x,a);
1632:       x->offloadmask = PETSC_OFFLOAD_GPU; /* Change the mask once GPU gets write access, don't wait until restore array */
1633:       return(0);
1634:     }
1635:   }
1636: #endif
1637:   VecGetArray(x,a);
1638:   return(0);
1639: }

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

1646:    Logically Collective on Vec

1648:    Input Parameter:
1649: .  x - the vector

1651:    Output Parameter:
1652: .  a - location to put pointer to the array

1654:    Level: intermediate

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

1659:    Concepts: vector^accessing local values

1661: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1662:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArray(), VecRestoreArrayWrite()
1663: @*/
1664: PetscErrorCode VecGetArrayWrite(Vec x,PetscScalar **a)
1665: {

1670:   VecSetErrorIfLocked(x,1);
1671:   if (!x->ops->getarraywrite) {
1672:     VecGetArray(x,a);
1673:   } else {
1674:     (*x->ops->getarraywrite)(x,a);
1675:   }
1676:   return(0);
1677: }

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

1682:    Not Collective

1684:    Input Parameters:
1685: .  x - the vector

1687:    Output Parameter:
1688: .  a - the array

1690:    Level: beginner

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

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

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

1701: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1702: @*/
1703: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1704: {
1706: #if defined(PETSC_HAVE_VIENNACL)
1707:   PetscBool      is_viennacltype = PETSC_FALSE;
1708: #endif

1712:   if (x->petscnative) {
1713: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1714:     if (x->offloadmask == PETSC_OFFLOAD_GPU) {
1715: #if defined(PETSC_HAVE_VIENNACL)
1716:       PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1717:       if (is_viennacltype) {
1718:         VecViennaCLCopyFromGPU(x);
1719:       } else
1720: #endif
1721:       {
1722: #if defined(PETSC_HAVE_CUDA)
1723:         VecCUDACopyFromGPU(x);
1724: #endif
1725:       }
1726:     }
1727: #endif
1728:     *a = *((PetscScalar **)x->data);
1729:   } else if (x->ops->getarrayread) {
1730:     (*x->ops->getarrayread)(x,a);
1731:   } else {
1732:     (*x->ops->getarray)(x,(PetscScalar**)a);
1733:   }
1734:   return(0);
1735: }

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

1742:    Not Collective

1744:    Input Parameters:
1745: .  x - the vector

1747:    Output Parameter:
1748: .  a - the array

1750:    Level: beginner

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


1756: .seealso: VecRestoreArrayReadInPlace(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayInPlace()
1757: @*/
1758: PetscErrorCode VecGetArrayReadInPlace(Vec x,const PetscScalar **a)
1759: {

1764: #if defined(PETSC_HAVE_CUDA)
1765:   if (x->petscnative && x->offloadmask & PETSC_OFFLOAD_GPU) {
1766:     PetscBool is_cudatype = PETSC_FALSE;
1767:     PetscObjectTypeCompareAny((PetscObject)x,&is_cudatype,VECSEQCUDA,VECMPICUDA,VECCUDA,"");
1768:     if (is_cudatype) {
1769:       VecCUDAGetArrayRead(x,a);
1770:       return(0);
1771:     }
1772:   }
1773: #endif
1774:   VecGetArrayRead(x,a);
1775:   return(0);
1776: }

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

1783:    Logically Collective on Vec

1785:    Input Parameter:
1786: +  x - the vectors
1787: -  n - the number of vectors

1789:    Output Parameter:
1790: .  a - location to put pointer to the array

1792:    Fortran Note:
1793:    This routine is not supported in Fortran.

1795:    Level: intermediate

1797: .seealso: VecGetArray(), VecRestoreArrays()
1798: @*/
1799: PetscErrorCode  VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1800: {
1802:   PetscInt       i;
1803:   PetscScalar    **q;

1809:   if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1810:   PetscMalloc1(n,&q);
1811:   for (i=0; i<n; ++i) {
1812:     VecGetArray(x[i],&q[i]);
1813:   }
1814:   *a = q;
1815:   return(0);
1816: }

1818: /*@C
1819:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1820:    has been called.

1822:    Logically Collective on Vec

1824:    Input Parameters:
1825: +  x - the vector
1826: .  n - the number of vectors
1827: -  a - location of pointer to arrays obtained from VecGetArrays()

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

1835:    Fortran Note:
1836:    This routine is not supported in Fortran.

1838:    Level: intermediate

1840: .seealso: VecGetArrays(), VecRestoreArray()
1841: @*/
1842: PetscErrorCode  VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1843: {
1845:   PetscInt       i;
1846:   PetscScalar    **q = *a;


1853:   for (i=0; i<n; ++i) {
1854:     VecRestoreArray(x[i],&q[i]);
1855:   }
1856:   PetscFree(q);
1857:   return(0);
1858: }

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

1863:    Logically Collective on Vec

1865:    Input Parameters:
1866: +  x - the vector
1867: -  a - location of pointer to array obtained from VecGetArray()

1869:    Level: beginner

1871: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1872:           VecGetArrayPair(), VecRestoreArrayPair()
1873: @*/
1874: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1875: {

1880:   if (x->petscnative) {
1881: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1882:     x->offloadmask = PETSC_OFFLOAD_CPU;
1883: #endif
1884:   } else {
1885:     (*x->ops->restorearray)(x,a);
1886:   }
1887:   if (a) *a = NULL;
1888:   PetscObjectStateIncrease((PetscObject)x);
1889:   return(0);
1890: }

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

1895:    Logically Collective on Vec

1897:    Input Parameters:
1898: +  x - the vector
1899: -  a - location of pointer to array obtained from VecGetArrayInPlace()

1901:    Level: beginner

1903: .seealso: VecGetArrayInPlace(), VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(),
1904:           VecPlaceArray(), VecRestoreArray2d(), VecGetArrayPair(), VecRestoreArrayPair()
1905: @*/
1906: PetscErrorCode VecRestoreArrayInPlace(Vec x,PetscScalar **a)
1907: {

1912: #if defined(PETSC_HAVE_CUDA)
1913:   if (x->petscnative && x->offloadmask == PETSC_OFFLOAD_GPU) {
1914:     PetscBool is_cudatype = PETSC_FALSE;
1915:     PetscObjectTypeCompareAny((PetscObject)x,&is_cudatype,VECSEQCUDA,VECMPICUDA,VECCUDA,"");
1916:     if (is_cudatype) {
1917:       VecCUDARestoreArray(x,a);
1918:       return(0);
1919:     }
1920:   }
1921: #endif
1922:   VecRestoreArray(x,a);
1923:   return(0);
1924: }


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

1930:    Logically Collective on Vec

1932:    Input Parameters:
1933: +  x - the vector
1934: -  a - location of pointer to array obtained from VecGetArray()

1936:    Level: beginner

1938: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1939:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite()
1940: @*/
1941: PetscErrorCode VecRestoreArrayWrite(Vec x,PetscScalar **a)
1942: {

1947:   if (x->petscnative) {
1948: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1949:     x->offloadmask = PETSC_OFFLOAD_CPU;
1950: #endif
1951:   } else {
1952:     if (x->ops->restorearraywrite) {
1953:       (*x->ops->restorearraywrite)(x,a);
1954:     } else {
1955:       (*x->ops->restorearray)(x,a);
1956:     }
1957:   }
1958:   if (a) *a = NULL;
1959:   PetscObjectStateIncrease((PetscObject)x);
1960:   return(0);
1961: }

1963: /*@C
1964:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

1966:    Not Collective

1968:    Input Parameters:
1969: +  vec - the vector
1970: -  array - the array

1972:    Level: beginner

1974: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1975: @*/
1976: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1977: {

1982:   if (x->petscnative) {
1983:     /* nothing */
1984:   } else if (x->ops->restorearrayread) {
1985:     (*x->ops->restorearrayread)(x,a);
1986:   } else {
1987:     (*x->ops->restorearray)(x,(PetscScalar**)a);
1988:   }
1989:   if (a) *a = NULL;
1990:   return(0);
1991: }

1993: /*@C
1994:    VecRestoreArrayReadInPlace - Restore array obtained with VecGetArrayReadInPlace()

1996:    Not Collective

1998:    Input Parameters:
1999: +  vec - the vector
2000: -  array - the array

2002:    Level: beginner

2004: .seealso: VecGetArrayReadInPlace(), VecRestoreArrayInPlace(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
2005: @*/
2006: PetscErrorCode VecRestoreArrayReadInPlace(Vec x,const PetscScalar **a)
2007: {

2011:   VecRestoreArrayRead(x,a);
2012:   return(0);
2013: }

2015: /*@
2016:    VecPlaceArray - Allows one to replace the array in a vector with an
2017:    array provided by the user. This is useful to avoid copying an array
2018:    into a vector.

2020:    Not Collective

2022:    Input Parameters:
2023: +  vec - the vector
2024: -  array - the array

2026:    Notes:
2027:    You can return to the original array with a call to VecResetArray()

2029:    Level: developer

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

2033: @*/
2034: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
2035: {

2042:   if (vec->ops->placearray) {
2043:     (*vec->ops->placearray)(vec,array);
2044:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
2045:   PetscObjectStateIncrease((PetscObject)vec);
2046:   return(0);
2047: }

2049: /*@C
2050:    VecReplaceArray - Allows one to replace the array in a vector with an
2051:    array provided by the user. This is useful to avoid copying an array
2052:    into a vector.

2054:    Not Collective

2056:    Input Parameters:
2057: +  vec - the vector
2058: -  array - the array

2060:    Notes:
2061:    This permanently replaces the array and frees the memory associated
2062:    with the old array.

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

2067:    Not supported from Fortran

2069:    Level: developer

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

2073: @*/
2074: PetscErrorCode  VecReplaceArray(Vec vec,const PetscScalar array[])
2075: {

2081:   if (vec->ops->replacearray) {
2082:     (*vec->ops->replacearray)(vec,array);
2083:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
2084:   PetscObjectStateIncrease((PetscObject)vec);
2085:   return(0);
2086: }


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

2092:    This function has semantics similar to VecGetArray():  the pointer
2093:    returned by this function points to a consistent view of the vector
2094:    data.  This may involve a copy operation of data from the host to the
2095:    device if the data on the device is out of date.  If the device
2096:    memory hasn't been allocated previously it will be allocated as part
2097:    of this function call.  VecCUDAGetArray() assumes that
2098:    the user will modify the vector data.  This is similar to
2099:    intent(inout) in fortran.

2101:    The CUDA device pointer has to be released by calling
2102:    VecCUDARestoreArray().  Upon restoring the vector data
2103:    the data on the host will be marked as out of date.  A subsequent
2104:    access of the host data will thus incur a data transfer from the
2105:    device to the host.


2108:    Input Parameter:
2109: .  v - the vector

2111:    Output Parameter:
2112: .  a - the CUDA device pointer

2114:    Fortran note:
2115:    This function is not currently available from Fortran.

2117:    Level: intermediate

2119: .seealso: VecCUDARestoreArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2120: @*/
2121: PETSC_EXTERN PetscErrorCode VecCUDAGetArray(Vec v, PetscScalar **a)
2122: {
2123: #if defined(PETSC_HAVE_CUDA)
2125: #endif

2129: #if defined(PETSC_HAVE_CUDA)
2130:   *a   = 0;
2131:   VecCUDACopyToGPU(v);
2132:   *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2133: #endif
2134:   return(0);
2135: }

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

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

2144:    Input Parameter:
2145: +  v - the vector
2146: -  a - the CUDA device pointer.  This pointer is invalid after
2147:        VecCUDARestoreArray() returns.

2149:    Fortran note:
2150:    This function is not currently available from Fortran.

2152:    Level: intermediate

2154: .seealso: VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2155: @*/
2156: PETSC_EXTERN PetscErrorCode VecCUDARestoreArray(Vec v, PetscScalar **a)
2157: {

2162: #if defined(PETSC_HAVE_CUDA)
2163:   v->offloadmask = PETSC_OFFLOAD_GPU;
2164: #endif

2166:   PetscObjectStateIncrease((PetscObject)v);
2167:   return(0);
2168: }

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

2173:    This function is analogous to VecGetArrayRead():  The pointer
2174:    returned by this function points to a consistent view of the vector
2175:    data.  This may involve a copy operation of data from the host to the
2176:    device if the data on the device is out of date.  If the device
2177:    memory hasn't been allocated previously it will be allocated as part
2178:    of this function call.  VecCUDAGetArrayRead() assumes that the
2179:    user will not modify the vector data.  This is analgogous to
2180:    intent(in) in Fortran.

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

2188:    Input Parameter:
2189: .  v - the vector

2191:    Output Parameter:
2192: .  a - the CUDA pointer.

2194:    Fortran note:
2195:    This function is not currently available from Fortran.

2197:    Level: intermediate

2199: .seealso: VecCUDARestoreArrayRead(), VecCUDAGetArray(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2200: @*/
2201: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayRead(Vec v, const PetscScalar **a)
2202: {
2203: #if defined(PETSC_HAVE_CUDA)
2205: #endif

2209: #if defined(PETSC_HAVE_CUDA)
2210:   *a   = 0;
2211:   VecCUDACopyToGPU(v);
2212:   *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2213: #endif
2214:   return(0);
2215: }

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

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

2225:    Input Parameter:
2226: +  v - the vector
2227: -  a - the CUDA device pointer.  This pointer is invalid after
2228:        VecCUDARestoreArrayRead() returns.

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

2233:    Level: intermediate

2235: .seealso: VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2236: @*/
2237: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayRead(Vec v, const PetscScalar **a)
2238: {
2241:   *a = NULL;
2242:   return(0);
2243: }

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

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

2254:    The device pointer needs to be released with
2255:    VecCUDARestoreArrayWrite().  When the pointer is released the
2256:    host data of the vector is marked as out of data.  Subsequent access
2257:    of the host data with e.g. VecGetArray() incurs a device to host data
2258:    transfer.


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

2264:    Output Parameter:
2265: .  a - the CUDA pointer

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

2270:    Level: advanced

2272: .seealso: VecCUDARestoreArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2273: @*/
2274: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayWrite(Vec v, PetscScalar **a)
2275: {
2276: #if defined(PETSC_HAVE_CUDA)
2278: #endif

2282: #if defined(PETSC_HAVE_CUDA)
2283:   *a   = 0;
2284:   VecCUDAAllocateCheck(v);
2285:   *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2286: #endif
2287:   return(0);
2288: }

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

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

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

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

2305:    Level: intermediate

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

2315: #if defined(PETSC_HAVE_CUDA)
2316:   v->offloadmask = PETSC_OFFLOAD_GPU;
2317: #endif

2319:   PetscObjectStateIncrease((PetscObject)v);
2320:   return(0);
2321: }

2323: /*@C
2324:    VecCUDAPlaceArray - Allows one to replace the GPU array in a vector with a
2325:    GPU array provided by the user. This is useful to avoid copying an
2326:    array into a vector.

2328:    Not Collective

2330:    Input Parameters:
2331: +  vec - the vector
2332: -  array - the GPU array

2334:    Notes:
2335:    You can return to the original GPU array with a call to VecCUDAResetArray()
2336:    It is not possible to use VecCUDAPlaceArray() and VecPlaceArray() at the
2337:    same time on the same vector.

2339:    Level: developer

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

2343: @*/
2344: PetscErrorCode VecCUDAPlaceArray(Vec vin,const PetscScalar a[])
2345: {

2350: #if defined(PETSC_HAVE_CUDA)
2351:   VecCUDACopyToGPU(vin);
2352:   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()");
2353:   ((Vec_Seq*)vin->data)->unplacedarray  = (PetscScalar *) ((Vec_CUDA*)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2354:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2355:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2356: #endif
2357:   PetscObjectStateIncrease((PetscObject)vin);
2358:   return(0);
2359: }

2361: /*@C
2362:    VecCUDAReplaceArray - Allows one to replace the GPU array in a vector
2363:    with a GPU array provided by the user. This is useful to avoid copying
2364:    a GPU array into a vector.

2366:    Not Collective

2368:    Input Parameters:
2369: +  vec - the vector
2370: -  array - the GPU array

2372:    Notes:
2373:    This permanently replaces the GPU array and frees the memory associated
2374:    with the old GPU array.

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

2379:    Not supported from Fortran

2381:    Level: developer

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

2385: @*/
2386: PetscErrorCode VecCUDAReplaceArray(Vec vin,const PetscScalar a[])
2387: {
2388: #if defined(PETSC_HAVE_CUDA)
2389:   cudaError_t err;
2390: #endif

2395: #if defined(PETSC_HAVE_CUDA)
2396:   err = cudaFree(((Vec_CUDA*)vin->spptr)->GPUarray);CHKERRCUDA(err);
2397:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2398:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2399: #endif
2400:   PetscObjectStateIncrease((PetscObject)vin);
2401:   return(0);
2402: }

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

2408:    Not Collective

2410:    Input Parameters:
2411: .  vec - the vector

2413:    Level: developer

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

2417: @*/
2418: PetscErrorCode VecCUDAResetArray(Vec vin)
2419: {

2424: #if defined(PETSC_HAVE_CUDA)
2425:   VecCUDACopyToGPU(vin);
2426:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar *) ((Vec_Seq*)vin->data)->unplacedarray;
2427:   ((Vec_Seq*)vin->data)->unplacedarray = 0;
2428:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2429: #endif
2430:   PetscObjectStateIncrease((PetscObject)vin);
2431:   return(0);
2432: }

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

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

2441:     Collective on Vec

2443:     Input Parameters:
2444: +   x - a vector to mimic
2445: -   n - the number of vectors to obtain

2447:     Output Parameters:
2448: +   y - Fortran90 pointer to the array of vectors
2449: -   ierr - error code

2451:     Example of Usage:
2452: .vb
2453:  #include <petsc/finclude/petscvec.h>
2454:     use petscvec

2456:     Vec x
2457:     Vec, pointer :: y(:)
2458:     ....
2459:     call VecDuplicateVecsF90(x,2,y,ierr)
2460:     call VecSet(y(2),alpha,ierr)
2461:     call VecSet(y(2),alpha,ierr)
2462:     ....
2463:     call VecDestroyVecsF90(2,y,ierr)
2464: .ve

2466:     Notes:
2467:     Not yet supported for all F90 compilers

2469:     Use VecDestroyVecsF90() to free the space.

2471:     Level: beginner

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

2475: M*/

2477: /*MC
2478:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2479:     VecGetArrayF90().

2481:     Synopsis:
2482:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2484:     Logically Collective on Vec

2486:     Input Parameters:
2487: +   x - vector
2488: -   xx_v - the Fortran90 pointer to the array

2490:     Output Parameter:
2491: .   ierr - error code

2493:     Example of Usage:
2494: .vb
2495:  #include <petsc/finclude/petscvec.h>
2496:     use petscvec

2498:     PetscScalar, pointer :: xx_v(:)
2499:     ....
2500:     call VecGetArrayF90(x,xx_v,ierr)
2501:     xx_v(3) = a
2502:     call VecRestoreArrayF90(x,xx_v,ierr)
2503: .ve

2505:     Level: beginner

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

2509: M*/

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

2514:     Synopsis:
2515:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

2517:     Collective on Vec

2519:     Input Parameters:
2520: +   n - the number of vectors previously obtained
2521: -   x - pointer to array of vector pointers

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

2526:     Notes:
2527:     Not yet supported for all F90 compilers

2529:     Level: beginner

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

2533: M*/

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

2541:     Synopsis:
2542:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2544:     Logically Collective on Vec

2546:     Input Parameter:
2547: .   x - vector

2549:     Output Parameters:
2550: +   xx_v - the Fortran90 pointer to the array
2551: -   ierr - error code

2553:     Example of Usage:
2554: .vb
2555:  #include <petsc/finclude/petscvec.h>
2556:     use petscvec

2558:     PetscScalar, pointer :: xx_v(:)
2559:     ....
2560:     call VecGetArrayF90(x,xx_v,ierr)
2561:     xx_v(3) = a
2562:     call VecRestoreArrayF90(x,xx_v,ierr)
2563: .ve

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

2567:     Level: beginner

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

2571: M*/

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

2579:     Synopsis:
2580:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2582:     Logically Collective on Vec

2584:     Input Parameter:
2585: .   x - vector

2587:     Output Parameters:
2588: +   xx_v - the Fortran90 pointer to the array
2589: -   ierr - error code

2591:     Example of Usage:
2592: .vb
2593:  #include <petsc/finclude/petscvec.h>
2594:     use petscvec

2596:     PetscScalar, pointer :: xx_v(:)
2597:     ....
2598:     call VecGetArrayReadF90(x,xx_v,ierr)
2599:     a = xx_v(3)
2600:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2601: .ve

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

2605:     Level: beginner

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

2609: M*/

2611: /*MC
2612:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2613:     VecGetArrayReadF90().

2615:     Synopsis:
2616:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2618:     Logically Collective on Vec

2620:     Input Parameters:
2621: +   x - vector
2622: -   xx_v - the Fortran90 pointer to the array

2624:     Output Parameter:
2625: .   ierr - error code

2627:     Example of Usage:
2628: .vb
2629:  #include <petsc/finclude/petscvec.h>
2630:     use petscvec

2632:     PetscScalar, pointer :: xx_v(:)
2633:     ....
2634:     call VecGetArrayReadF90(x,xx_v,ierr)
2635:     a = xx_v(3)
2636:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2637: .ve

2639:     Level: beginner

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

2643: M*/

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

2650:    Logically Collective

2652:    Input Parameter:
2653: +  x - the vector
2654: .  m - first dimension of two dimensional array
2655: .  n - second dimension of two dimensional array
2656: .  mstart - first index you will use in first coordinate direction (often 0)
2657: -  nstart - first index in the second coordinate direction (often 0)

2659:    Output Parameter:
2660: .  a - location to put pointer to the array

2662:    Level: developer

2664:   Notes:
2665:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2666:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2667:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2668:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

2672: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2673:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2674:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2675: @*/
2676: PetscErrorCode  VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2677: {
2679:   PetscInt       i,N;
2680:   PetscScalar    *aa;

2686:   VecGetLocalSize(x,&N);
2687:   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);
2688:   VecGetArray(x,&aa);

2690:   PetscMalloc1(m,a);
2691:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2692:   *a -= mstart;
2693:   return(0);
2694: }

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

2701:    Logically Collective

2703:    Input Parameter:
2704: +  x - the vector
2705: .  m - first dimension of two dimensional array
2706: .  n - second dimension of two dimensional array
2707: .  mstart - first index you will use in first coordinate direction (often 0)
2708: -  nstart - first index in the second coordinate direction (often 0)

2710:    Output Parameter:
2711: .  a - location to put pointer to the array

2713:    Level: developer

2715:   Notes:
2716:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2717:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2718:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2719:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

2725: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2726:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2727:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2728: @*/
2729: PetscErrorCode  VecGetArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2730: {
2732:   PetscInt       i,N;
2733:   PetscScalar    *aa;

2739:   VecGetLocalSize(x,&N);
2740:   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);
2741:   VecGetArrayWrite(x,&aa);

2743:   PetscMalloc1(m,a);
2744:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2745:   *a -= mstart;
2746:   return(0);
2747: }

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

2752:    Logically Collective

2754:    Input Parameters:
2755: +  x - the vector
2756: .  m - first dimension of two dimensional array
2757: .  n - second dimension of the two dimensional array
2758: .  mstart - first index you will use in first coordinate direction (often 0)
2759: .  nstart - first index in the second coordinate direction (often 0)
2760: -  a - location of pointer to array obtained from VecGetArray2d()

2762:    Level: developer

2764:    Notes:
2765:    For regular PETSc vectors this routine does not involve any copies. For
2766:    any special vectors that do not store local vector data in a contiguous
2767:    array, this routine will copy the data back into the underlying
2768:    vector data structure from the array obtained with VecGetArray().

2770:    This routine actually zeros out the a pointer.

2772: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2773:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2774:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2775: @*/
2776: PetscErrorCode  VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2777: {
2779:   void           *dummy;

2785:   dummy = (void*)(*a + mstart);
2786:   PetscFree(dummy);
2787:   VecRestoreArray(x,NULL);
2788:   return(0);
2789: }

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

2794:    Logically Collective

2796:    Input Parameters:
2797: +  x - the vector
2798: .  m - first dimension of two dimensional array
2799: .  n - second dimension of the two dimensional array
2800: .  mstart - first index you will use in first coordinate direction (often 0)
2801: .  nstart - first index in the second coordinate direction (often 0)
2802: -  a - location of pointer to array obtained from VecGetArray2d()

2804:    Level: developer

2806:    Notes:
2807:    For regular PETSc vectors this routine does not involve any copies. For
2808:    any special vectors that do not store local vector data in a contiguous
2809:    array, this routine will copy the data back into the underlying
2810:    vector data structure from the array obtained with VecGetArray().

2812:    This routine actually zeros out the a pointer.

2814: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2815:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2816:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2817: @*/
2818: PetscErrorCode  VecRestoreArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2819: {
2821:   void           *dummy;

2827:   dummy = (void*)(*a + mstart);
2828:   PetscFree(dummy);
2829:   VecRestoreArrayWrite(x,NULL);
2830:   return(0);
2831: }

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

2838:    Logically Collective

2840:    Input Parameter:
2841: +  x - the vector
2842: .  m - first dimension of two dimensional array
2843: -  mstart - first index you will use in first coordinate direction (often 0)

2845:    Output Parameter:
2846: .  a - location to put pointer to the array

2848:    Level: developer

2850:   Notes:
2851:    For a vector obtained from DMCreateLocalVector() mstart are likely
2852:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2853:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2857: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2858:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2859:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2860: @*/
2861: PetscErrorCode  VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2862: {
2864:   PetscInt       N;

2870:   VecGetLocalSize(x,&N);
2871:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2872:   VecGetArray(x,a);
2873:   *a  -= mstart;
2874:   return(0);
2875: }

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

2882:    Logically Collective

2884:    Input Parameter:
2885: +  x - the vector
2886: .  m - first dimension of two dimensional array
2887: -  mstart - first index you will use in first coordinate direction (often 0)

2889:    Output Parameter:
2890: .  a - location to put pointer to the array

2892:    Level: developer

2894:   Notes:
2895:    For a vector obtained from DMCreateLocalVector() mstart are likely
2896:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2897:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2901: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2902:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2903:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2904: @*/
2905: PetscErrorCode  VecGetArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2906: {
2908:   PetscInt       N;

2914:   VecGetLocalSize(x,&N);
2915:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2916:   VecGetArrayWrite(x,a);
2917:   *a  -= mstart;
2918:   return(0);
2919: }

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

2924:    Logically Collective

2926:    Input Parameters:
2927: +  x - the vector
2928: .  m - first dimension of two dimensional array
2929: .  mstart - first index you will use in first coordinate direction (often 0)
2930: -  a - location of pointer to array obtained from VecGetArray21()

2932:    Level: developer

2934:    Notes:
2935:    For regular PETSc vectors this routine does not involve any copies. For
2936:    any special vectors that do not store local vector data in a contiguous
2937:    array, this routine will copy the data back into the underlying
2938:    vector data structure from the array obtained with VecGetArray1d().

2940:    This routine actually zeros out the a pointer.

2942: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2943:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2944:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2945: @*/
2946: PetscErrorCode  VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2947: {

2953:   VecRestoreArray(x,NULL);
2954:   return(0);
2955: }

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

2960:    Logically Collective

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

2968:    Level: developer

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

2976:    This routine actually zeros out the a pointer.

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

2980: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2981:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2982:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2983: @*/
2984: PetscErrorCode  VecRestoreArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2985: {

2991:   VecRestoreArrayWrite(x,NULL);
2992:   return(0);
2993: }

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

3000:    Logically Collective

3002:    Input Parameter:
3003: +  x - the vector
3004: .  m - first dimension of three dimensional array
3005: .  n - second dimension of three dimensional array
3006: .  p - third dimension of three dimensional array
3007: .  mstart - first index you will use in first coordinate direction (often 0)
3008: .  nstart - first index in the second coordinate direction (often 0)
3009: -  pstart - first index in the third coordinate direction (often 0)

3011:    Output Parameter:
3012: .  a - location to put pointer to the array

3014:    Level: developer

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

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

3024: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3025:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3026:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3027: @*/
3028: PetscErrorCode  VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3029: {
3031:   PetscInt       i,N,j;
3032:   PetscScalar    *aa,**b;

3038:   VecGetLocalSize(x,&N);
3039:   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);
3040:   VecGetArray(x,&aa);

3042:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3043:   b    = (PetscScalar**)((*a) + m);
3044:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3045:   for (i=0; i<m; i++)
3046:     for (j=0; j<n; j++)
3047:       b[i*n+j] = aa + i*n*p + j*p - pstart;

3049:   *a -= mstart;
3050:   return(0);
3051: }

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

3058:    Logically Collective

3060:    Input Parameter:
3061: +  x - the vector
3062: .  m - first dimension of three dimensional array
3063: .  n - second dimension of three dimensional array
3064: .  p - third dimension of three dimensional array
3065: .  mstart - first index you will use in first coordinate direction (often 0)
3066: .  nstart - first index in the second coordinate direction (often 0)
3067: -  pstart - first index in the third coordinate direction (often 0)

3069:    Output Parameter:
3070: .  a - location to put pointer to the array

3072:    Level: developer

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

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

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

3084: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3085:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3086:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3087: @*/
3088: PetscErrorCode  VecGetArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3089: {
3091:   PetscInt       i,N,j;
3092:   PetscScalar    *aa,**b;

3098:   VecGetLocalSize(x,&N);
3099:   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);
3100:   VecGetArrayWrite(x,&aa);

3102:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3103:   b    = (PetscScalar**)((*a) + m);
3104:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3105:   for (i=0; i<m; i++)
3106:     for (j=0; j<n; j++)
3107:       b[i*n+j] = aa + i*n*p + j*p - pstart;

3109:   *a -= mstart;
3110:   return(0);
3111: }

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

3116:    Logically Collective

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

3128:    Level: developer

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

3136:    This routine actually zeros out the a pointer.

3138: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3139:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3140:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3141: @*/
3142: PetscErrorCode  VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3143: {
3145:   void           *dummy;

3151:   dummy = (void*)(*a + mstart);
3152:   PetscFree(dummy);
3153:   VecRestoreArray(x,NULL);
3154:   return(0);
3155: }

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

3160:    Logically Collective

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

3172:    Level: developer

3174:    Notes:
3175:    For regular PETSc vectors this routine does not involve any copies. For
3176:    any special vectors that do not store local vector data in a contiguous
3177:    array, this routine will copy the data back into the underlying
3178:    vector data structure from the array obtained with VecGetArray().

3180:    This routine actually zeros out the a pointer.

3182: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3183:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3184:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3185: @*/
3186: PetscErrorCode  VecRestoreArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3187: {
3189:   void           *dummy;

3195:   dummy = (void*)(*a + mstart);
3196:   PetscFree(dummy);
3197:   VecRestoreArrayWrite(x,NULL);
3198:   return(0);
3199: }

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

3206:    Logically Collective

3208:    Input Parameter:
3209: +  x - the vector
3210: .  m - first dimension of four dimensional array
3211: .  n - second dimension of four dimensional array
3212: .  p - third dimension of four dimensional array
3213: .  q - fourth dimension of four dimensional array
3214: .  mstart - first index you will use in first coordinate direction (often 0)
3215: .  nstart - first index in the second coordinate direction (often 0)
3216: .  pstart - first index in the third coordinate direction (often 0)
3217: -  qstart - first index in the fourth coordinate direction (often 0)

3219:    Output Parameter:
3220: .  a - location to put pointer to the array

3222:    Level: beginner

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

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

3232: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3233:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3234:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3235: @*/
3236: PetscErrorCode  VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3237: {
3239:   PetscInt       i,N,j,k;
3240:   PetscScalar    *aa,***b,**c;

3246:   VecGetLocalSize(x,&N);
3247:   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);
3248:   VecGetArray(x,&aa);

3250:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3251:   b    = (PetscScalar***)((*a) + m);
3252:   c    = (PetscScalar**)(b + m*n);
3253:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3254:   for (i=0; i<m; i++)
3255:     for (j=0; j<n; j++)
3256:       b[i*n+j] = c + i*n*p + j*p - pstart;
3257:   for (i=0; i<m; i++)
3258:     for (j=0; j<n; j++)
3259:       for (k=0; k<p; k++)
3260:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3261:   *a -= mstart;
3262:   return(0);
3263: }

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

3270:    Logically Collective

3272:    Input Parameter:
3273: +  x - the vector
3274: .  m - first dimension of four dimensional array
3275: .  n - second dimension of four dimensional array
3276: .  p - third dimension of four dimensional array
3277: .  q - fourth dimension of four dimensional array
3278: .  mstart - first index you will use in first coordinate direction (often 0)
3279: .  nstart - first index in the second coordinate direction (often 0)
3280: .  pstart - first index in the third coordinate direction (often 0)
3281: -  qstart - first index in the fourth coordinate direction (often 0)

3283:    Output Parameter:
3284: .  a - location to put pointer to the array

3286:    Level: beginner

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

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

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

3298: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3299:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3300:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3301: @*/
3302: PetscErrorCode  VecGetArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3303: {
3305:   PetscInt       i,N,j,k;
3306:   PetscScalar    *aa,***b,**c;

3312:   VecGetLocalSize(x,&N);
3313:   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);
3314:   VecGetArrayWrite(x,&aa);

3316:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3317:   b    = (PetscScalar***)((*a) + m);
3318:   c    = (PetscScalar**)(b + m*n);
3319:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3320:   for (i=0; i<m; i++)
3321:     for (j=0; j<n; j++)
3322:       b[i*n+j] = c + i*n*p + j*p - pstart;
3323:   for (i=0; i<m; i++)
3324:     for (j=0; j<n; j++)
3325:       for (k=0; k<p; k++)
3326:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3327:   *a -= mstart;
3328:   return(0);
3329: }

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

3334:    Logically Collective

3336:    Input Parameters:
3337: +  x - the vector
3338: .  m - first dimension of four dimensional array
3339: .  n - second dimension of the four dimensional array
3340: .  p - third dimension of the four dimensional array
3341: .  q - fourth dimension of the four dimensional array
3342: .  mstart - first index you will use in first coordinate direction (often 0)
3343: .  nstart - first index in the second coordinate direction (often 0)
3344: .  pstart - first index in the third coordinate direction (often 0)
3345: .  qstart - first index in the fourth coordinate direction (often 0)
3346: -  a - location of pointer to array obtained from VecGetArray4d()

3348:    Level: beginner

3350:    Notes:
3351:    For regular PETSc vectors this routine does not involve any copies. For
3352:    any special vectors that do not store local vector data in a contiguous
3353:    array, this routine will copy the data back into the underlying
3354:    vector data structure from the array obtained with VecGetArray().

3356:    This routine actually zeros out the a pointer.

3358: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3359:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3360:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3361: @*/
3362: PetscErrorCode  VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3363: {
3365:   void           *dummy;

3371:   dummy = (void*)(*a + mstart);
3372:   PetscFree(dummy);
3373:   VecRestoreArray(x,NULL);
3374:   return(0);
3375: }

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

3380:    Logically Collective

3382:    Input Parameters:
3383: +  x - the vector
3384: .  m - first dimension of four dimensional array
3385: .  n - second dimension of the four dimensional array
3386: .  p - third dimension of the four dimensional array
3387: .  q - fourth dimension of the four dimensional array
3388: .  mstart - first index you will use in first coordinate direction (often 0)
3389: .  nstart - first index in the second coordinate direction (often 0)
3390: .  pstart - first index in the third coordinate direction (often 0)
3391: .  qstart - first index in the fourth coordinate direction (often 0)
3392: -  a - location of pointer to array obtained from VecGetArray4d()

3394:    Level: beginner

3396:    Notes:
3397:    For regular PETSc vectors this routine does not involve any copies. For
3398:    any special vectors that do not store local vector data in a contiguous
3399:    array, this routine will copy the data back into the underlying
3400:    vector data structure from the array obtained with VecGetArray().

3402:    This routine actually zeros out the a pointer.

3404: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3405:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3406:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3407: @*/
3408: PetscErrorCode  VecRestoreArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3409: {
3411:   void           *dummy;

3417:   dummy = (void*)(*a + mstart);
3418:   PetscFree(dummy);
3419:   VecRestoreArrayWrite(x,NULL);
3420:   return(0);
3421: }

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

3428:    Logically Collective

3430:    Input Parameter:
3431: +  x - the vector
3432: .  m - first dimension of two dimensional array
3433: .  n - second dimension of two dimensional array
3434: .  mstart - first index you will use in first coordinate direction (often 0)
3435: -  nstart - first index in the second coordinate direction (often 0)

3437:    Output Parameter:
3438: .  a - location to put pointer to the array

3440:    Level: developer

3442:   Notes:
3443:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
3444:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3445:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3446:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

3450: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3451:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3452:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3453: @*/
3454: PetscErrorCode  VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3455: {
3456:   PetscErrorCode    ierr;
3457:   PetscInt          i,N;
3458:   const PetscScalar *aa;

3464:   VecGetLocalSize(x,&N);
3465:   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);
3466:   VecGetArrayRead(x,&aa);

3468:   PetscMalloc1(m,a);
3469:   for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
3470:   *a -= mstart;
3471:   return(0);
3472: }

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

3477:    Logically Collective

3479:    Input Parameters:
3480: +  x - the vector
3481: .  m - first dimension of two dimensional array
3482: .  n - second dimension of the two dimensional array
3483: .  mstart - first index you will use in first coordinate direction (often 0)
3484: .  nstart - first index in the second coordinate direction (often 0)
3485: -  a - location of pointer to array obtained from VecGetArray2d()

3487:    Level: developer

3489:    Notes:
3490:    For regular PETSc vectors this routine does not involve any copies. For
3491:    any special vectors that do not store local vector data in a contiguous
3492:    array, this routine will copy the data back into the underlying
3493:    vector data structure from the array obtained with VecGetArray().

3495:    This routine actually zeros out the a pointer.

3497: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3498:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3499:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3500: @*/
3501: PetscErrorCode  VecRestoreArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3502: {
3504:   void           *dummy;

3510:   dummy = (void*)(*a + mstart);
3511:   PetscFree(dummy);
3512:   VecRestoreArrayRead(x,NULL);
3513:   return(0);
3514: }

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

3521:    Logically Collective

3523:    Input Parameter:
3524: +  x - the vector
3525: .  m - first dimension of two dimensional array
3526: -  mstart - first index you will use in first coordinate direction (often 0)

3528:    Output Parameter:
3529: .  a - location to put pointer to the array

3531:    Level: developer

3533:   Notes:
3534:    For a vector obtained from DMCreateLocalVector() mstart are likely
3535:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3536:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3540: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3541:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3542:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3543: @*/
3544: PetscErrorCode  VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3545: {
3547:   PetscInt       N;

3553:   VecGetLocalSize(x,&N);
3554:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3555:   VecGetArrayRead(x,(const PetscScalar**)a);
3556:   *a  -= mstart;
3557:   return(0);
3558: }

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

3563:    Logically Collective

3565:    Input Parameters:
3566: +  x - the vector
3567: .  m - first dimension of two dimensional array
3568: .  mstart - first index you will use in first coordinate direction (often 0)
3569: -  a - location of pointer to array obtained from VecGetArray21()

3571:    Level: developer

3573:    Notes:
3574:    For regular PETSc vectors this routine does not involve any copies. For
3575:    any special vectors that do not store local vector data in a contiguous
3576:    array, this routine will copy the data back into the underlying
3577:    vector data structure from the array obtained with VecGetArray1dRead().

3579:    This routine actually zeros out the a pointer.

3581: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3582:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3583:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3584: @*/
3585: PetscErrorCode  VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3586: {

3592:   VecRestoreArrayRead(x,NULL);
3593:   return(0);
3594: }


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

3602:    Logically Collective

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

3613:    Output Parameter:
3614: .  a - location to put pointer to the array

3616:    Level: developer

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

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

3626: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3627:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3628:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3629: @*/
3630: PetscErrorCode  VecGetArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3631: {
3632:   PetscErrorCode    ierr;
3633:   PetscInt          i,N,j;
3634:   const PetscScalar *aa;
3635:   PetscScalar       **b;

3641:   VecGetLocalSize(x,&N);
3642:   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);
3643:   VecGetArrayRead(x,&aa);

3645:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3646:   b    = (PetscScalar**)((*a) + m);
3647:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3648:   for (i=0; i<m; i++)
3649:     for (j=0; j<n; j++)
3650:       b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;

3652:   *a -= mstart;
3653:   return(0);
3654: }

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

3659:    Logically Collective

3661:    Input Parameters:
3662: +  x - the vector
3663: .  m - first dimension of three dimensional array
3664: .  n - second dimension of the three dimensional array
3665: .  p - third dimension of the three dimensional array
3666: .  mstart - first index you will use in first coordinate direction (often 0)
3667: .  nstart - first index in the second coordinate direction (often 0)
3668: .  pstart - first index in the third coordinate direction (often 0)
3669: -  a - location of pointer to array obtained from VecGetArray3dRead()

3671:    Level: developer

3673:    Notes:
3674:    For regular PETSc vectors this routine does not involve any copies. For
3675:    any special vectors that do not store local vector data in a contiguous
3676:    array, this routine will copy the data back into the underlying
3677:    vector data structure from the array obtained with VecGetArray().

3679:    This routine actually zeros out the a pointer.

3681: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3682:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3683:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3684: @*/
3685: PetscErrorCode  VecRestoreArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3686: {
3688:   void           *dummy;

3694:   dummy = (void*)(*a + mstart);
3695:   PetscFree(dummy);
3696:   VecRestoreArrayRead(x,NULL);
3697:   return(0);
3698: }

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

3705:    Logically Collective

3707:    Input Parameter:
3708: +  x - the vector
3709: .  m - first dimension of four dimensional array
3710: .  n - second dimension of four dimensional array
3711: .  p - third dimension of four dimensional array
3712: .  q - fourth dimension of four dimensional array
3713: .  mstart - first index you will use in first coordinate direction (often 0)
3714: .  nstart - first index in the second coordinate direction (often 0)
3715: .  pstart - first index in the third coordinate direction (often 0)
3716: -  qstart - first index in the fourth coordinate direction (often 0)

3718:    Output Parameter:
3719: .  a - location to put pointer to the array

3721:    Level: beginner

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

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

3731: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3732:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3733:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3734: @*/
3735: PetscErrorCode  VecGetArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3736: {
3737:   PetscErrorCode    ierr;
3738:   PetscInt          i,N,j,k;
3739:   const PetscScalar *aa;
3740:   PetscScalar       ***b,**c;

3746:   VecGetLocalSize(x,&N);
3747:   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);
3748:   VecGetArrayRead(x,&aa);

3750:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3751:   b    = (PetscScalar***)((*a) + m);
3752:   c    = (PetscScalar**)(b + m*n);
3753:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3754:   for (i=0; i<m; i++)
3755:     for (j=0; j<n; j++)
3756:       b[i*n+j] = c + i*n*p + j*p - pstart;
3757:   for (i=0; i<m; i++)
3758:     for (j=0; j<n; j++)
3759:       for (k=0; k<p; k++)
3760:         c[i*n*p+j*p+k] = (PetscScalar*) aa + i*n*p*q + j*p*q + k*q - qstart;
3761:   *a -= mstart;
3762:   return(0);
3763: }

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

3768:    Logically Collective

3770:    Input Parameters:
3771: +  x - the vector
3772: .  m - first dimension of four dimensional array
3773: .  n - second dimension of the four dimensional array
3774: .  p - third dimension of the four dimensional array
3775: .  q - fourth dimension of the four dimensional array
3776: .  mstart - first index you will use in first coordinate direction (often 0)
3777: .  nstart - first index in the second coordinate direction (often 0)
3778: .  pstart - first index in the third coordinate direction (often 0)
3779: .  qstart - first index in the fourth coordinate direction (often 0)
3780: -  a - location of pointer to array obtained from VecGetArray4dRead()

3782:    Level: beginner

3784:    Notes:
3785:    For regular PETSc vectors this routine does not involve any copies. For
3786:    any special vectors that do not store local vector data in a contiguous
3787:    array, this routine will copy the data back into the underlying
3788:    vector data structure from the array obtained with VecGetArray().

3790:    This routine actually zeros out the a pointer.

3792: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3793:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3794:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3795: @*/
3796: PetscErrorCode  VecRestoreArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3797: {
3799:   void           *dummy;

3805:   dummy = (void*)(*a + mstart);
3806:   PetscFree(dummy);
3807:   VecRestoreArrayRead(x,NULL);
3808:   return(0);
3809: }

3811: #if defined(PETSC_USE_DEBUG)

3813: /*@
3814:    VecLockGet  - Gets the current lock status of a vector

3816:    Logically Collective on Vec

3818:    Input Parameter:
3819: .  x - the vector

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

3825:    Level: beginner

3827: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop()
3828: @*/
3829: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
3830: {
3833:   *state = x->lock;
3834:   return(0);
3835: }

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

3840:    Logically Collective on Vec

3842:    Input Parameter:
3843: .  x - the vector

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

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

3851:    Level: beginner

3853: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPop(), VecLockGet()
3854: @*/
3855: PetscErrorCode VecLockReadPush(Vec x)
3856: {
3859:   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");
3860:   x->lock++;
3861:   return(0);
3862: }

3864: /*@
3865:    VecLockReadPop  - Pops a read-only lock from a vector

3867:    Logically Collective on Vec

3869:    Input Parameter:
3870: .  x - the vector

3872:    Level: beginner

3874: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockGet()
3875: @*/
3876: PetscErrorCode VecLockReadPop(Vec x)
3877: {
3880:   x->lock--;
3881:   if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector has been unlocked from read-only access too many times");
3882:   return(0);
3883: }

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

3888:    Logically Collective on Vec

3890:    Input Parameter:
3891: +  x   - the vector
3892: -  flg - PETSC_TRUE to lock the vector for writing; PETSC_FALSE to unlock it.

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


3901:        VecGetArray(x,&xdata); // begin phase
3902:        VecLockWriteSet_Private(v,PETSC_TRUE);

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

3906:        VecRestoreArray(x,&vdata); // end phase
3907:        VecLockWriteSet_Private(v,PETSC_FALSE);

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

3912:    Level: beginner

3914: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop(), VecLockGet()
3915: @*/
3916: PetscErrorCode VecLockWriteSet_Private(Vec x,PetscBool flg)
3917: {
3920:   if (flg) {
3921:     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");
3922:     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");
3923:     else x->lock = -1;
3924:   } else {
3925:     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");
3926:     x->lock = 0;
3927:   }
3928:   return(0);
3929: }

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

3934:    Level: deprecated

3936: .seealso: VecLockReadPush()
3937: @*/
3938: PetscErrorCode VecLockPush(Vec x)
3939: {
3942:   VecLockReadPush(x);
3943:   return(0);
3944: }

3946: /*@
3947:    VecLockPop  - Pops a read-only lock from a vector

3949:    Level: deprecated

3951: .seealso: VecLockReadPop()
3952: @*/
3953: PetscErrorCode VecLockPop(Vec x)
3954: {
3957:   VecLockReadPop(x);
3958:   return(0);
3959: }

3961: #endif