Actual source code: rvector.c

petsc-master 2019-12-03
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:     /* get current stashed norms */
473:     for (i=0; i<4; i++) {
474:       PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
475:     }
476:     (*x->ops->scale)(x,alpha);
477:     PetscObjectStateIncrease((PetscObject)x);
478:     /* put the scaled stashed norms back into the Vec */
479:     for (i=0; i<4; i++) {
480:       if (flgs[i]) {
481:         PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],PetscAbsScalar(alpha)*norms[i]);
482:       }
483:     }
484:   }
485:   PetscLogEventEnd(VEC_Scale,x,0,0,0);
486:   return(0);
487: }

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

492:    Logically Collective on Vec

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

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

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

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

511:    Level: beginner

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

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

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

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

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


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

556:    Logically Collective on Vec

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

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

565:    Level: intermediate

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

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


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

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

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

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

608:    Logically Collective on Vec

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

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

617:    Level: intermediate

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


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

636:   VecCheckSameSize(y,1,x,4);
637:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
640:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
641:   (*y->ops->axpby)(y,alpha,beta,x);
642:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
643:   PetscObjectStateIncrease((PetscObject)y);
644:   return(0);
645: }

647: /*@
648:    VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z

650:    Logically Collective on Vec

652:    Input Parameters:
653: +  alpha,beta, gamma - the scalars
654: -  x, y, z  - the vectors

656:    Output Parameter:
657: .  z - output vector

659:    Level: intermediate

661:    Notes:
662:     x, y and z must be different vectors
663:     The implementation is optimized for alpha of 1.0 and gamma of 1.0 or 0.0


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

681:   VecCheckSameSize(x,1,y,5);
682:   VecCheckSameSize(x,1,z,6);
683:   if (x == y || x == z) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
684:   if (y == z) SETERRQ(PetscObjectComm((PetscObject)y),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");

689:   PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
690:   (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
691:   PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
692:   PetscObjectStateIncrease((PetscObject)z);
693:   return(0);
694: }

696: /*@
697:    VecAYPX - Computes y = x + beta y.

699:    Logically Collective on Vec

701:    Input Parameters:
702: +  beta - the scalar
703: -  x, y  - the vectors

705:    Output Parameter:
706: .  y - output vector

708:    Level: intermediate

710:    Notes:
711:     x and y MUST be different vectors
712:     The implementation is optimized for beta of -1.0, 0.0, and 1.0


715: .seealso:  VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
716: @*/
717: PetscErrorCode  VecAYPX(Vec y,PetscScalar beta,Vec x)
718: {

727:   VecCheckSameSize(x,1,y,3);
728:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y must be different vectors");

731:   PetscLogEventBegin(VEC_AYPX,x,y,0,0);
732:    (*y->ops->aypx)(y,beta,x);
733:   PetscLogEventEnd(VEC_AYPX,x,y,0,0);
734:   PetscObjectStateIncrease((PetscObject)y);
735:   return(0);
736: }


