Actual source code: rvector.c

petsc-master 2020-01-26
Report Typos and Errors

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

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

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

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

 46:    Logically Collective on Vec

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

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

 54:    Level: advanced

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

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

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

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

 81:    Collective on Vec

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

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

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

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

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

105:    Level: intermediate


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

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

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

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

132:    Collective on Vec

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

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

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

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

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

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

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

155:    Level: intermediate


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

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

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

174:    Collective on Vec

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

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

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

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

195:    Level: intermediate

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


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

206: @*/

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


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

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

237:    Not Collective

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

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

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

254:    Level: intermediate

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

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


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

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


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

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

290:    Collective on Vec

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

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

299:    Level: intermediate


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

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

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

327:    Collective on Vec

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

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

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

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


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

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

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

362:    Collective on Vec

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

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

371:    Level: intermediate

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

376:    This returns the smallest index with the minumum value


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

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

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

399:    Collective on Vec

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

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

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

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

416:    Level: intermediate

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

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

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

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

442:    Not collective on Vec

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

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

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

455:    Level: intermediate


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

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

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

493:    Logically Collective on Vec

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

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

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

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

512:    Level: beginner

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

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

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

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

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


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

557:    Logically Collective on Vec

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

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

566:    Level: intermediate

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

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


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

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

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

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

610:    Logically Collective on Vec

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

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

619:    Level: intermediate

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


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

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

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

654:    Logically Collective on Vec

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

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

663:    Level: intermediate

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


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

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

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

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

705:    Logically Collective on Vec

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

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

714:    Level: intermediate

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


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

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

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


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

749:    Logically Collective on Vec

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

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

758:    Level: intermediate

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


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

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

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


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

798:    Not Collective

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

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

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

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

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

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

827:    Level: beginner

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

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

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

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

854:     Not Collective

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

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

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

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

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

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

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

877:    Level: beginner

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

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

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

898:    Not Collective

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

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

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

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

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

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

927:    Level: intermediate

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

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

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


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

955:    Not Collective

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

966:    Level: intermediate

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

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

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

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

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

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

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

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

1018:    Not Collective

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

1029:    Level: intermediate

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

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

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

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


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

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

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

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

1078:    Collective on Vec

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

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

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

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

1097:    Level: intermediate


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

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

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

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

1127:    Collective on Vec

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

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

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

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

1146:    Level: intermediate


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

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

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

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

1177:    Logically Collective on Vec

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

1185:    Level: intermediate

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

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

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

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

1224:    Collective on IS

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

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

1233:    Level: advanced

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

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

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

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

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

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

1314:    Collective on IS

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

1321:    Level: advanced

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

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

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

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

1359:    Not collective.

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

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

1367:    Level: beginner

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

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

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

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

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

1407:    Not collective.

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

1413:    Level: beginner

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

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

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

1439:    Collective on v, not collective on w.

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

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

1447:    Level: beginner

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

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

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

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

1484:    Logically collective.

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

1490:    Level: beginner

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

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

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

1521:    Logically Collective on Vec

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

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

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

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

1547:    Level: beginner

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

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

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

1604:    Logically Collective on Vec

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

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

1612:    Level: beginner

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

1623:   VecSetErrorIfLocked(x,1);

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

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

1645:    Logically Collective on Vec

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

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

1653:    Level: intermediate

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

1658:    Concepts: vector^accessing local values

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

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

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

1681:    Not Collective

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

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

1689:    Level: beginner

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

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

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

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

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

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

1741:    Not Collective

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

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

1749:    Level: beginner

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


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

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

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

1782:    Logically Collective on Vec

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

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

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

1794:    Level: intermediate

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

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

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

1821:    Logically Collective on Vec

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

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

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

1837:    Level: intermediate

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


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

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

1862:    Logically Collective on Vec

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

1868:    Level: beginner

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

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

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

1894:    Logically Collective on Vec

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

1900:    Level: beginner

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

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


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

1929:    Logically Collective on Vec

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

1935:    Level: beginner

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

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

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

1965:    Not Collective

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

1971:    Level: beginner

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

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

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

1995:    Not Collective

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

2001:    Level: beginner

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

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

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

2019:    Not Collective

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

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

2028:    Level: developer

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

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

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

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

2053:    Not Collective

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

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

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

2066:    Not supported from Fortran

2068:    Level: developer

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

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

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


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

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

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


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

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

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

2116:    Level: intermediate

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

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

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

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

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

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

2151:    Level: intermediate

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

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

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

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

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

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

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

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

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

2196:    Level: intermediate

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

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

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

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

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

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

2232:    Level: intermediate

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

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

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

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


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

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

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