739: /*@
740:    VecWAXPY - Computes w = alpha x + y.

742:    Logically Collective on Vec

744:    Input Parameters:
745: +  alpha - the scalar
746: -  x, y  - the vectors

748:    Output Parameter:
749: .  w - the result

751:    Level: intermediate

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


758: .seealso: VecAXPY(), VecAYPX(), VecAXPBY(), VecMAXPY(), VecAXPBYPCZ()
759: @*/
760: PetscErrorCode  VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
761: {

773:   VecCheckSameSize(x,3,y,4);
774:   VecCheckSameSize(x,3,w,1);
775:   if (w == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector y, suggest VecAXPY()");
776:   if (w == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector x, suggest VecAYPX()");

779:   PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
780:    (*w->ops->waxpy)(w,alpha,x,y);
781:   PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
782:   PetscObjectStateIncrease((PetscObject)w);
783:   return(0);
784: }


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

790:    Not Collective

792:    Input Parameters:
793: +  x - vector to insert in
794: .  ni - number of elements to add
795: .  ix - indices where to add
796: .  y - array of values
797: -  iora - either INSERT_VALUES or ADD_VALUES, where
798:    ADD_VALUES adds values to any existing entries, and
799:    INSERT_VALUES replaces existing entries with new values

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

804:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
805:    options cannot be mixed without intervening calls to the assembly
806:    routines.

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

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

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

819:    Level: beginner

821: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
822:            VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
823: @*/
824: PetscErrorCode  VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
825: {

830:   if (!ni) return(0);
834:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
835:   (*x->ops->setvalues)(x,ni,ix,y,iora);
836:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
837:   PetscObjectStateIncrease((PetscObject)x);
838:   return(0);
839: }

841: /*@
842:    VecGetValues - Gets values from certain locations of a vector. Currently
843:           can only get values on the same processor

845:     Not Collective

847:    Input Parameters:
848: +  x - vector to get values from
849: .  ni - number of elements to get
850: -  ix - indices where to get them from (in global 1d numbering)

852:    Output Parameter:
853: .   y - array of values

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

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

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

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

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

868:    Level: beginner

870: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues()
871: @*/
872: PetscErrorCode  VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
873: {

878:   if (!ni) return(0);
882:   (*x->ops->getvalues)(x,ni,ix,y);
883:   return(0);
884: }

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

889:    Not Collective

891:    Input Parameters:
892: +  x - vector to insert in
893: .  ni - number of blocks to add
894: .  ix - indices where to add in block count, rather than element count
895: .  y - array of values
896: -  iora - either INSERT_VALUES or ADD_VALUES, where
897:    ADD_VALUES adds values to any existing entries, and
898:    INSERT_VALUES replaces existing entries with new values

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

904:    Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
905:    options cannot be mixed without intervening calls to the assembly
906:    routines.

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

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

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

918:    Level: intermediate

920: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
921:            VecSetValues()
922: @*/
923: PetscErrorCode  VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
924: {

929:   if (!ni) return(0);
933:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
934:   (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
935:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
936:   PetscObjectStateIncrease((PetscObject)x);
937:   return(0);
938: }


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

945:    Not Collective

947:    Input Parameters:
948: +  x - vector to insert in
949: .  ni - number of elements to add
950: .  ix - indices where to add
951: .  y - array of values
952: -  iora - either INSERT_VALUES or ADD_VALUES, where
953:    ADD_VALUES adds values to any existing entries, and
954:    INSERT_VALUES replaces existing entries with new values

956:    Level: intermediate

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

961:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
962:    options cannot be mixed without intervening calls to the assembly
963:    routines.

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

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

970: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetLocalToGlobalMapping(),
971:            VecSetValuesBlockedLocal()
972: @*/
973: PetscErrorCode  VecSetValuesLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
974: {
976:   PetscInt       lixp[128],*lix = lixp;

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

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

1004: /*@
1005:    VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1006:    using a local ordering of the nodes.

1008:    Not Collective

1010:    Input Parameters:
1011: +  x - vector to insert in
1012: .  ni - number of blocks to add
1013: .  ix - indices where to add in block count, not element count
1014: .  y - array of values
1015: -  iora - either INSERT_VALUES or ADD_VALUES, where
1016:    ADD_VALUES adds values to any existing entries, and
1017:    INSERT_VALUES replaces existing entries with new values

1019:    Level: intermediate

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

1025:    Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1026:    options cannot be mixed without intervening calls to the assembly
1027:    routines.

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

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


1035: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
1036:            VecSetLocalToGlobalMapping()
1037: @*/
1038: PetscErrorCode  VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1039: {
1041:   PetscInt       lixp[128],*lix = lixp;

1045:   if (!ni) return(0);
1049:   if (ni > 128) {
1050:     PetscMalloc1(ni,&lix);
1051:   }

1053:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1054:   ISLocalToGlobalMappingApplyBlock(x->map->mapping,ni,(PetscInt*)ix,lix);
1055:   (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1056:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1057:   if (ni > 128) {
1058:     PetscFree(lix);
1059:   }
1060:   PetscObjectStateIncrease((PetscObject)x);
1061:   return(0);
1062: }

1064: /*@
1065:    VecMTDot - Computes indefinite vector multiple dot products.
1066:    That is, it does NOT use the complex conjugate.

1068:    Collective on Vec

1070:    Input Parameters:
1071: +  x - one vector
1072: .  nv - number of vectors
1073: -  y - array of vectors.  Note that vectors are pointers

1075:    Output Parameter:
1076: .  val - array of the dot products

1078:    Notes for Users of Complex Numbers:
1079:    For complex vectors, VecMTDot() computes the indefinite form
1080: $      val = (x,y) = y^T x,
1081:    where y^T denotes the transpose of y.

1083:    Use VecMDot() for the inner product
1084: $      val = (x,y) = y^H x,
1085:    where y^H denotes the conjugate transpose of y.

1087:    Level: intermediate


1090: .seealso: VecMDot(), VecTDot()
1091: @*/
1092: PetscErrorCode  VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1093: {

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

1108:   PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1109:   (*x->ops->mtdot)(x,nv,y,val);
1110:   PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1111:   return(0);
1112: }

1114: /*@
1115:    VecMDot - Computes vector multiple dot products.

1117:    Collective on Vec

1119:    Input Parameters:
1120: +  x - one vector
1121: .  nv - number of vectors
1122: -  y - array of vectors.

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

1127:    Notes for Users of Complex Numbers:
1128:    For complex vectors, VecMDot() computes
1129: $     val = (x,y) = y^H x,
1130:    where y^H denotes the conjugate transpose of y.

1132:    Use VecMTDot() for the indefinite form
1133: $     val = (x,y) = y^T x,
1134:    where y^T denotes the transpose of y.

1136:    Level: intermediate


1139: .seealso: VecMTDot(), VecDot()
1140: @*/
1141: PetscErrorCode  VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1142: {

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

1158:   PetscLogEventBegin(VEC_MDot,x,*y,0,0);
1159:   (*x->ops->mdot)(x,nv,y,val);
1160:   PetscLogEventEnd(VEC_MDot,x,*y,0,0);
1161:   return(0);
1162: }

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

1167:    Logically Collective on Vec

1169:    Input Parameters:
1170: +  nv - number of scalars and x-vectors
1171: .  alpha - array of scalars
1172: .  y - one vector
1173: -  x - array of vectors

1175:    Level: intermediate

1177:    Notes:
1178:     y cannot be any of the x vectors

1180: .seealso:  VecAYPX(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
1181: @*/
1182: PetscErrorCode  VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1183: {
1185:   PetscInt       i;

1190:   if (!nv) return(0);
1191:   if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1198:   VecCheckSameSize(y,1,*x,4);

1201:   PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1202:   (*y->ops->maxpy)(y,nv,alpha,x);
1203:   PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1204:   PetscObjectStateIncrease((PetscObject)y);
1205:   return(0);
1206: }

1208: /*@
1209:    VecGetSubVector - Gets a vector representing part of another vector

1211:    Collective on IS

1213:    Input Arguments:
1214: + X - vector from which to extract a subvector
1215: - is - index set representing portion of X to extract

1217:    Output Arguments:
1218: . Y - subvector corresponding to is

1220:    Level: advanced

1222:    Notes:
1223:    The subvector Y should be returned with VecRestoreSubVector().

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

1228: .seealso: MatCreateSubMatrix()
1229: @*/
1230: PetscErrorCode  VecGetSubVector(Vec X,IS is,Vec *Y)
1231: {
1232:   PetscErrorCode   ierr;
1233:   Vec              Z;

1239:   if (X->ops->getsubvector) {
1240:     (*X->ops->getsubvector)(X,is,&Z);
1241:   } else { /* Default implementation currently does no caching */
1242:     PetscInt  gstart,gend,start;
1243:     PetscBool contiguous,gcontiguous;
1244:     VecGetOwnershipRange(X,&gstart,&gend);
1245:     ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1246:     MPIU_Allreduce(&contiguous,&gcontiguous,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1247:     if (gcontiguous) { /* We can do a no-copy implementation */
1248:       PetscInt n,N,bs;
1249:       PetscInt state;

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

1298: /*@
1299:    VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()

1301:    Collective on IS

1303:    Input Arguments:
1304: + X - vector from which subvector was obtained
1305: . is - index set representing the subset of X
1306: - Y - subvector being restored

1308:    Level: advanced

1310: .seealso: VecGetSubVector()
1311: @*/
1312: PetscErrorCode  VecRestoreSubVector(Vec X,IS is,Vec *Y)
1313: {

1321:   if (X->ops->restoresubvector) {
1322:     (*X->ops->restoresubvector)(X,is,Y);
1323:   } else {
1324:     PETSC_UNUSED PetscObjectState dummystate = 0;
1325:     PetscBool valid;
1326:     PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,dummystate,valid);
1327:     if (!valid) {
1328:       VecScatter scatter;

1330:       PetscObjectQuery((PetscObject)*Y,"VecGetSubVector_Scatter",(PetscObject*)&scatter);
1331:       if (scatter) {
1332:         VecScatterBegin(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1333:         VecScatterEnd(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1334:       }
1335:     }
1336:     VecDestroy(Y);
1337:   }
1338:   return(0);
1339: }

1341: /*@
1342:    VecGetLocalVectorRead - Maps the local portion of a vector into a
1343:    vector.  You must call VecRestoreLocalVectorRead() when the local
1344:    vector is no longer needed.

1346:    Not collective.

1348:    Input parameter:
1349: .  v - The vector for which the local vector is desired.

1351:    Output parameter:
1352: .  w - Upon exit this contains the local vector.

1354:    Level: beginner

1356:    Notes:
1357:    This function is similar to VecGetArrayRead() which maps the local
1358:    portion into a raw pointer.  VecGetLocalVectorRead() is usually
1359:    almost as efficient as VecGetArrayRead() but in certain circumstances
1360:    VecGetLocalVectorRead() can be much more efficient than
1361:    VecGetArrayRead().  This is because the construction of a contiguous
1362:    array representing the vector data required by VecGetArrayRead() can
1363:    be an expensive operation for certain vector types.  For example, for
1364:    GPU vectors VecGetArrayRead() requires that the data between device
1365:    and host is synchronized.

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

1370: .seealso: VecRestoreLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1371: @*/
1372: PetscErrorCode VecGetLocalVectorRead(Vec v,Vec w)
1373: {
1375:   PetscScalar    *a;

1380:   VecCheckSameLocalSize(v,1,w,2);
1381:   if (v->ops->getlocalvectorread) {
1382:     (*v->ops->getlocalvectorread)(v,w);
1383:   } else {
1384:     VecGetArrayRead(v,(const PetscScalar**)&a);
1385:     VecPlaceArray(w,a);
1386:   }
1387:   return(0);
1388: }

1390: /*@
1391:    VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1392:    previously mapped into a vector using VecGetLocalVectorRead().

1394:    Not collective.

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

1400:    Level: beginner

1402: .seealso: VecGetLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1403: @*/
1404: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec w)
1405: {
1407:   PetscScalar    *a;

1412:   if (v->ops->restorelocalvectorread) {
1413:     (*v->ops->restorelocalvectorread)(v,w);
1414:   } else {
1415:     VecGetArrayRead(w,(const PetscScalar**)&a);
1416:     VecRestoreArrayRead(v,(const PetscScalar**)&a);
1417:     VecResetArray(w);
1418:   }
1419:   return(0);
1420: }

1422: /*@
1423:    VecGetLocalVector - Maps the local portion of a vector into a
1424:    vector.

1426:    Collective on v, not collective on w.

1428:    Input parameter:
1429: .  v - The vector for which the local vector is desired.

1431:    Output parameter:
1432: .  w - Upon exit this contains the local vector.

1434:    Level: beginner

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

1447: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1448: @*/
1449: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1450: {
1452:   PetscScalar    *a;

1457:   VecCheckSameLocalSize(v,1,w,2);
1458:   if (v->ops->getlocalvector) {
1459:     (*v->ops->getlocalvector)(v,w);
1460:   } else {
1461:     VecGetArray(v,&a);
1462:     VecPlaceArray(w,a);
1463:   }
1464:   return(0);
1465: }

1467: /*@
1468:    VecRestoreLocalVector - Unmaps the local portion of a vector
1469:    previously mapped into a vector using VecGetLocalVector().

1471:    Logically collective.

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

1477:    Level: beginner

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

1489:   if (v->ops->restorelocalvector) {
1490:     (*v->ops->restorelocalvector)(v,w);
1491:   } else {
1492:     VecGetArray(w,&a);
1493:     VecRestoreArray(v,&a);
1494:     VecResetArray(w);
1495:   }
1496:   return(0);
1497: }

1499: /*@C
1500:    VecGetArray - Returns a pointer to a contiguous array that contains this
1501:    processor's portion of the vector data. For the standard PETSc
1502:    vectors, VecGetArray() returns a pointer to the local data array and
1503:    does not use any copies. If the underlying vector data is not stored
1504:    in a contiguous array this routine will copy the data to a contiguous
1505:    array and return a pointer to that. You MUST call VecRestoreArray()
1506:    when you no longer need access to the array.

1508:    Logically Collective on Vec

1510:    Input Parameter:
1511: .  x - the vector

1513:    Output Parameter:
1514: .  a - location to put pointer to the array

1516:    Fortran Note:
1517:    This routine is used differently from Fortran 77
1518: $    Vec         x
1519: $    PetscScalar x_array(1)
1520: $    PetscOffset i_x
1521: $    PetscErrorCode ierr
1522: $       call VecGetArray(x,x_array,i_x,ierr)
1523: $
1524: $   Access first local entry in vector with
1525: $      value = x_array(i_x + 1)
1526: $
1527: $      ...... other code
1528: $       call VecRestoreArray(x,x_array,i_x,ierr)
1529:    For Fortran 90 see VecGetArrayF90()

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

1534:    Level: beginner

1536: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1537:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1538: @*/
1539: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1540: {
1542: #if defined(PETSC_HAVE_VIENNACL)
1543:   PetscBool      is_viennacltype = PETSC_FALSE;
1544: #endif

1548:   VecSetErrorIfLocked(x,1);
1549:   if (x->petscnative) {
1550: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1551:     if (x->offloadmask == PETSC_OFFLOAD_GPU) {
1552: #if defined(PETSC_HAVE_VIENNACL)
1553:       PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1554:       if (is_viennacltype) {
1555:         VecViennaCLCopyFromGPU(x);
1556:       } else
1557: #endif
1558:       {
1559: #if defined(PETSC_HAVE_CUDA)
1560:         VecCUDACopyFromGPU(x);
1561: #endif
1562:       }
1563:     } else if (x->offloadmask == PETSC_OFFLOAD_UNALLOCATED) {
1564: #if defined(PETSC_HAVE_VIENNACL)
1565:       PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1566:       if (is_viennacltype) {
1567:         VecViennaCLAllocateCheckHost(x);
1568:       } else
1569: #endif
1570:       {
1571: #if defined(PETSC_HAVE_CUDA)
1572:         VecCUDAAllocateCheckHost(x);
1573: #endif
1574:       }
1575:     }
1576: #endif
1577:     *a = *((PetscScalar**)x->data);
1578:   } else {
1579:     if (x->ops->getarray) {
1580:       (*x->ops->getarray)(x,a);
1581:     } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array for vector type \"%s\"",((PetscObject)x)->type_name);
1582:   }
1583:   return(0);
1584: }

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

1591:    Logically Collective on Vec

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

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

1599:    Level: beginner

1601: .seealso: VecRestoreArrayInPlace(), VecRestoreArrayInPlace(), VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(),
1602:           VecPlaceArray(), VecGetArray2d(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1603: @*/
1604: PetscErrorCode VecGetArrayInPlace(Vec x,PetscScalar **a)
1605: {

1610:   VecSetErrorIfLocked(x,1);

1612: #if defined(PETSC_HAVE_CUDA)
1613:   if (x->petscnative && (x->offloadmask & PETSC_OFFLOAD_GPU)) { /* Prefer working on GPU when offloadmask is PETSC_OFFLOAD_BOTH */
1614:     PetscBool is_cudatype = PETSC_FALSE;
1615:     PetscObjectTypeCompareAny((PetscObject)x,&is_cudatype,VECSEQCUDA,VECMPICUDA,VECCUDA,"");
1616:     if (is_cudatype) {
1617:       VecCUDAGetArray(x,a);
1618:       x->offloadmask = PETSC_OFFLOAD_GPU; /* Change the mask once GPU gets write access, don't wait until restore array */
1619:       return(0);
1620:     }
1621:   }
1622: #endif
1623:   VecGetArray(x,a);
1624:   return(0);
1625: }

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

1632:    Logically Collective on Vec

1634:    Input Parameter:
1635: .  x - the vector

1637:    Output Parameter:
1638: .  a - location to put pointer to the array

1640:    Level: intermediate

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

1645:    Concepts: vector^accessing local values

1647: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1648:           VecGetArrayPair(), VecRestoreArrayPair(), VecGetArray(), VecRestoreArrayWrite()
1649: @*/
1650: PetscErrorCode VecGetArrayWrite(Vec x,PetscScalar **a)
1651: {

1656:   VecSetErrorIfLocked(x,1);
1657:   if (!x->ops->getarraywrite) {
1658:     VecGetArray(x,a);
1659:   } else {
1660:     (*x->ops->getarraywrite)(x,a);
1661:   }
1662:   return(0);
1663: }

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

1668:    Not Collective

1670:    Input Parameters:
1671: .  x - the vector

1673:    Output Parameter:
1674: .  a - the array

1676:    Level: beginner

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

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

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

1687: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1688: @*/
1689: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1690: {
1692: #if defined(PETSC_HAVE_VIENNACL)
1693:   PetscBool      is_viennacltype = PETSC_FALSE;
1694: #endif

1698:   if (x->petscnative) {
1699: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1700:     if (x->offloadmask == PETSC_OFFLOAD_GPU) {
1701: #if defined(PETSC_HAVE_VIENNACL)
1702:       PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1703:       if (is_viennacltype) {
1704:         VecViennaCLCopyFromGPU(x);
1705:       } else
1706: #endif
1707:       {
1708: #if defined(PETSC_HAVE_CUDA)
1709:         VecCUDACopyFromGPU(x);
1710: #endif
1711:       }
1712:     }
1713: #endif
1714:     *a = *((PetscScalar **)x->data);
1715:   } else if (x->ops->getarrayread) {
1716:     (*x->ops->getarrayread)(x,a);
1717:   } else {
1718:     (*x->ops->getarray)(x,(PetscScalar**)a);
1719:   }
1720:   return(0);
1721: }

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

1728:    Not Collective

1730:    Input Parameters:
1731: .  x - the vector

1733:    Output Parameter:
1734: .  a - the array

1736:    Level: beginner

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


1742: .seealso: VecRestoreArrayReadInPlace(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayInPlace()
1743: @*/
1744: PetscErrorCode VecGetArrayReadInPlace(Vec x,const PetscScalar **a)
1745: {

1750: #if defined(PETSC_HAVE_CUDA)
1751:   if (x->petscnative && x->offloadmask & PETSC_OFFLOAD_GPU) {
1752:     PetscBool is_cudatype = PETSC_FALSE;
1753:     PetscObjectTypeCompareAny((PetscObject)x,&is_cudatype,VECSEQCUDA,VECMPICUDA,VECCUDA,"");
1754:     if (is_cudatype) {
1755:       VecCUDAGetArrayRead(x,a);
1756:       return(0);
1757:     }
1758:   }
1759: #endif
1760:   VecGetArrayRead(x,a);
1761:   return(0);
1762: }

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

1769:    Logically Collective on Vec

1771:    Input Parameter:
1772: +  x - the vectors
1773: -  n - the number of vectors

1775:    Output Parameter:
1776: .  a - location to put pointer to the array

1778:    Fortran Note:
1779:    This routine is not supported in Fortran.

1781:    Level: intermediate

1783: .seealso: VecGetArray(), VecRestoreArrays()
1784: @*/
1785: PetscErrorCode  VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1786: {
1788:   PetscInt       i;
1789:   PetscScalar    **q;

1795:   if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1796:   PetscMalloc1(n,&q);
1797:   for (i=0; i<n; ++i) {
1798:     VecGetArray(x[i],&q[i]);
1799:   }
1800:   *a = q;
1801:   return(0);
1802: }

1804: /*@C
1805:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1806:    has been called.

1808:    Logically Collective on Vec

1810:    Input Parameters:
1811: +  x - the vector
1812: .  n - the number of vectors
1813: -  a - location of pointer to arrays obtained from VecGetArrays()

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

1821:    Fortran Note:
1822:    This routine is not supported in Fortran.

1824:    Level: intermediate

1826: .seealso: VecGetArrays(), VecRestoreArray()
1827: @*/
1828: PetscErrorCode  VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1829: {
1831:   PetscInt       i;
1832:   PetscScalar    **q = *a;


1839:   for (i=0; i<n; ++i) {
1840:     VecRestoreArray(x[i],&q[i]);
1841:   }
1842:   PetscFree(q);
1843:   return(0);
1844: }

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

1849:    Logically Collective on Vec

1851:    Input Parameters:
1852: +  x - the vector
1853: -  a - location of pointer to array obtained from VecGetArray()

1855:    Level: beginner

1857: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1858:           VecGetArrayPair(), VecRestoreArrayPair()
1859: @*/
1860: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1861: {

1866:   if (x->petscnative) {
1867: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1868:     x->offloadmask = PETSC_OFFLOAD_CPU;
1869: #endif
1870:   } else {
1871:     (*x->ops->restorearray)(x,a);
1872:   }
1873:   if (a) *a = NULL;
1874:   PetscObjectStateIncrease((PetscObject)x);
1875:   return(0);
1876: }

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

1881:    Logically Collective on Vec

1883:    Input Parameters:
1884: +  x - the vector
1885: -  a - location of pointer to array obtained from VecGetArrayInPlace()

1887:    Level: beginner

1889: .seealso: VecGetArrayInPlace(), VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(),
1890:           VecPlaceArray(), VecRestoreArray2d(), VecGetArrayPair(), VecRestoreArrayPair()
1891: @*/
1892: PetscErrorCode VecRestoreArrayInPlace(Vec x,PetscScalar **a)
1893: {

1898: #if defined(PETSC_HAVE_CUDA)
1899:   if (x->petscnative && x->offloadmask == PETSC_OFFLOAD_GPU) {
1900:     PetscBool is_cudatype = PETSC_FALSE;
1901:     PetscObjectTypeCompareAny((PetscObject)x,&is_cudatype,VECSEQCUDA,VECMPICUDA,VECCUDA,"");
1902:     if (is_cudatype) {
1903:       VecCUDARestoreArray(x,a);
1904:       return(0);
1905:     }
1906:   }
1907: #endif
1908:   VecRestoreArray(x,a);
1909:   return(0);
1910: }


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

1916:    Logically Collective on Vec

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

1922:    Level: beginner

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

1933:   if (x->petscnative) {
1934: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1935:     x->offloadmask = PETSC_OFFLOAD_CPU;
1936: #endif
1937:   } else {
1938:     if (x->ops->restorearraywrite) {
1939:       (*x->ops->restorearraywrite)(x,a);
1940:     } else {
1941:       (*x->ops->restorearray)(x,a);
1942:     }
1943:   }
1944:   if (a) *a = NULL;
1945:   PetscObjectStateIncrease((PetscObject)x);
1946:   return(0);
1947: }

1949: /*@C
1950:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

1952:    Not Collective

1954:    Input Parameters:
1955: +  vec - the vector
1956: -  array - the array

1958:    Level: beginner

1960: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1961: @*/
1962: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1963: {

1968:   if (x->petscnative) {
1969:     /* nothing */
1970:   } else if (x->ops->restorearrayread) {
1971:     (*x->ops->restorearrayread)(x,a);
1972:   } else {
1973:     (*x->ops->restorearray)(x,(PetscScalar**)a);
1974:   }
1975:   if (a) *a = NULL;
1976:   return(0);
1977: }

1979: /*@C
1980:    VecRestoreArrayReadInPlace - Restore array obtained with VecGetArrayReadInPlace()

1982:    Not Collective

1984:    Input Parameters:
1985: +  vec - the vector
1986: -  array - the array

1988:    Level: beginner

1990: .seealso: VecGetArrayReadInPlace(), VecRestoreArrayInPlace(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1991: @*/
1992: PetscErrorCode VecRestoreArrayReadInPlace(Vec x,const PetscScalar **a)
1993: {

1997:   VecRestoreArrayRead(x,a);
1998:   return(0);
1999: }

2001: /*@
2002:    VecPlaceArray - Allows one to replace the array in a vector with an
2003:    array provided by the user. This is useful to avoid copying an array
2004:    into a vector.

2006:    Not Collective

2008:    Input Parameters:
2009: +  vec - the vector
2010: -  array - the array

2012:    Notes:
2013:    You can return to the original array with a call to VecResetArray()

2015:    Level: developer

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

2019: @*/
2020: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
2021: {

2028:   if (vec->ops->placearray) {
2029:     (*vec->ops->placearray)(vec,array);
2030:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
2031:   PetscObjectStateIncrease((PetscObject)vec);
2032:   return(0);
2033: }

2035: /*@C
2036:    VecReplaceArray - Allows one to replace the array in a vector with an
2037:    array provided by the user. This is useful to avoid copying an array
2038:    into a vector.

2040:    Not Collective

2042:    Input Parameters:
2043: +  vec - the vector
2044: -  array - the array

2046:    Notes:
2047:    This permanently replaces the array and frees the memory associated
2048:    with the old array.

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

2053:    Not supported from Fortran

2055:    Level: developer

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

2059: @*/
2060: PetscErrorCode  VecReplaceArray(Vec vec,const PetscScalar array[])
2061: {

2067:   if (vec->ops->replacearray) {
2068:     (*vec->ops->replacearray)(vec,array);
2069:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
2070:   PetscObjectStateIncrease((PetscObject)vec);
2071:   return(0);
2072: }


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

2078:    This function has semantics similar to VecGetArray():  the pointer
2079:    returned by this function points to a consistent view of the vector
2080:    data.  This may involve a copy operation of data from the host to the
2081:    device if the data on the device is out of date.  If the device
2082:    memory hasn't been allocated previously it will be allocated as part
2083:    of this function call.  VecCUDAGetArray() assumes that
2084:    the user will modify the vector data.  This is similar to
2085:    intent(inout) in fortran.

2087:    The CUDA device pointer has to be released by calling
2088:    VecCUDARestoreArray().  Upon restoring the vector data
2089:    the data on the host will be marked as out of date.  A subsequent
2090:    access of the host data will thus incur a data transfer from the
2091:    device to the host.


2094:    Input Parameter:
2095: .  v - the vector

2097:    Output Parameter:
2098: .  a - the CUDA device pointer

2100:    Fortran note:
2101:    This function is not currently available from Fortran.

2103:    Level: intermediate

2105: .seealso: VecCUDARestoreArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2106: @*/
2107: PETSC_EXTERN PetscErrorCode VecCUDAGetArray(Vec v, PetscScalar **a)
2108: {
2109: #if defined(PETSC_HAVE_CUDA)
2111: #endif

2115: #if defined(PETSC_HAVE_CUDA)
2116:   *a   = 0;
2117:   VecCUDACopyToGPU(v);
2118:   *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2119: #endif
2120:   return(0);
2121: }

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

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

2130:    Input Parameter:
2131: +  v - the vector
2132: -  a - the CUDA device pointer.  This pointer is invalid after
2133:        VecCUDARestoreArray() returns.

2135:    Fortran note:
2136:    This function is not currently available from Fortran.

2138:    Level: intermediate

2140: .seealso: VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2141: @*/
2142: PETSC_EXTERN PetscErrorCode VecCUDARestoreArray(Vec v, PetscScalar **a)
2143: {

2148: #if defined(PETSC_HAVE_CUDA)
2149:   v->offloadmask = PETSC_OFFLOAD_GPU;
2150: #endif

2152:   PetscObjectStateIncrease((PetscObject)v);
2153:   return(0);
2154: }

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

2159:    This function is analogous to VecGetArrayRead():  The pointer
2160:    returned by this function points to a consistent view of the vector
2161:    data.  This may involve a copy operation of data from the host to the
2162:    device if the data on the device is out of date.  If the device
2163:    memory hasn't been allocated previously it will be allocated as part
2164:    of this function call.  VecCUDAGetArrayRead() assumes that the
2165:    user will not modify the vector data.  This is analgogous to
2166:    intent(in) in Fortran.

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

2174:    Input Parameter:
2175: .  v - the vector

2177:    Output Parameter:
2178: .  a - the CUDA pointer.

2180:    Fortran note:
2181:    This function is not currently available from Fortran.

2183:    Level: intermediate

2185: .seealso: VecCUDARestoreArrayRead(), VecCUDAGetArray(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2186: @*/
2187: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayRead(Vec v, const PetscScalar **a)
2188: {
2189: #if defined(PETSC_HAVE_CUDA)
2191: #endif

2195: #if defined(PETSC_HAVE_CUDA)
2196:   *a   = 0;
2197:   VecCUDACopyToGPU(v);
2198:   *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2199: #endif
2200:   return(0);
2201: }

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

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

2211:    Input Parameter:
2212: +  v - the vector
2213: -  a - the CUDA device pointer.  This pointer is invalid after
2214:        VecCUDARestoreArrayRead() returns.

2216:    Fortran note:
2217:    This function is not currently available from Fortran.

2219:    Level: intermediate

2221: .seealso: VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2222: @*/
2223: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayRead(Vec v, const PetscScalar **a)
2224: {
2227:   *a = NULL;
2228:   return(0);
2229: }

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

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

2240:    The device pointer needs to be released with
2241:    VecCUDARestoreArrayWrite().  When the pointer is released the
2242:    host data of the vector is marked as out of data.  Subsequent access
2243:    of the host data with e.g. VecGetArray() incurs a device to host data
2244:    transfer.


2247:    Input Parameter:
2248: .  v - the vector

2250:    Output Parameter:
2251: .  a - the CUDA pointer

2253:    Fortran note:
2254:    This function is not currently available from Fortran.

2256:    Level: advanced

2258: .seealso: VecCUDARestoreArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2259: @*/
2260: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayWrite(Vec v, PetscScalar **a)
2261: {
2262: #if defined(PETSC_HAVE_CUDA)
2264: #endif

2268: #if defined(PETSC_HAVE_CUDA)
2269:   *a   = 0;
2270:   VecCUDAAllocateCheck(v);
2271:   *a   = ((Vec_CUDA*)v->spptr)->GPUarray;
2272: #endif
2273:   return(0);
2274: }

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

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

2283:    Input Parameter:
2284: +  v - the vector
2285: -  a - the CUDA device pointer.  This pointer is invalid after
2286:        VecCUDARestoreArrayWrite() returns.

2288:    Fortran note:
2289:    This function is not currently available from Fortran.

2291:    Level: intermediate

2293: .seealso: VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2294: @*/
2295: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayWrite(Vec v, PetscScalar **a)
2296: {

2301: #if defined(PETSC_HAVE_CUDA)
2302:   v->offloadmask = PETSC_OFFLOAD_GPU;
2303: #endif

2305:   PetscObjectStateIncrease((PetscObject)v);
2306:   return(0);
2307: }

2309: /*@C
2310:    VecCUDAPlaceArray - Allows one to replace the GPU array in a vector with a
2311:    GPU array provided by the user. This is useful to avoid copying an
2312:    array into a vector.

2314:    Not Collective

2316:    Input Parameters:
2317: +  vec - the vector
2318: -  array - the GPU array

2320:    Notes:
2321:    You can return to the original GPU array with a call to VecCUDAResetArray()
2322:    It is not possible to use VecCUDAPlaceArray() and VecPlaceArray() at the
2323:    same time on the same vector.

2325:    Level: developer

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

2329: @*/
2330: PetscErrorCode VecCUDAPlaceArray(Vec vin,PetscScalar *a)
2331: {

2336: #if defined(PETSC_HAVE_CUDA)
2337:   VecCUDACopyToGPU(vin);
2338:   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()");
2339:   ((Vec_Seq*)vin->data)->unplacedarray  = (PetscScalar *) ((Vec_CUDA*)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2340:   ((Vec_CUDA*)vin->spptr)->GPUarray = a;
2341:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2342: #endif
2343:   PetscObjectStateIncrease((PetscObject)vin);
2344:   return(0);
2345: }

2347: /*@C
2348:    VecCUDAReplaceArray - Allows one to replace the GPU array in a vector
2349:    with a GPU array provided by the user. This is useful to avoid copying
2350:    a GPU array into a vector.

2352:    Not Collective

2354:    Input Parameters:
2355: +  vec - the vector
2356: -  array - the GPU array

2358:    Notes:
2359:    This permanently replaces the GPU array and frees the memory associated
2360:    with the old GPU array.

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

2365:    Not supported from Fortran

2367:    Level: developer

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

2371: @*/
2372: PetscErrorCode VecCUDAReplaceArray(Vec vin,PetscScalar *a)
2373: {
2374: #if defined(PETSC_HAVE_CUDA)
2375:   cudaError_t err;
2376: #endif

2381: #if defined(PETSC_HAVE_CUDA)
2382:   err = cudaFree(((Vec_CUDA*)vin->spptr)->GPUarray);CHKERRCUDA(err);
2383:   ((Vec_CUDA*)vin->spptr)->GPUarray = a;
2384:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2385: #endif
2386:   PetscObjectStateIncrease((PetscObject)vin);
2387:   return(0);
2388: }

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

2394:    Not Collective

2396:    Input Parameters:
2397: .  vec - the vector

2399:    Level: developer

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

2403: @*/
2404: PetscErrorCode VecCUDAResetArray(Vec vin)
2405: {

2410: #if defined(PETSC_HAVE_CUDA)
2411:   VecCUDACopyToGPU(vin);
2412:   ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar *) ((Vec_Seq*)vin->data)->unplacedarray;
2413:   ((Vec_Seq*)vin->data)->unplacedarray = 0;
2414:   vin->offloadmask = PETSC_OFFLOAD_GPU;
2415: #endif
2416:   PetscObjectStateIncrease((PetscObject)vin);
2417:   return(0);
2418: }




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

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

2430:     Collective on Vec

2432:     Input Parameters:
2433: +   x - a vector to mimic
2434: -   n - the number of vectors to obtain

2436:     Output Parameters:
2437: +   y - Fortran90 pointer to the array of vectors
2438: -   ierr - error code

2440:     Example of Usage:
2441: .vb
2442: #include <petsc/finclude/petscvec.h>
2443:     use petscvec

2445:     Vec x
2446:     Vec, pointer :: y(:)
2447:     ....
2448:     call VecDuplicateVecsF90(x,2,y,ierr)
2449:     call VecSet(y(2),alpha,ierr)
2450:     call VecSet(y(2),alpha,ierr)
2451:     ....
2452:     call VecDestroyVecsF90(2,y,ierr)
2453: .ve

2455:     Notes:
2456:     Not yet supported for all F90 compilers

2458:     Use VecDestroyVecsF90() to free the space.

2460:     Level: beginner

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

2464: M*/

2466: /*MC
2467:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2468:     VecGetArrayF90().

2470:     Synopsis:
2471:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2473:     Logically Collective on Vec

2475:     Input Parameters:
2476: +   x - vector
2477: -   xx_v - the Fortran90 pointer to the array

2479:     Output Parameter:
2480: .   ierr - error code

2482:     Example of Usage:
2483: .vb
2484: #include <petsc/finclude/petscvec.h>
2485:     use petscvec

2487:     PetscScalar, pointer :: xx_v(:)
2488:     ....
2489:     call VecGetArrayF90(x,xx_v,ierr)
2490:     xx_v(3) = a
2491:     call VecRestoreArrayF90(x,xx_v,ierr)
2492: .ve

2494:     Level: beginner

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

2498: M*/

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

2503:     Synopsis:
2504:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

2506:     Collective on Vec

2508:     Input Parameters:
2509: +   n - the number of vectors previously obtained
2510: -   x - pointer to array of vector pointers

2512:     Output Parameter:
2513: .   ierr - error code

2515:     Notes:
2516:     Not yet supported for all F90 compilers

2518:     Level: beginner

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

2522: M*/

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

2530:     Synopsis:
2531:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2533:     Logically Collective on Vec

2535:     Input Parameter:
2536: .   x - vector

2538:     Output Parameters:
2539: +   xx_v - the Fortran90 pointer to the array
2540: -   ierr - error code

2542:     Example of Usage:
2543: .vb
2544: #include <petsc/finclude/petscvec.h>
2545:     use petscvec

2547:     PetscScalar, pointer :: xx_v(:)
2548:     ....
2549:     call VecGetArrayF90(x,xx_v,ierr)
2550:     xx_v(3) = a
2551:     call VecRestoreArrayF90(x,xx_v,ierr)
2552: .ve

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

2556:     Level: beginner

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

2560: M*/

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

2568:     Synopsis:
2569:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2571:     Logically Collective on Vec

2573:     Input Parameter:
2574: .   x - vector

2576:     Output Parameters:
2577: +   xx_v - the Fortran90 pointer to the array
2578: -   ierr - error code

2580:     Example of Usage:
2581: .vb
2582: #include <petsc/finclude/petscvec.h>
2583:     use petscvec

2585:     PetscScalar, pointer :: xx_v(:)
2586:     ....
2587:     call VecGetArrayReadF90(x,xx_v,ierr)
2588:     a = xx_v(3)
2589:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2590: .ve

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

2594:     Level: beginner

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

2598: M*/

2600: /*MC
2601:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2602:     VecGetArrayReadF90().

2604:     Synopsis:
2605:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2607:     Logically Collective on Vec

2609:     Input Parameters:
2610: +   x - vector
2611: -   xx_v - the Fortran90 pointer to the array

2613:     Output Parameter:
2614: .   ierr - error code

2616:     Example of Usage:
2617: .vb
2618: #include <petsc/finclude/petscvec.h>
2619:     use petscvec

2621:     PetscScalar, pointer :: xx_v(:)
2622:     ....
2623:     call VecGetArrayReadF90(x,xx_v,ierr)
2624:     a = xx_v(3)
2625:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2626: .ve

2628:     Level: beginner

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

2632: M*/

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

2639:    Logically Collective

2641:    Input Parameter:
2642: +  x - the vector
2643: .  m - first dimension of two dimensional array
2644: .  n - second dimension of two dimensional array
2645: .  mstart - first index you will use in first coordinate direction (often 0)
2646: -  nstart - first index in the second coordinate direction (often 0)

2648:    Output Parameter:
2649: .  a - location to put pointer to the array

2651:    Level: developer

2653:   Notes:
2654:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2655:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2656:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2657:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

2661: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2662:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2663:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2664: @*/
2665: PetscErrorCode  VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2666: {
2668:   PetscInt       i,N;
2669:   PetscScalar    *aa;

2675:   VecGetLocalSize(x,&N);
2676:   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);
2677:   VecGetArray(x,&aa);

2679:   PetscMalloc1(m,a);
2680:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2681:   *a -= mstart;
2682:   return(0);
2683: }

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

2690:    Logically Collective

2692:    Input Parameter:
2693: +  x - the vector
2694: .  m - first dimension of two dimensional array
2695: .  n - second dimension of two dimensional array
2696: .  mstart - first index you will use in first coordinate direction (often 0)
2697: -  nstart - first index in the second coordinate direction (often 0)

2699:    Output Parameter:
2700: .  a - location to put pointer to the array

2702:    Level: developer

2704:   Notes:
2705:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2706:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2707:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2708:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

2714: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2715:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2716:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2717: @*/
2718: PetscErrorCode  VecGetArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2719: {
2721:   PetscInt       i,N;
2722:   PetscScalar    *aa;

2728:   VecGetLocalSize(x,&N);
2729:   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);
2730:   VecGetArrayWrite(x,&aa);

2732:   PetscMalloc1(m,a);
2733:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2734:   *a -= mstart;
2735:   return(0);
2736: }

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

2741:    Logically Collective

2743:    Input Parameters:
2744: +  x - the vector
2745: .  m - first dimension of two dimensional array
2746: .  n - second dimension of the two dimensional array
2747: .  mstart - first index you will use in first coordinate direction (often 0)
2748: .  nstart - first index in the second coordinate direction (often 0)
2749: -  a - location of pointer to array obtained from VecGetArray2d()

2751:    Level: developer

2753:    Notes:
2754:    For regular PETSc vectors this routine does not involve any copies. For
2755:    any special vectors that do not store local vector data in a contiguous
2756:    array, this routine will copy the data back into the underlying
2757:    vector data structure from the array obtained with VecGetArray().

2759:    This routine actually zeros out the a pointer.

2761: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2762:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2763:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2764: @*/
2765: PetscErrorCode  VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2766: {
2768:   void           *dummy;

2774:   dummy = (void*)(*a + mstart);
2775:   PetscFree(dummy);
2776:   VecRestoreArray(x,NULL);
2777:   return(0);
2778: }

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

2783:    Logically Collective

2785:    Input Parameters:
2786: +  x - the vector
2787: .  m - first dimension of two dimensional array
2788: .  n - second dimension of the two dimensional array
2789: .  mstart - first index you will use in first coordinate direction (often 0)
2790: .  nstart - first index in the second coordinate direction (often 0)
2791: -  a - location of pointer to array obtained from VecGetArray2d()

2793:    Level: developer

2795:    Notes:
2796:    For regular PETSc vectors this routine does not involve any copies. For
2797:    any special vectors that do not store local vector data in a contiguous
2798:    array, this routine will copy the data back into the underlying
2799:    vector data structure from the array obtained with VecGetArray().

2801:    This routine actually zeros out the a pointer.

2803: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2804:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2805:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2806: @*/
2807: PetscErrorCode  VecRestoreArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2808: {
2810:   void           *dummy;

2816:   dummy = (void*)(*a + mstart);
2817:   PetscFree(dummy);
2818:   VecRestoreArrayWrite(x,NULL);
2819:   return(0);
2820: }

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

2827:    Logically Collective

2829:    Input Parameter:
2830: +  x - the vector
2831: .  m - first dimension of two dimensional array
2832: -  mstart - first index you will use in first coordinate direction (often 0)

2834:    Output Parameter:
2835: .  a - location to put pointer to the array

2837:    Level: developer

2839:   Notes:
2840:    For a vector obtained from DMCreateLocalVector() mstart are likely
2841:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2842:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2846: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2847:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2848:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2849: @*/
2850: PetscErrorCode  VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2851: {
2853:   PetscInt       N;

2859:   VecGetLocalSize(x,&N);
2860:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2861:   VecGetArray(x,a);
2862:   *a  -= mstart;
2863:   return(0);
2864: }

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

2871:    Logically Collective

2873:    Input Parameter:
2874: +  x - the vector
2875: .  m - first dimension of two dimensional array
2876: -  mstart - first index you will use in first coordinate direction (often 0)

2878:    Output Parameter:
2879: .  a - location to put pointer to the array

2881:    Level: developer

2883:   Notes:
2884:    For a vector obtained from DMCreateLocalVector() mstart are likely
2885:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2886:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2890: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2891:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2892:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2893: @*/
2894: PetscErrorCode  VecGetArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2895: {
2897:   PetscInt       N;

2903:   VecGetLocalSize(x,&N);
2904:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2905:   VecGetArrayWrite(x,a);
2906:   *a  -= mstart;
2907:   return(0);
2908: }

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

2913:    Logically Collective

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

2921:    Level: developer

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

2929:    This routine actually zeros out the a pointer.

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

2942:   VecRestoreArray(x,NULL);
2943:   return(0);
2944: }

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

2949:    Logically Collective

2951:    Input Parameters:
2952: +  x - the vector
2953: .  m - first dimension of two dimensional array
2954: .  mstart - first index you will use in first coordinate direction (often 0)
2955: -  a - location of pointer to array obtained from VecGetArray21()

2957:    Level: developer

2959:    Notes:
2960:    For regular PETSc vectors this routine does not involve any copies. For
2961:    any special vectors that do not store local vector data in a contiguous
2962:    array, this routine will copy the data back into the underlying
2963:    vector data structure from the array obtained with VecGetArray1d().

2965:    This routine actually zeros out the a pointer.

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

2969: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2970:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2971:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2972: @*/
2973: PetscErrorCode  VecRestoreArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2974: {

2980:   VecRestoreArrayWrite(x,NULL);
2981:   return(0);
2982: }

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

2989:    Logically Collective

2991:    Input Parameter:
2992: +  x - the vector
2993: .  m - first dimension of three dimensional array
2994: .  n - second dimension of three dimensional array
2995: .  p - third dimension of three dimensional array
2996: .  mstart - first index you will use in first coordinate direction (often 0)
2997: .  nstart - first index in the second coordinate direction (often 0)
2998: -  pstart - first index in the third coordinate direction (often 0)

3000:    Output Parameter:
3001: .  a - location to put pointer to the array

3003:    Level: developer

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

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

3013: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3014:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3015:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3016: @*/
3017: PetscErrorCode  VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3018: {
3020:   PetscInt       i,N,j;
3021:   PetscScalar    *aa,**b;

3027:   VecGetLocalSize(x,&N);
3028:   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);
3029:   VecGetArray(x,&aa);

3031:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3032:   b    = (PetscScalar**)((*a) + m);
3033:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3034:   for (i=0; i<m; i++)
3035:     for (j=0; j<n; j++)
3036:       b[i*n+j] = aa + i*n*p + j*p - pstart;

3038:   *a -= mstart;
3039:   return(0);
3040: }

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

3047:    Logically Collective

3049:    Input Parameter:
3050: +  x - the vector
3051: .  m - first dimension of three dimensional array
3052: .  n - second dimension of three dimensional array
3053: .  p - third dimension of three dimensional array
3054: .  mstart - first index you will use in first coordinate direction (often 0)
3055: .  nstart - first index in the second coordinate direction (often 0)
3056: -  pstart - first index in the third coordinate direction (often 0)

3058:    Output Parameter:
3059: .  a - location to put pointer to the array

3061:    Level: developer

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

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

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

3073: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3074:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3075:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3076: @*/
3077: PetscErrorCode  VecGetArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3078: {
3080:   PetscInt       i,N,j;
3081:   PetscScalar    *aa,**b;

3087:   VecGetLocalSize(x,&N);
3088:   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);
3089:   VecGetArrayWrite(x,&aa);

3091:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3092:   b    = (PetscScalar**)((*a) + m);
3093:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3094:   for (i=0; i<m; i++)
3095:     for (j=0; j<n; j++)
3096:       b[i*n+j] = aa + i*n*p + j*p - pstart;

3098:   *a -= mstart;
3099:   return(0);
3100: }

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

3105:    Logically Collective

3107:    Input Parameters:
3108: +  x - the vector
3109: .  m - first dimension of three dimensional array
3110: .  n - second dimension of the three dimensional array
3111: .  p - third dimension of the three dimensional array
3112: .  mstart - first index you will use in first coordinate direction (often 0)
3113: .  nstart - first index in the second coordinate direction (often 0)
3114: .  pstart - first index in the third coordinate direction (often 0)
3115: -  a - location of pointer to array obtained from VecGetArray3d()

3117:    Level: developer

3119:    Notes:
3120:    For regular PETSc vectors this routine does not involve any copies. For
3121:    any special vectors that do not store local vector data in a contiguous
3122:    array, this routine will copy the data back into the underlying
3123:    vector data structure from the array obtained with VecGetArray().

3125:    This routine actually zeros out the a pointer.

3127: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3128:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3129:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3130: @*/
3131: PetscErrorCode  VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3132: {
3134:   void           *dummy;

3140:   dummy = (void*)(*a + mstart);
3141:   PetscFree(dummy);
3142:   VecRestoreArray(x,NULL);
3143:   return(0);
3144: }

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

3149:    Logically Collective

3151:    Input Parameters:
3152: +  x - the vector
3153: .  m - first dimension of three dimensional array
3154: .  n - second dimension of the three dimensional array
3155: .  p - third dimension of the three dimensional array
3156: .  mstart - first index you will use in first coordinate direction (often 0)
3157: .  nstart - first index in the second coordinate direction (often 0)
3158: .  pstart - first index in the third coordinate direction (often 0)
3159: -  a - location of pointer to array obtained from VecGetArray3d()

3161:    Level: developer

3163:    Notes:
3164:    For regular PETSc vectors this routine does not involve any copies. For
3165:    any special vectors that do not store local vector data in a contiguous
3166:    array, this routine will copy the data back into the underlying
3167:    vector data structure from the array obtained with VecGetArray().

3169:    This routine actually zeros out the a pointer.

3171: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3172:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3173:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3174: @*/
3175: PetscErrorCode  VecRestoreArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3176: {
3178:   void           *dummy;

3184:   dummy = (void*)(*a + mstart);
3185:   PetscFree(dummy);
3186:   VecRestoreArrayWrite(x,NULL);
3187:   return(0);
3188: }

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

3195:    Logically Collective

3197:    Input Parameter:
3198: +  x - the vector
3199: .  m - first dimension of four dimensional array
3200: .  n - second dimension of four dimensional array
3201: .  p - third dimension of four dimensional array
3202: .  q - fourth dimension of four dimensional array
3203: .  mstart - first index you will use in first coordinate direction (often 0)
3204: .  nstart - first index in the second coordinate direction (often 0)
3205: .  pstart - first index in the third coordinate direction (often 0)
3206: -  qstart - first index in the fourth coordinate direction (often 0)

3208:    Output Parameter:
3209: .  a - location to put pointer to the array

3211:    Level: beginner

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

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

3221: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3222:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3223:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3224: @*/
3225: PetscErrorCode  VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3226: {
3228:   PetscInt       i,N,j,k;
3229:   PetscScalar    *aa,***b,**c;

3235:   VecGetLocalSize(x,&N);
3236:   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);
3237:   VecGetArray(x,&aa);

3239:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3240:   b    = (PetscScalar***)((*a) + m);
3241:   c    = (PetscScalar**)(b + m*n);
3242:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3243:   for (i=0; i<m; i++)
3244:     for (j=0; j<n; j++)
3245:       b[i*n+j] = c + i*n*p + j*p - pstart;
3246:   for (i=0; i<m; i++)
3247:     for (j=0; j<n; j++)
3248:       for (k=0; k<p; k++)
3249:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3250:   *a -= mstart;
3251:   return(0);
3252: }

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

3259:    Logically Collective

3261:    Input Parameter:
3262: +  x - the vector
3263: .  m - first dimension of four dimensional array
3264: .  n - second dimension of four dimensional array
3265: .  p - third dimension of four dimensional array
3266: .  q - fourth dimension of four dimensional array
3267: .  mstart - first index you will use in first coordinate direction (often 0)
3268: .  nstart - first index in the second coordinate direction (often 0)
3269: .  pstart - first index in the third coordinate direction (often 0)
3270: -  qstart - first index in the fourth coordinate direction (often 0)

3272:    Output Parameter:
3273: .  a - location to put pointer to the array

3275:    Level: beginner

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

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

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

3287: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3288:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3289:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3290: @*/
3291: PetscErrorCode  VecGetArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3292: {
3294:   PetscInt       i,N,j,k;
3295:   PetscScalar    *aa,***b,**c;

3301:   VecGetLocalSize(x,&N);
3302:   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);
3303:   VecGetArrayWrite(x,&aa);

3305:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3306:   b    = (PetscScalar***)((*a) + m);
3307:   c    = (PetscScalar**)(b + m*n);
3308:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3309:   for (i=0; i<m; i++)
3310:     for (j=0; j<n; j++)
3311:       b[i*n+j] = c + i*n*p + j*p - pstart;
3312:   for (i=0; i<m; i++)
3313:     for (j=0; j<n; j++)
3314:       for (k=0; k<p; k++)
3315:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3316:   *a -= mstart;
3317:   return(0);
3318: }

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

3323:    Logically Collective

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

3337:    Level: beginner

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

3345:    This routine actually zeros out the a pointer.

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

3360:   dummy = (void*)(*a + mstart);
3361:   PetscFree(dummy);
3362:   VecRestoreArray(x,NULL);
3363:   return(0);
3364: }

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

3369:    Logically Collective

3371:    Input Parameters:
3372: +  x - the vector
3373: .  m - first dimension of four dimensional array
3374: .  n - second dimension of the four dimensional array
3375: .  p - third dimension of the four dimensional array
3376: .  q - fourth dimension of the four dimensional array
3377: .  mstart - first index you will use in first coordinate direction (often 0)
3378: .  nstart - first index in the second coordinate direction (often 0)
3379: .  pstart - first index in the third coordinate direction (often 0)
3380: .  qstart - first index in the fourth coordinate direction (often 0)
3381: -  a - location of pointer to array obtained from VecGetArray4d()

3383:    Level: beginner

3385:    Notes:
3386:    For regular PETSc vectors this routine does not involve any copies. For
3387:    any special vectors that do not store local vector data in a contiguous
3388:    array, this routine will copy the data back into the underlying
3389:    vector data structure from the array obtained with VecGetArray().

3391:    This routine actually zeros out the a pointer.

3393: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3394:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3395:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3396: @*/
3397: PetscErrorCode  VecRestoreArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3398: {
3400:   void           *dummy;

3406:   dummy = (void*)(*a + mstart);
3407:   PetscFree(dummy);
3408:   VecRestoreArrayWrite(x,NULL);
3409:   return(0);
3410: }

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

3417:    Logically Collective

3419:    Input Parameter:
3420: +  x - the vector
3421: .  m - first dimension of two dimensional array
3422: .  n - second dimension of two dimensional array
3423: .  mstart - first index you will use in first coordinate direction (often 0)
3424: -  nstart - first index in the second coordinate direction (often 0)

3426:    Output Parameter:
3427: .  a - location to put pointer to the array

3429:    Level: developer

3431:   Notes:
3432:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
3433:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3434:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3435:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

3439: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3440:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3441:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3442: @*/
3443: PetscErrorCode  VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3444: {
3445:   PetscErrorCode    ierr;
3446:   PetscInt          i,N;
3447:   const PetscScalar *aa;

3453:   VecGetLocalSize(x,&N);
3454:   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);
3455:   VecGetArrayRead(x,&aa);

3457:   PetscMalloc1(m,a);
3458:   for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
3459:   *a -= mstart;
3460:   return(0);
3461: }

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

3466:    Logically Collective

3468:    Input Parameters:
3469: +  x - the vector
3470: .  m - first dimension of two dimensional array
3471: .  n - second dimension of the two dimensional array
3472: .  mstart - first index you will use in first coordinate direction (often 0)
3473: .  nstart - first index in the second coordinate direction (often 0)
3474: -  a - location of pointer to array obtained from VecGetArray2d()

3476:    Level: developer

3478:    Notes:
3479:    For regular PETSc vectors this routine does not involve any copies. For
3480:    any special vectors that do not store local vector data in a contiguous
3481:    array, this routine will copy the data back into the underlying
3482:    vector data structure from the array obtained with VecGetArray().

3484:    This routine actually zeros out the a pointer.

3486: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3487:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3488:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3489: @*/
3490: PetscErrorCode  VecRestoreArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3491: {
3493:   void           *dummy;

3499:   dummy = (void*)(*a + mstart);
3500:   PetscFree(dummy);
3501:   VecRestoreArrayRead(x,NULL);
3502:   return(0);
3503: }

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

3510:    Logically Collective

3512:    Input Parameter:
3513: +  x - the vector
3514: .  m - first dimension of two dimensional array
3515: -  mstart - first index you will use in first coordinate direction (often 0)

3517:    Output Parameter:
3518: .  a - location to put pointer to the array

3520:    Level: developer

3522:   Notes:
3523:    For a vector obtained from DMCreateLocalVector() mstart are likely
3524:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3525:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

3529: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3530:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3531:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3532: @*/
3533: PetscErrorCode  VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3534: {
3536:   PetscInt       N;

3542:   VecGetLocalSize(x,&N);
3543:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3544:   VecGetArrayRead(x,(const PetscScalar**)a);
3545:   *a  -= mstart;
3546:   return(0);
3547: }

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

3552:    Logically Collective

3554:    Input Parameters:
3555: +  x - the vector
3556: .  m - first dimension of two dimensional array
3557: .  mstart - first index you will use in first coordinate direction (often 0)
3558: -  a - location of pointer to array obtained from VecGetArray21()

3560:    Level: developer

3562:    Notes:
3563:    For regular PETSc vectors this routine does not involve any copies. For
3564:    any special vectors that do not store local vector data in a contiguous
3565:    array, this routine will copy the data back into the underlying
3566:    vector data structure from the array obtained with VecGetArray1dRead().

3568:    This routine actually zeros out the a pointer.

3570: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3571:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3572:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3573: @*/
3574: PetscErrorCode  VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3575: {

3581:   VecRestoreArrayRead(x,NULL);
3582:   return(0);
3583: }


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

3591:    Logically Collective

3593:    Input Parameter:
3594: +  x - the vector
3595: .  m - first dimension of three dimensional array
3596: .  n - second dimension of three dimensional array
3597: .  p - third dimension of three dimensional array
3598: .  mstart - first index you will use in first coordinate direction (often 0)
3599: .  nstart - first index in the second coordinate direction (often 0)
3600: -  pstart - first index in the third coordinate direction (often 0)

3602:    Output Parameter:
3603: .  a - location to put pointer to the array

3605:    Level: developer

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

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

3615: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3616:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3617:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3618: @*/
3619: PetscErrorCode  VecGetArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3620: {
3621:   PetscErrorCode    ierr;
3622:   PetscInt          i,N,j;
3623:   const PetscScalar *aa;
3624:   PetscScalar       **b;

3630:   VecGetLocalSize(x,&N);
3631:   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);
3632:   VecGetArrayRead(x,&aa);

3634:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3635:   b    = (PetscScalar**)((*a) + m);
3636:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3637:   for (i=0; i<m; i++)
3638:     for (j=0; j<n; j++)
3639:       b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;

3641:   *a -= mstart;
3642:   return(0);
3643: }

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

3648:    Logically Collective

3650:    Input Parameters:
3651: +  x - the vector
3652: .  m - first dimension of three dimensional array
3653: .  n - second dimension of the three dimensional array
3654: .  p - third dimension of the three dimensional array
3655: .  mstart - first index you will use in first coordinate direction (often 0)
3656: .  nstart - first index in the second coordinate direction (often 0)
3657: .  pstart - first index in the third coordinate direction (often 0)
3658: -  a - location of pointer to array obtained from VecGetArray3dRead()

3660:    Level: developer

3662:    Notes:
3663:    For regular PETSc vectors this routine does not involve any copies. For
3664:    any special vectors that do not store local vector data in a contiguous
3665:    array, this routine will copy the data back into the underlying
3666:    vector data structure from the array obtained with VecGetArray().

3668:    This routine actually zeros out the a pointer.

3670: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3671:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3672:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3673: @*/
3674: PetscErrorCode  VecRestoreArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3675: {
3677:   void           *dummy;

3683:   dummy = (void*)(*a + mstart);
3684:   PetscFree(dummy);
3685:   VecRestoreArrayRead(x,NULL);
3686:   return(0);
3687: }

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

3694:    Logically Collective

3696:    Input Parameter:
3697: +  x - the vector
3698: .  m - first dimension of four dimensional array
3699: .  n - second dimension of four dimensional array
3700: .  p - third dimension of four dimensional array
3701: .  q - fourth dimension of four dimensional array
3702: .  mstart - first index you will use in first coordinate direction (often 0)
3703: .  nstart - first index in the second coordinate direction (often 0)
3704: .  pstart - first index in the third coordinate direction (often 0)
3705: -  qstart - first index in the fourth coordinate direction (often 0)

3707:    Output Parameter:
3708: .  a - location to put pointer to the array

3710:    Level: beginner

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

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

3720: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3721:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3722:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3723: @*/
3724: PetscErrorCode  VecGetArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3725: {
3726:   PetscErrorCode    ierr;
3727:   PetscInt          i,N,j,k;
3728:   const PetscScalar *aa;
3729:   PetscScalar       ***b,**c;

3735:   VecGetLocalSize(x,&N);
3736:   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);
3737:   VecGetArrayRead(x,&aa);

3739:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3740:   b    = (PetscScalar***)((*a) + m);
3741:   c    = (PetscScalar**)(b + m*n);
3742:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3743:   for (i=0; i<m; i++)
3744:     for (j=0; j<n; j++)
3745:       b[i*n+j] = c + i*n*p + j*p - pstart;
3746:   for (i=0; i<m; i++)
3747:     for (j=0; j<n; j++)
3748:       for (k=0; k<p; k++)
3749:         c[i*n*p+j*p+k] = (PetscScalar*) aa + i*n*p*q + j*p*q + k*q - qstart;
3750:   *a -= mstart;
3751:   return(0);
3752: }

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

3757:    Logically Collective

3759:    Input Parameters:
3760: +  x - the vector
3761: .  m - first dimension of four dimensional array
3762: .  n - second dimension of the four dimensional array
3763: .  p - third dimension of the four dimensional array
3764: .  q - fourth dimension of the four dimensional array
3765: .  mstart - first index you will use in first coordinate direction (often 0)
3766: .  nstart - first index in the second coordinate direction (often 0)
3767: .  pstart - first index in the third coordinate direction (often 0)
3768: .  qstart - first index in the fourth coordinate direction (often 0)
3769: -  a - location of pointer to array obtained from VecGetArray4dRead()

3771:    Level: beginner

3773:    Notes:
3774:    For regular PETSc vectors this routine does not involve any copies. For
3775:    any special vectors that do not store local vector data in a contiguous
3776:    array, this routine will copy the data back into the underlying
3777:    vector data structure from the array obtained with VecGetArray().

3779:    This routine actually zeros out the a pointer.

3781: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3782:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3783:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3784: @*/
3785: PetscErrorCode  VecRestoreArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3786: {
3788:   void           *dummy;

3794:   dummy = (void*)(*a + mstart);
3795:   PetscFree(dummy);
3796:   VecRestoreArrayRead(x,NULL);
3797:   return(0);
3798: }

3800: #if defined(PETSC_USE_DEBUG)

3802: /*@
3803:    VecLockGet  - Gets the current lock status of a vector

3805:    Logically Collective on Vec

3807:    Input Parameter:
3808: .  x - the vector

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

3814:    Level: beginner

3816: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop()
3817: @*/
3818: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
3819: {
3822:   *state = x->lock;
3823:   return(0);
3824: }

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

3829:    Logically Collective on Vec

3831:    Input Parameter:
3832: .  x - the vector

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

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

3840:    Level: beginner

3842: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPop(), VecLockGet()
3843: @*/
3844: PetscErrorCode VecLockReadPush(Vec x)
3845: {
3848:   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");
3849:   x->lock++;
3850:   return(0);
3851: }

3853: /*@
3854:    VecLockReadPop  - Pops a read-only lock from a vector

3856:    Logically Collective on Vec

3858:    Input Parameter:
3859: .  x - the vector

3861:    Level: beginner

3863: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockGet()
3864: @*/
3865: PetscErrorCode VecLockReadPop(Vec x)
3866: {
3869:   x->lock--;
3870:   if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector has been unlocked from read-only access too many times");
3871:   return(0);
3872: }

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

3877:    Logically Collective on Vec

3879:    Input Parameter:
3880: +  x   - the vector
3881: -  flg - PETSC_TRUE to lock the vector for writing; PETSC_FALSE to unlock it.

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


3890:        VecGetArray(x,&xdata); // begin phase
3891:        VecLockWriteSet_Private(v,PETSC_TRUE);

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

3895:        VecRestoreArray(x,&vdata); // end phase
3896:        VecLockWriteSet_Private(v,PETSC_FALSE);

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

3901:    Level: beginner

3903: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop(), VecLockGet()
3904: @*/
3905: PetscErrorCode VecLockWriteSet_Private(Vec x,PetscBool flg)
3906: {
3909:   if (flg) {
3910:     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");
3911:     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");
3912:     else x->lock = -1;
3913:   } else {
3914:     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");
3915:     x->lock = 0;
3916:   }
3917:   return(0);
3918: }

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

3923:    Level: deprecated

3925: .seealso: VecLockReadPush()
3926: @*/
3927: PetscErrorCode VecLockPush(Vec x)
3928: {
3931:   VecLockReadPush(x);
3932:   return(0);
3933: }

3935: /*@
3936:    VecLockPop  - Pops a read-only lock from a vector

3938:    Level: deprecated

3940: .seealso: VecLockReadPop()
3941: @*/
3942: PetscErrorCode VecLockPop(Vec x)
3943: {
3946:   VecLockReadPop(x);
3947:   return(0);
3948: }

3950: #endif