2269:    Level: advanced

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

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

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

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

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

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

2304:    Level: intermediate

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

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

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

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

2327:    Not Collective

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

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

2338:    Level: developer

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

2342: @*/
2343: PetscErrorCode VecCUDAPlaceArray(Vec vin,PetscScalar *a)
2344: {

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

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

2365:    Not Collective

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

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

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

2378:    Not supported from Fortran

2380:    Level: developer

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

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

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

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

2407:    Not Collective

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

2412:    Level: developer

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

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

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

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

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

2440:     Collective on Vec

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

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

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

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

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

2468:     Use VecDestroyVecsF90() to free the space.

2470:     Level: beginner

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

2474: M*/

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

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

2483:     Logically Collective on Vec

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

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

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

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

2504:     Level: beginner

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

2508: M*/

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

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

2516:     Collective on Vec

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

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

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

2528:     Level: beginner

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

2532: M*/

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

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

2543:     Logically Collective on Vec

2545:     Input Parameter:
2546: .   x - vector

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

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

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

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

2566:     Level: beginner

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

2570: M*/

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

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

2581:     Logically Collective on Vec

2583:     Input Parameter:
2584: .   x - vector

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

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

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

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

2604:     Level: beginner

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

2608: M*/

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

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

2617:     Logically Collective on Vec

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

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

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

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

2638:     Level: beginner

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

2642: M*/

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

2649:    Logically Collective

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

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

2661:    Level: developer

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

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

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

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

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

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

2700:    Logically Collective

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

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

2712:    Level: developer

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

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

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

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

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

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

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

2751:    Logically Collective

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

2761:    Level: developer

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

2769:    This routine actually zeros out the a pointer.

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

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

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

2793:    Logically Collective

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

2803:    Level: developer

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

2811:    This routine actually zeros out the a pointer.

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

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

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

2837:    Logically Collective

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

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

2847:    Level: developer

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

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

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

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

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

2881:    Logically Collective

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

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

2891:    Level: developer

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

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

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

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

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

2923:    Logically Collective

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

2931:    Level: developer

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

2939:    This routine actually zeros out the a pointer.

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

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

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

2959:    Logically Collective

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

2967:    Level: developer

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

2975:    This routine actually zeros out the a pointer.

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

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

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

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

2999:    Logically Collective

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

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

3013:    Level: developer

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

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

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

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

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

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

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

3057:    Logically Collective

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

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

3071:    Level: developer

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

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

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

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

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

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

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

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

3115:    Logically Collective

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

3127:    Level: developer

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

3135:    This routine actually zeros out the a pointer.

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

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

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

3159:    Logically Collective

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

3171:    Level: developer

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

3179:    This routine actually zeros out the a pointer.

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

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

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

3205:    Logically Collective

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

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

3221:    Level: beginner

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

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

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

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

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

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

3269:    Logically Collective

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

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

3285:    Level: beginner

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

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

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

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

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

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

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

3333:    Logically Collective

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

3347:    Level: beginner

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

3355:    This routine actually zeros out the a pointer.

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

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

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

3379:    Logically Collective

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

3393:    Level: beginner

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

3401:    This routine actually zeros out the a pointer.

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

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

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

3427:    Logically Collective

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

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

3439:    Level: developer

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

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

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

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

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

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

3476:    Logically Collective

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

3486:    Level: developer

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

3494:    This routine actually zeros out the a pointer.

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

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

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

3520:    Logically Collective

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

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

3530:    Level: developer

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

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

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

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

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

3562:    Logically Collective

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

3570:    Level: developer

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

3578:    This routine actually zeros out the a pointer.

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

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


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

3601:    Logically Collective

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

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

3615:    Level: developer

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

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

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

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

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

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

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

3658:    Logically Collective

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

3670:    Level: developer

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

3678:    This routine actually zeros out the a pointer.

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

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

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

3704:    Logically Collective

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

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

3720:    Level: beginner

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

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

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

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

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

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

3767:    Logically Collective

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

3781:    Level: beginner

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

3789:    This routine actually zeros out the a pointer.

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

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

3810: #if defined(PETSC_USE_DEBUG)

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

3815:    Logically Collective on Vec

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

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

3824:    Level: beginner

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

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

3839:    Logically Collective on Vec

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

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

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

3850:    Level: beginner

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

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

3866:    Logically Collective on Vec

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

3871:    Level: beginner

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

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

3887:    Logically Collective on Vec

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

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


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

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

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

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

3911:    Level: beginner

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

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

3933:    Level: deprecated

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

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

3948:    Level: deprecated

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

3960: #endif