Actual source code: rvector.c

petsc-3.7.4 2016-10-02
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>       /*I  "petscvec.h"   I*/
  7: static PetscInt VecGetSubVectorSavedStateId = -1;

 10:   if ((x)->map->N != (y)->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths %d != %d", (x)->map->N, (y)->map->N); \
 11:   if ((x)->map->n != (y)->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths %d != %d", (x)->map->n, (y)->map->n);

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

 23: #if defined(PETSC_HAVE_CUSP)
 24:   if ((vec->petscnative || vec->ops->getarray) && (vec->valid_GPU_array == PETSC_CUSP_CPU || vec->valid_GPU_array == PETSC_CUSP_BOTH)) {
 25: #elif defined(PETSC_HAVE_VECCUDA)
 26:   if ((vec->petscnative || vec->ops->getarray) && (vec->valid_GPU_array == PETSC_CUDA_CPU || vec->valid_GPU_array == PETSC_CUDA_BOTH)) {
 27: #else
 28:   if (vec->petscnative || vec->ops->getarray) {
 29: #endif
 30:     VecGetLocalSize(vec,&n);
 31:     VecGetArrayRead(vec,&x);
 32:     for (i=0; i<n; i++) {
 33:       if (begin) {
 34:         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);
 35:       } else {
 36:         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);
 37:       }
 38:     }
 39:     VecRestoreArrayRead(vec,&x);
 40:   }
 41:   return(0);
 42: #else
 43:   return 0;
 44: #endif
 45: }

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

 52:    Logically Collective on Vec

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

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

 60:    Level: advanced

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

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


 80:   (*x->ops->maxpointwisedivide)(x,y,max);
 81:   return(0);
 82: }

 86: /*@
 87:    VecDot - Computes the vector dot product.

 89:    Collective on Vec

 91:    Input Parameters:
 92: .  x, y - the vectors

 94:    Output Parameter:
 95: .  val - the dot product

 97:    Performance Issues:
 98: $    per-processor memory bandwidth
 99: $    interprocessor latency
100: $    work load inbalance that causes certain processes to arrive much earlier than others

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

109:    Use VecTDot() for the indefinite form
110: $     val = (x,y) = y^T x,
111:    where y^T denotes the transpose of y.

113:    Level: intermediate

115:    Concepts: inner product
116:    Concepts: vector^inner product

118: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDotRealPart()
119: @*/
120: PetscErrorCode  VecDot(Vec x,Vec y,PetscScalar *val)
121: {


133:   PetscLogEventBarrierBegin(VEC_DotBarrier,x,y,0,0,PetscObjectComm((PetscObject)x));
134:   (*x->ops->dot)(x,y,val);
135:   PetscLogEventBarrierEnd(VEC_DotBarrier,x,y,0,0,PetscObjectComm((PetscObject)x));
136:   return(0);
137: }

141: /*@
142:    VecDotRealPart - Computes the real part of the vector dot product.

144:    Collective on Vec

146:    Input Parameters:
147: .  x, y - the vectors

149:    Output Parameter:
150: .  val - the real part of the dot product;

152:    Performance Issues:
153: $    per-processor memory bandwidth
154: $    interprocessor latency
155: $    work load inbalance that causes certain processes to arrive much earlier than others

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

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

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

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

167:    Level: intermediate

169:    Concepts: inner product
170:    Concepts: vector^inner product

172: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDot(), VecDotNorm2()
173: @*/
174: PetscErrorCode  VecDotRealPart(Vec x,Vec y,PetscReal *val)
175: {
177:   PetscScalar    fdot;

180:   VecDot(x,y,&fdot);
181:   *val = PetscRealPart(fdot);
182:   return(0);
183: }

187: /*@
188:    VecNorm  - Computes the vector norm.

190:    Collective on Vec

192:    Input Parameters:
193: +  x - the vector
194: -  type - one of NORM_1, NORM_2, NORM_INFINITY.  Also available
195:           NORM_1_AND_2, which computes both norms and stores them
196:           in a two element array.

198:    Output Parameter:
199: .  val - the norm

201:    Notes:
202: $     NORM_1 denotes sum_i |x_i|
203: $     NORM_2 denotes sqrt(sum_i (x_i)^2)
204: $     NORM_INFINITY denotes max_i |x_i|

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

211:    Level: intermediate

213:    Performance Issues:
214: $    per-processor memory bandwidth
215: $    interprocessor latency
216: $    work load inbalance that causes certain processes to arrive much earlier than others

218:    Concepts: norm
219:    Concepts: vector^norm

221: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNormAvailable(),
222:           VecNormBegin(), VecNormEnd()

224: @*/
225: PetscErrorCode  VecNorm(Vec x,NormType type,PetscReal *val)
226: {
227:   PetscBool      flg;

234:   if (((PetscObject)x)->precision != sizeof(PetscReal)) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Wrong precision of input argument");

236:   /*
237:    * Cached data?
238:    */
239:   if (type!=NORM_1_AND_2) {
240:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,flg);
241:     if (flg) return(0);
242:   }
243:   PetscLogEventBarrierBegin(VEC_NormBarrier,x,0,0,0,PetscObjectComm((PetscObject)x));
244:   (*x->ops->norm)(x,type,val);
245:   PetscLogEventBarrierEnd(VEC_NormBarrier,x,0,0,0,PetscObjectComm((PetscObject)x));

247:   if (type!=NORM_1_AND_2) {
248:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[type],*val);
249:   }
250:   return(0);
251: }

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

258:    Not Collective

260:    Input Parameters:
261: +  x - the vector
262: -  type - one of NORM_1, NORM_2, NORM_INFINITY.  Also available
263:           NORM_1_AND_2, which computes both norms and stores them
264:           in a two element array.

266:    Output Parameter:
267: +  available - PETSC_TRUE if the val returned is valid
268: -  val - the norm

270:    Notes:
271: $     NORM_1 denotes sum_i |x_i|
272: $     NORM_2 denotes sqrt(sum_i (x_i)^2)
273: $     NORM_INFINITY denotes max_i |x_i|

275:    Level: intermediate

277:    Performance Issues:
278: $    per-processor memory bandwidth
279: $    interprocessor latency
280: $    work load inbalance that causes certain processes to arrive much earlier than others

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

287:    Concepts: norm
288:    Concepts: vector^norm

290: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNorm()
291:           VecNormBegin(), VecNormEnd()

293: @*/
294: PetscErrorCode  VecNormAvailable(Vec x,NormType type,PetscBool  *available,PetscReal *val)
295: {


303:   *available = PETSC_FALSE;
304:   if (type!=NORM_1_AND_2) {
305:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,*available);
306:   }
307:   return(0);
308: }

312: /*@
313:    VecNormalize - Normalizes a vector by 2-norm.

315:    Collective on Vec

317:    Input Parameters:
318: +  x - the vector

320:    Output Parameter:
321: .  x - the normalized vector
322: -  val - the vector norm before normalization

324:    Level: intermediate

326:    Concepts: vector^normalizing
327:    Concepts: normalizing^vector

329: @*/
330: PetscErrorCode  VecNormalize(Vec x,PetscReal *val)
331: {
333:   PetscReal      norm;

338:   PetscLogEventBegin(VEC_Normalize,x,0,0,0);
339:   VecNorm(x,NORM_2,&norm);
340:   if (norm == 0.0) {
341:     PetscInfo(x,"Vector of zero norm can not be normalized; Returning only the zero norm\n");
342:   } else if (norm != 1.0) {
343:     PetscScalar tmp = 1.0/norm;
344:     VecScale(x,tmp);
345:   }
346:   if (val) *val = norm;
347:   PetscLogEventEnd(VEC_Normalize,x,0,0,0);
348:   return(0);
349: }

353: /*@C
354:    VecMax - Determines the maximum vector component and its location.

356:    Collective on Vec

358:    Input Parameter:
359: .  x - the vector

361:    Output Parameters:
362: +  val - the maximum component
363: -  p - the location of val (pass NULL if you don't want this)

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

368:    Returns the smallest index with the maximum value
369:    Level: intermediate

371:    Concepts: maximum^of vector
372:    Concepts: vector^maximum value

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

384:   PetscLogEventBegin(VEC_Max,x,0,0,0);
385:   (*x->ops->max)(x,p,val);
386:   PetscLogEventEnd(VEC_Max,x,0,0,0);
387:   return(0);
388: }

392: /*@
393:    VecMin - Determines the minimum vector component and its location.

395:    Collective on Vec

397:    Input Parameters:
398: .  x - the vector

400:    Output Parameter:
401: +  val - the minimum component
402: -  p - the location of val (pass NULL if you don't want this location)

404:    Level: intermediate

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

409:    This returns the smallest index with the minumum value

411:    Concepts: minimum^of vector
412:    Concepts: vector^minimum entry

414: .seealso: VecMax()
415: @*/
416: PetscErrorCode  VecMin(Vec x,PetscInt *p,PetscReal *val)
417: {

424:   PetscLogEventBegin(VEC_Min,x,0,0,0);
425:   (*x->ops->min)(x,p,val);
426:   PetscLogEventEnd(VEC_Min,x,0,0,0);
427:   return(0);
428: }

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

436:    Collective on Vec

438:    Input Parameters:
439: .  x, y - the vectors

441:    Output Parameter:
442: .  val - the dot product

444:    Notes for Users of Complex Numbers:
445:    For complex vectors, VecTDot() computes the indefinite form
446: $     val = (x,y) = y^T x,
447:    where y^T denotes the transpose of y.

449:    Use VecDot() for the inner product
450: $     val = (x,y) = y^H x,
451:    where y^H denotes the conjugate transpose of y.

453:    Level: intermediate

455:    Concepts: inner product^non-Hermitian
456:    Concepts: vector^inner product
457:    Concepts: non-Hermitian inner product

459: .seealso: VecDot(), VecMTDot()
460: @*/
461: PetscErrorCode  VecTDot(Vec x,Vec y,PetscScalar *val)
462: {


474:   PetscLogEventBegin(VEC_TDot,x,y,0,0);
475:   (*x->ops->tdot)(x,y,val);
476:   PetscLogEventEnd(VEC_TDot,x,y,0,0);
477:   return(0);
478: }

482: /*@
483:    VecScale - Scales a vector.

485:    Not collective on Vec

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

491:    Output Parameter:
492: .  x - the scaled vector

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

498:    Level: intermediate

500:    Concepts: vector^scaling
501:    Concepts: scaling^vector

503: @*/
504: PetscErrorCode  VecScale(Vec x, PetscScalar alpha)
505: {
506:   PetscReal      norms[4] = {0.0,0.0,0.0, 0.0};
507:   PetscBool      flgs[4];
509:   PetscInt       i;

514:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
515:   PetscLogEventBegin(VEC_Scale,x,0,0,0);
516:   if (alpha != (PetscScalar)1.0) {
517:     /* get current stashed norms */
518:     for (i=0; i<4; i++) {
519:       PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
520:     }
521:     (*x->ops->scale)(x,alpha);
522:     PetscObjectStateIncrease((PetscObject)x);
523:     /* put the scaled stashed norms back into the Vec */
524:     for (i=0; i<4; i++) {
525:       if (flgs[i]) {
526:         PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],PetscAbsScalar(alpha)*norms[i]);
527:       }
528:     }
529:   }
530:   PetscLogEventEnd(VEC_Scale,x,0,0,0);
531:   return(0);
532: }

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

539:    Logically Collective on Vec

541:    Input Parameters:
542: +  x  - the vector
543: -  alpha - the scalar

545:    Output Parameter:
546: .  x  - the vector

548:    Note:
549:    For a vector of dimension n, VecSet() computes
550: $     x[i] = alpha, for i=1,...,n,
551:    so that all vector entries then equal the identical
552:    scalar value, alpha.  Use the more general routine
553:    VecSetValues() to set different vector entries.

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

558:    Level: beginner

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

562:    Concepts: vector^setting to constant

564: @*/
565: PetscErrorCode  VecSet(Vec x,PetscScalar alpha)
566: {
567:   PetscReal      val;

573:   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()");
575:   VecLocked(x,1);

577:   PetscLogEventBegin(VEC_Set,x,0,0,0);
578:   (*x->ops->set)(x,alpha);
579:   PetscLogEventEnd(VEC_Set,x,0,0,0);
580:   PetscObjectStateIncrease((PetscObject)x);

582:   /*  norms can be simply set (if |alpha|*N not too large) */
583:   val  = PetscAbsScalar(alpha);
584:   if (x->map->N == 0) {
585:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],0.0l);
586:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],0.0);
587:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],0.0);
588:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],0.0);
589:   } else if (val > PETSC_MAX_REAL/x->map->N) {
590:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
591:   } else {
592:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],x->map->N * val);
593:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
594:     val  = PetscSqrtReal((PetscReal)x->map->N) * val;
595:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],val);
596:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],val);
597:   }
598:   return(0);
599: }


604: /*@
605:    VecAXPY - Computes y = alpha x + y.

607:    Logically Collective on Vec

609:    Input Parameters:
610: +  alpha - the scalar
611: -  x, y  - the vectors

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

616:    Level: intermediate

618:    Notes: x and y MUST be different vectors

620:    Concepts: vector^BLAS
621:    Concepts: BLAS

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

636:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
638:   VecLocked(y,1);

640:   VecLockPush(x);
641:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
642:   (*y->ops->axpy)(y,alpha,x);
643:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
644:   VecLockPop(x);
645:   PetscObjectStateIncrease((PetscObject)y);
646:   return(0);
647: }

651: /*@
652:    VecAXPBY - Computes y = alpha x + beta y.

654:    Logically Collective on Vec

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

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

663:    Level: intermediate

665:    Notes: x and y MUST be different vectors

667:    Concepts: BLAS
668:    Concepts: vector^BLAS

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

683:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");

687:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
688:   (*y->ops->axpby)(y,alpha,beta,x);
689:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
690:   PetscObjectStateIncrease((PetscObject)y);
691:   return(0);
692: }

696: /*@
697:    VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z

699:    Logically Collective on Vec

701:    Input Parameters:
702: +  alpha,beta, gamma - the scalars
703: -  x, y, z  - the vectors

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

708:    Level: intermediate

710:    Notes: x, y and z must be different vectors

712:    Developer Note:   alpha = 1 or gamma = 1 or gamma = 0.0 are handled as special cases

714:    Concepts: BLAS
715:    Concepts: vector^BLAS

717: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
718: @*/
719: PetscErrorCode  VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
720: {

734:   if (x == y || x == z) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
735:   if (y == z) SETERRQ(PetscObjectComm((PetscObject)y),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");

740:   PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
741:   (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
742:   PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
743:   PetscObjectStateIncrease((PetscObject)z);
744:   return(0);
745: }

749: /*@
750:    VecAYPX - Computes y = x + alpha y.

752:    Logically Collective on Vec

754:    Input Parameters:
755: +  alpha - the scalar
756: -  x, y  - the vectors

758:    Output Parameter:
759: .  y - output vector

761:    Level: intermediate

763:    Notes: x and y MUST be different vectors

765:    Concepts: vector^BLAS
766:    Concepts: BLAS

768: .seealso: VecAXPY(), VecWAXPY()
769: @*/
770: PetscErrorCode  VecAYPX(Vec y,PetscScalar alpha,Vec x)
771: {

779:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y must be different vectors");

782:   PetscLogEventBegin(VEC_AYPX,x,y,0,0);
783:    (*y->ops->aypx)(y,alpha,x);
784:   PetscLogEventEnd(VEC_AYPX,x,y,0,0);
785:   PetscObjectStateIncrease((PetscObject)y);
786:   return(0);
787: }


792: /*@
793:    VecWAXPY - Computes w = alpha x + y.

795:    Logically Collective on Vec

797:    Input Parameters:
798: +  alpha - the scalar
799: -  x, y  - the vectors

801:    Output Parameter:
802: .  w - the result

804:    Level: intermediate

806:    Notes: w cannot be either x or y, but x and y can be the same

808:    Concepts: vector^BLAS
809:    Concepts: BLAS

811: .seealso: VecAXPY(), VecAYPX(), VecAXPBY()
812: @*/
813: PetscErrorCode  VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
814: {

828:   if (w == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector y, suggest VecAXPY()");
829:   if (w == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector x, suggest VecAYPX()");

832:   PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
833:    (*w->ops->waxpy)(w,alpha,x,y);
834:   PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
835:   PetscObjectStateIncrease((PetscObject)w);
836:   return(0);
837: }


842: /*@
843:    VecSetValues - Inserts or adds values into certain locations of a vector.

845:    Not Collective

847:    Input Parameters:
848: +  x - vector to insert in
849: .  ni - number of elements to add
850: .  ix - indices where to add
851: .  y - array of values
852: -  iora - either INSERT_VALUES or ADD_VALUES, where
853:    ADD_VALUES adds values to any existing entries, and
854:    INSERT_VALUES replaces existing entries with new values

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

859:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
860:    options cannot be mixed without intervening calls to the assembly
861:    routines.

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

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

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

874:    Level: beginner

876:    Concepts: vector^setting values

878: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
879:            VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
880: @*/
881: PetscErrorCode  VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
882: {

887:   if (!ni) return(0);
891:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
892:   (*x->ops->setvalues)(x,ni,ix,y,iora);
893:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
894:   PetscObjectStateIncrease((PetscObject)x);
895:   return(0);
896: }

900: /*@
901:    VecGetValues - Gets values from certain locations of a vector. Currently
902:           can only get values on the same processor

904:     Not Collective

906:    Input Parameters:
907: +  x - vector to get values from
908: .  ni - number of elements to get
909: -  ix - indices where to get them from (in global 1d numbering)

911:    Output Parameter:
912: .   y - array of values

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

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

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

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

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

927:    Level: beginner

929:    Concepts: vector^getting values

931: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues()
932: @*/
933: PetscErrorCode  VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
934: {

939:   if (!ni) return(0);
943:   (*x->ops->getvalues)(x,ni,ix,y);
944:   return(0);
945: }

949: /*@
950:    VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.

952:    Not Collective

954:    Input Parameters:
955: +  x - vector to insert in
956: .  ni - number of blocks to add
957: .  ix - indices where to add in block count, rather than element count
958: .  y - array of values
959: -  iora - either INSERT_VALUES or ADD_VALUES, where
960:    ADD_VALUES adds values to any existing entries, and
961:    INSERT_VALUES replaces existing entries with new values

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

967:    Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
968:    options cannot be mixed without intervening calls to the assembly
969:    routines.

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

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

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

981:    Level: intermediate

983:    Concepts: vector^setting values blocked

985: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
986:            VecSetValues()
987: @*/
988: PetscErrorCode  VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
989: {

997:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
998:   (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
999:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1000:   PetscObjectStateIncrease((PetscObject)x);
1001:   return(0);
1002: }


1007: /*@
1008:    VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1009:    using a local ordering of the nodes.

1011:    Not Collective

1013:    Input Parameters:
1014: +  x - vector to insert in
1015: .  ni - number of elements to add
1016: .  ix - indices where to add
1017: .  y - array of values
1018: -  iora - either INSERT_VALUES or ADD_VALUES, where
1019:    ADD_VALUES adds values to any existing entries, and
1020:    INSERT_VALUES replaces existing entries with new values

1022:    Level: intermediate

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

1027:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
1028:    options cannot be mixed without intervening calls to the assembly
1029:    routines.

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

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

1036:    Concepts: vector^setting values with local numbering

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

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

1053:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1054:   if (!x->ops->setvalueslocal) {
1055:     if (!x->map->mapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMapping()");
1056:     if (ni > 128) {
1057:       PetscMalloc1(ni,&lix);
1058:     }
1059:     ISLocalToGlobalMappingApply(x->map->mapping,ni,(PetscInt*)ix,lix);
1060:     (*x->ops->setvalues)(x,ni,lix,y,iora);
1061:     if (ni > 128) {
1062:       PetscFree(lix);
1063:     }
1064:   } else {
1065:     (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
1066:   }
1067:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1068:   PetscObjectStateIncrease((PetscObject)x);
1069:   return(0);
1070: }

1074: /*@
1075:    VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1076:    using a local ordering of the nodes.

1078:    Not Collective

1080:    Input Parameters:
1081: +  x - vector to insert in
1082: .  ni - number of blocks to add
1083: .  ix - indices where to add in block count, not element count
1084: .  y - array of values
1085: -  iora - either INSERT_VALUES or ADD_VALUES, where
1086:    ADD_VALUES adds values to any existing entries, and
1087:    INSERT_VALUES replaces existing entries with new values

1089:    Level: intermediate

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

1095:    Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1096:    options cannot be mixed without intervening calls to the assembly
1097:    routines.

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

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


1105:    Concepts: vector^setting values blocked with local numbering

1107: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
1108:            VecSetLocalToGlobalMapping()
1109: @*/
1110: PetscErrorCode  VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1111: {
1113:   PetscInt       lixp[128],*lix = lixp;

1120:   if (ni > 128) {
1121:     PetscMalloc1(ni,&lix);
1122:   }

1124:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1125:   ISLocalToGlobalMappingApplyBlock(x->map->mapping,ni,(PetscInt*)ix,lix);
1126:   (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1127:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1128:   if (ni > 128) {
1129:     PetscFree(lix);
1130:   }
1131:   PetscObjectStateIncrease((PetscObject)x);
1132:   return(0);
1133: }

1137: /*@
1138:    VecMTDot - Computes indefinite vector multiple dot products.
1139:    That is, it does NOT use the complex conjugate.

1141:    Collective on Vec

1143:    Input Parameters:
1144: +  x - one vector
1145: .  nv - number of vectors
1146: -  y - array of vectors.  Note that vectors are pointers

1148:    Output Parameter:
1149: .  val - array of the dot products

1151:    Notes for Users of Complex Numbers:
1152:    For complex vectors, VecMTDot() computes the indefinite form
1153: $      val = (x,y) = y^T x,
1154:    where y^T denotes the transpose of y.

1156:    Use VecMDot() for the inner product
1157: $      val = (x,y) = y^H x,
1158:    where y^H denotes the conjugate transpose of y.

1160:    Level: intermediate

1162:    Concepts: inner product^multiple
1163:    Concepts: vector^multiple inner products

1165: .seealso: VecMDot(), VecTDot()
1166: @*/
1167: PetscErrorCode  VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1168: {


1181:   PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1182:   (*x->ops->mtdot)(x,nv,y,val);
1183:   PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1184:   return(0);
1185: }

1189: /*@
1190:    VecMDot - Computes vector multiple dot products.

1192:    Collective on Vec

1194:    Input Parameters:
1195: +  x - one vector
1196: .  nv - number of vectors
1197: -  y - array of vectors.

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

1202:    Notes for Users of Complex Numbers:
1203:    For complex vectors, VecMDot() computes
1204: $     val = (x,y) = y^H x,
1205:    where y^H denotes the conjugate transpose of y.

1207:    Use VecMTDot() for the indefinite form
1208: $     val = (x,y) = y^T x,
1209:    where y^T denotes the transpose of y.

1211:    Level: intermediate

1213:    Concepts: inner product^multiple
1214:    Concepts: vector^multiple inner products

1216: .seealso: VecMTDot(), VecDot()
1217: @*/
1218: PetscErrorCode  VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1219: {

1224:   if (!nv) return(0);
1225:   if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);

1234:   PetscLogEventBarrierBegin(VEC_MDotBarrier,x,*y,0,0,PetscObjectComm((PetscObject)x));
1235:   (*x->ops->mdot)(x,nv,y,val);
1236:   PetscLogEventBarrierEnd(VEC_MDotBarrier,x,*y,0,0,PetscObjectComm((PetscObject)x));
1237:   return(0);
1238: }

1242: /*@
1243:    VecMAXPY - Computes y = y + sum alpha[j] x[j]

1245:    Logically Collective on Vec

1247:    Input Parameters:
1248: +  nv - number of scalars and x-vectors
1249: .  alpha - array of scalars
1250: .  y - one vector
1251: -  x - array of vectors

1253:    Level: intermediate

1255:    Notes: y cannot be any of the x vectors

1257:    Concepts: BLAS

1259: .seealso: VecAXPY(), VecWAXPY(), VecAYPX()
1260: @*/
1261: PetscErrorCode  VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1262: {
1264:   PetscInt       i;

1268:   if (!nv) return(0);
1269:   if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);

1279:   PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1280:   (*y->ops->maxpy)(y,nv,alpha,x);
1281:   PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1282:   PetscObjectStateIncrease((PetscObject)y);
1283:   return(0);
1284: }

1288: /*@
1289:    VecGetSubVector - Gets a vector representing part of another vector

1291:    Collective on IS (and Vec if nonlocal entries are needed)

1293:    Input Arguments:
1294: + X - vector from which to extract a subvector
1295: - is - index set representing portion of X to extract

1297:    Output Arguments:
1298: . Y - subvector corresponding to is

1300:    Level: advanced

1302:    Notes:
1303:    The subvector Y should be returned with VecRestoreSubVector().

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

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

1319:   if (X->ops->getsubvector) {
1320:     (*X->ops->getsubvector)(X,is,&Z);
1321:   } else {                      /* Default implementation currently does no caching */
1322:     PetscInt  gstart,gend,start;
1323:     PetscBool contiguous,gcontiguous;
1324:     VecGetOwnershipRange(X,&gstart,&gend);
1325:     ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1326:     MPIU_Allreduce(&contiguous,&gcontiguous,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1327:     if (gcontiguous) {          /* We can do a no-copy implementation */
1328:       PetscInt    n,N,bs;
1329:       PetscMPIInt size;
1330:       PetscInt    state;

1332:       ISGetLocalSize(is,&n);
1333:       VecGetBlockSize(X,&bs);
1334:       if (n%bs || bs == 1) bs = -1; /* Do not decide block size if we do not have to */
1335:       MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1336:       VecLockGet(X,&state);
1337:       if (state) {
1338:         const PetscScalar *x;
1339:         VecGetArrayRead(X,&x);
1340:         if (size == 1) {
1341:           VecCreateSeqWithArray(PetscObjectComm((PetscObject)X),bs,n,(PetscScalar*)x+start,&Z);
1342:         } else {
1343:           ISGetSize(is,&N);
1344:           VecCreateMPIWithArray(PetscObjectComm((PetscObject)X),bs,n,N,(PetscScalar*)x+start,&Z);
1345:         }
1346:         VecRestoreArrayRead(X,&x);
1347:         VecLockPush(Z);
1348:       } else {
1349:         PetscScalar *x;
1350:         VecGetArray(X,&x);
1351:         if (size == 1) {
1352:           VecCreateSeqWithArray(PetscObjectComm((PetscObject)X),bs,n,x+start,&Z);
1353:         } else {
1354:           ISGetSize(is,&N);
1355:           VecCreateMPIWithArray(PetscObjectComm((PetscObject)X),bs,n,N,x+start,&Z);
1356:         }
1357:         VecRestoreArray(X,&x);
1358:       }
1359:     } else {                    /* Have to create a scatter and do a copy */
1360:       VecScatter scatter;
1361:       PetscInt   n,N;
1362:       ISGetLocalSize(is,&n);
1363:       ISGetSize(is,&N);
1364:       VecCreate(PetscObjectComm((PetscObject)is),&Z);
1365:       VecSetSizes(Z,n,N);
1366:       VecSetType(Z,((PetscObject)X)->type_name);
1367:       VecScatterCreate(X,is,Z,NULL,&scatter);
1368:       VecScatterBegin(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1369:       VecScatterEnd(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1370:       PetscObjectCompose((PetscObject)Z,"VecGetSubVector_Scatter",(PetscObject)scatter);
1371:       VecScatterDestroy(&scatter);
1372:     }
1373:   }
1374:   /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1375:   if (VecGetSubVectorSavedStateId < 0) {PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);}
1376:   PetscObjectComposedDataSetInt((PetscObject)Z,VecGetSubVectorSavedStateId,1);
1377:   *Y   = Z;
1378:   return(0);
1379: }

1383: /*@
1384:    VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()

1386:    Collective on IS (and Vec if nonlocal entries need to be written)

1388:    Input Arguments:
1389: + X - vector from which subvector was obtained
1390: . is - index set representing the subset of X
1391: - Y - subvector being restored

1393:    Level: advanced

1395: .seealso: VecGetSubVector()
1396: @*/
1397: PetscErrorCode  VecRestoreSubVector(Vec X,IS is,Vec *Y)
1398: {

1406:   if (X->ops->restoresubvector) {
1407:     (*X->ops->restoresubvector)(X,is,Y);
1408:   } else {
1409:     PETSC_UNUSED PetscObjectState dummystate = 0;
1410:     PetscBool valid;
1411:     PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,dummystate,valid);
1412:     if (!valid) {
1413:       VecScatter scatter;

1415:       PetscObjectQuery((PetscObject)*Y,"VecGetSubVector_Scatter",(PetscObject*)&scatter);
1416:       if (scatter) {
1417:         VecScatterBegin(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1418:         VecScatterEnd(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1419:       }
1420:     }
1421:     VecDestroy(Y);
1422:   }
1423:   return(0);
1424: }

1426: /*@
1427:    VecGetLocalVectorRead - Maps the local portion of a vector into a
1428:    vector.  You must call VecRestoreLocalVectorRead() when the local
1429:    vector is no longer needed.

1431:    Not collective.

1433:    Input parameter:
1434: .  v - The vector for which the local vector is desired.

1436:    Output parameter:
1437: .  w - Upon exit this contains the local vector.

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

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

1455: .seealso: VecRestoreLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1456: @*/
1459: PetscErrorCode VecGetLocalVectorRead(Vec v,Vec w)
1460: {
1462:   PetscScalar    *a;
1463:   PetscInt       m1,m2;

1468:   VecGetLocalSize(v,&m1);
1469:   VecGetLocalSize(w,&m2);
1470:   if (m1 != m2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vectors of different local sizes.");
1471:   if (v->ops->getlocalvectorread) {
1472:     (*v->ops->getlocalvectorread)(v,w);
1473:   } else {
1474:     VecGetArrayRead(v,(const PetscScalar**)&a);
1475:     VecPlaceArray(w,a);
1476:   }
1477:   return(0);
1478: }

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

1484:    Not collective.

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

1490:    Level: beginner

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

1504:   if (v->ops->restorelocalvectorread) {
1505:     (*v->ops->restorelocalvectorread)(v,w);
1506:   } else {
1507:     VecGetArrayRead(w,(const PetscScalar**)&a);
1508:     VecRestoreArrayRead(v,(const PetscScalar**)&a);
1509:     VecResetArray(w);
1510:   }
1511:   return(0);
1512: }

1514: /*@
1515:    VecGetLocalVector - Maps the local portion of a vector into a
1516:    vector.

1518:    Collective on v, not collective on w.

1520:    Input parameter:
1521: .  v - The vector for which the local vector is desired.

1523:    Output parameter:
1524: .  w - Upon exit this contains the local vector.

1526:    Level: beginner
1527:    
1528:    Notes:
1529:    This function is similar to VecGetArray() which maps the local
1530:    portion into a raw pointer.  VecGetLocalVector() is usually about as
1531:    efficient as VecGetArray() but in certain circumstances
1532:    VecGetLocalVector() can be much more efficient than VecGetArray().
1533:    This is because the construction of a contiguous array representing
1534:    the vector data required by VecGetArray() can be an expensive
1535:    operation for certain vector types.  For example, for GPU vectors
1536:    VecGetArray() requires that the data between device and host is
1537:    synchronized.

1539: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1540: @*/
1543: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1544: {
1546:   PetscScalar    *a;
1547:   PetscInt       m1,m2;

1552:   VecGetLocalSize(v,&m1);
1553:   VecGetLocalSize(w,&m2);
1554:   if (m1 != m2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vectors of different local sizes.");
1555:   if (v->ops->getlocalvector) {
1556:     (*v->ops->getlocalvector)(v,w);
1557:   } else {
1558:     VecGetArray(v,&a);
1559:     VecPlaceArray(w,a);
1560:   }
1561:   return(0);
1562: }

1564: /*@
1565:    VecRestoreLocalVector - Unmaps the local portion of a vector
1566:    previously mapped into a vector using VecGetLocalVector().

1568:    Logically collective.

1570:    Input parameter:
1571: .  v - The local portion of this vector was previously mapped into w using VecGetLocalVector().
1572: .  w - The vector into which the local portion of v was mapped.

1574:    Level: beginner

1576: .seealso: VecGetLocalVector(), VecGetLocalVectorRead(), VecRestoreLocalVectorRead(), LocalVectorRead(), VecGetArrayRead(), VecGetArray()
1577: @*/
1580: PetscErrorCode VecRestoreLocalVector(Vec v,Vec w)
1581: {
1583:   PetscScalar    *a;

1588:   if (v->ops->restorelocalvector) {
1589:     (*v->ops->restorelocalvector)(v,w);
1590:   } else {
1591:     VecGetArray(w,&a);
1592:     VecRestoreArray(v,&a);
1593:     VecResetArray(w);
1594:   }
1595:   return(0);
1596: }

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

1609:    Logically Collective on Vec

1611:    Input Parameter:
1612: .  x - the vector

1614:    Output Parameter:
1615: .  a - location to put pointer to the array

1617:    Fortran Note:
1618:    This routine is used differently from Fortran 77
1619: $    Vec         x
1620: $    PetscScalar x_array(1)
1621: $    PetscOffset i_x
1622: $    PetscErrorCode ierr
1623: $       call VecGetArray(x,x_array,i_x,ierr)
1624: $
1625: $   Access first local entry in vector with
1626: $      value = x_array(i_x + 1)
1627: $
1628: $      ...... other code
1629: $       call VecRestoreArray(x,x_array,i_x,ierr)
1630:    For Fortran 90 see VecGetArrayF90()

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

1635:    Level: beginner

1637:    Concepts: vector^accessing local values

1639: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d()
1640: @*/
1641: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1642: {

1647:   VecLocked(x,1);
1648:   if (x->petscnative) {
1649: #if defined(PETSC_HAVE_CUSP)
1650:     if (x->valid_GPU_array == PETSC_CUSP_GPU) {
1651:       VecCUSPCopyFromGPU(x);
1652:     } else if (x->valid_GPU_array == PETSC_CUSP_UNALLOCATED) {
1653:       VecCUSPAllocateCheckHost(x);
1654:     }
1655: #elif defined(PETSC_HAVE_VIENNACL)
1656:     if (x->valid_GPU_array == PETSC_VIENNACL_GPU) {
1657:       VecViennaCLCopyFromGPU(x);
1658:     }
1659: #elif defined(PETSC_HAVE_VECCUDA)
1660:     if (x->valid_GPU_array == PETSC_CUDA_GPU) {
1661:       VecCUDACopyFromGPU(x);
1662:     } else if (x->valid_GPU_array == PETSC_CUDA_UNALLOCATED) {
1663:       VecCUDAAllocateCheckHost(x);
1664:     }
1665: #endif
1666:     *a = *((PetscScalar**)x->data);
1667:   } else {
1668:     (*x->ops->getarray)(x,a);
1669:   }
1670:   return(0);
1671: }

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

1678:    Not Collective

1680:    Input Parameters:
1681: .  x - the vector

1683:    Output Parameter:
1684: .  a - the array

1686:    Level: beginner

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

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

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

1697: .seealso: VecGetArray(), VecRestoreArray()
1698: @*/
1699: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1700: {

1705:   if (x->petscnative) {
1706: #if defined(PETSC_HAVE_CUSP)
1707:     if (x->valid_GPU_array == PETSC_CUSP_GPU) {
1708:       VecCUSPCopyFromGPU(x);
1709:     }
1710: #elif defined(PETSC_HAVE_VIENNACL)
1711:     if (x->valid_GPU_array == PETSC_VIENNACL_GPU) {
1712:       VecViennaCLCopyFromGPU(x);
1713:     }
1714: #elif defined(PETSC_HAVE_VECCUDA)
1715:     if (x->valid_GPU_array == PETSC_CUDA_GPU) {
1716:       VecCUDACopyFromGPU(x);
1717:     }
1718: #endif
1719:     *a = *((PetscScalar **)x->data);
1720:   } else if (x->ops->getarrayread) {
1721:     (*x->ops->getarrayread)(x,a);
1722:   } else {
1723:     (*x->ops->getarray)(x,(PetscScalar**)a);
1724:   }
1725:   return(0);
1726: }

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

1735:    Logically Collective on Vec

1737:    Input Parameter:
1738: +  x - the vectors
1739: -  n - the number of vectors

1741:    Output Parameter:
1742: .  a - location to put pointer to the array

1744:    Fortran Note:
1745:    This routine is not supported in Fortran.

1747:    Level: intermediate

1749: .seealso: VecGetArray(), VecRestoreArrays()
1750: @*/
1751: PetscErrorCode  VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1752: {
1754:   PetscInt       i;
1755:   PetscScalar    **q;

1761:   if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1762:   PetscMalloc1(n,&q);
1763:   for (i=0; i<n; ++i) {
1764:     VecGetArray(x[i],&q[i]);
1765:   }
1766:   *a = q;
1767:   return(0);
1768: }

1772: /*@C
1773:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1774:    has been called.

1776:    Logically Collective on Vec

1778:    Input Parameters:
1779: +  x - the vector
1780: .  n - the number of vectors
1781: -  a - location of pointer to arrays obtained from VecGetArrays()

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

1789:    Fortran Note:
1790:    This routine is not supported in Fortran.

1792:    Level: intermediate

1794: .seealso: VecGetArrays(), VecRestoreArray()
1795: @*/
1796: PetscErrorCode  VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1797: {
1799:   PetscInt       i;
1800:   PetscScalar    **q = *a;


1807:   for (i=0; i<n; ++i) {
1808:     VecRestoreArray(x[i],&q[i]);
1809:   }
1810:   PetscFree(q);
1811:   return(0);
1812: }

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

1819:    Logically Collective on Vec

1821:    Input Parameters:
1822: +  x - the vector
1823: -  a - location of pointer to array obtained from VecGetArray()

1825:    Level: beginner

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

1833:    This routine actually zeros out the a pointer. This is to prevent accidental
1834:    us of the array after it has been restored. If you pass null for a it will
1835:    not zero the array pointer a.

1837:    Fortran Note:
1838:    This routine is used differently from Fortran 77
1839: $    Vec         x
1840: $    PetscScalar x_array(1)
1841: $    PetscOffset i_x
1842: $    PetscErrorCode ierr
1843: $       call VecGetArray(x,x_array,i_x,ierr)
1844: $
1845: $   Access first local entry in vector with
1846: $      value = x_array(i_x + 1)
1847: $
1848: $      ...... other code
1849: $       call VecRestoreArray(x,x_array,i_x,ierr)

1851:    See the Fortran chapter of the users manual and
1852:    petsc/src/snes/examples/tutorials/ex5f.F for details.
1853:    For Fortran 90 see VecRestoreArrayF90()

1855: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d()
1856: @*/
1857: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1858: {

1863:   if (x->petscnative) {
1864: #if defined(PETSC_HAVE_CUSP)
1865:     x->valid_GPU_array = PETSC_CUSP_CPU;
1866: #elif defined(PETSC_HAVE_VIENNACL)
1867:     x->valid_GPU_array = PETSC_VIENNACL_CPU;
1868: #elif defined(PETSC_HAVE_VECCUDA)
1869:     x->valid_GPU_array = PETSC_CUDA_CPU;
1870: #endif
1871:   } else {
1872:     (*x->ops->restorearray)(x,a);
1873:   }
1874:   if (a) *a = NULL;
1875:   PetscObjectStateIncrease((PetscObject)x);
1876:   return(0);
1877: }

1881: /*@C
1882:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

1884:    Not Collective

1886:    Input Parameters:
1887: +  vec - the vector
1888: -  array - the array

1890:    Level: beginner

1892: .seealso: VecGetArray(), VecRestoreArray()
1893: @*/
1894: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1895: {

1900:   if (x->petscnative) {
1901: #if defined(PETSC_HAVE_VIENNACL)
1902:     x->valid_GPU_array = PETSC_VIENNACL_CPU;
1903: #endif
1904:   } else if (x->ops->restorearrayread) {
1905:     (*x->ops->restorearrayread)(x,a);
1906:   } else {
1907:     (*x->ops->restorearray)(x,(PetscScalar**)a);
1908:   }
1909:   if (a) *a = NULL;
1910:   return(0);
1911: }

1915: /*@
1916:    VecPlaceArray - Allows one to replace the array in a vector with an
1917:    array provided by the user. This is useful to avoid copying an array
1918:    into a vector.

1920:    Not Collective

1922:    Input Parameters:
1923: +  vec - the vector
1924: -  array - the array

1926:    Notes:
1927:    You can return to the original array with a call to VecResetArray()

1929:    Level: developer

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

1933: @*/
1934: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
1935: {

1942:   if (vec->ops->placearray) {
1943:     (*vec->ops->placearray)(vec,array);
1944:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
1945:   PetscObjectStateIncrease((PetscObject)vec);
1946:   return(0);
1947: }

1951: /*@C
1952:    VecReplaceArray - Allows one to replace the array in a vector with an
1953:    array provided by the user. This is useful to avoid copying an array
1954:    into a vector.

1956:    Not Collective

1958:    Input Parameters:
1959: +  vec - the vector
1960: -  array - the array

1962:    Notes:
1963:    This permanently replaces the array and frees the memory associated
1964:    with the old array.

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

1969:    Not supported from Fortran

1971:    Level: developer

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

1975: @*/
1976: PetscErrorCode  VecReplaceArray(Vec vec,const PetscScalar array[])
1977: {

1983:   if (vec->ops->replacearray) {
1984:     (*vec->ops->replacearray)(vec,array);
1985:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
1986:   PetscObjectStateIncrease((PetscObject)vec);
1987:   return(0);
1988: }

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

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

1997:     Collective on Vec

1999:     Input Parameters:
2000: +   x - a vector to mimic
2001: -   n - the number of vectors to obtain

2003:     Output Parameters:
2004: +   y - Fortran90 pointer to the array of vectors
2005: -   ierr - error code

2007:     Example of Usage:
2008: .vb
2009:     Vec x
2010:     Vec, pointer :: y(:)
2011:     ....
2012:     call VecDuplicateVecsF90(x,2,y,ierr)
2013:     call VecSet(y(2),alpha,ierr)
2014:     call VecSet(y(2),alpha,ierr)
2015:     ....
2016:     call VecDestroyVecsF90(2,y,ierr)
2017: .ve

2019:     Notes:
2020:     Not yet supported for all F90 compilers

2022:     Use VecDestroyVecsF90() to free the space.

2024:     Level: beginner

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

2028: M*/

2030: /*MC
2031:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2032:     VecGetArrayF90().

2034:     Synopsis:
2035:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2037:     Logically Collective on Vec

2039:     Input Parameters:
2040: +   x - vector
2041: -   xx_v - the Fortran90 pointer to the array

2043:     Output Parameter:
2044: .   ierr - error code

2046:     Example of Usage:
2047: .vb
2048:     PetscScalar, pointer :: xx_v(:)
2049:     ....
2050:     call VecGetArrayF90(x,xx_v,ierr)
2051:     xx_v(3) = a
2052:     call VecRestoreArrayF90(x,xx_v,ierr)
2053: .ve

2055:     Level: beginner

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

2059: M*/

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

2064:     Synopsis:
2065:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

2067:     Collective on Vec

2069:     Input Parameters:
2070: +   n - the number of vectors previously obtained
2071: -   x - pointer to array of vector pointers

2073:     Output Parameter:
2074: .   ierr - error code

2076:     Notes:
2077:     Not yet supported for all F90 compilers

2079:     Level: beginner

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

2083: M*/

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

2091:     Synopsis:
2092:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2094:     Logically Collective on Vec

2096:     Input Parameter:
2097: .   x - vector

2099:     Output Parameters:
2100: +   xx_v - the Fortran90 pointer to the array
2101: -   ierr - error code

2103:     Example of Usage:
2104: .vb
2105:     PetscScalar, pointer :: xx_v(:)
2106:     ....
2107:     call VecGetArrayF90(x,xx_v,ierr)
2108:     xx_v(3) = a
2109:     call VecRestoreArrayF90(x,xx_v,ierr)
2110: .ve

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

2114:     Level: beginner

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

2118: M*/

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

2126:     Synopsis:
2127:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2129:     Logically Collective on Vec

2131:     Input Parameter:
2132: .   x - vector

2134:     Output Parameters:
2135: +   xx_v - the Fortran90 pointer to the array
2136: -   ierr - error code

2138:     Example of Usage:
2139: .vb
2140:     PetscScalar, pointer :: xx_v(:)
2141:     ....
2142:     call VecGetArrayReadF90(x,xx_v,ierr)
2143:     a = xx_v(3)
2144:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2145: .ve

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

2149:     Level: beginner

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

2153: M*/

2155: /*MC
2156:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2157:     VecGetArrayReadF90().

2159:     Synopsis:
2160:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2162:     Logically Collective on Vec

2164:     Input Parameters:
2165: +   x - vector
2166: -   xx_v - the Fortran90 pointer to the array

2168:     Output Parameter:
2169: .   ierr - error code

2171:     Example of Usage:
2172: .vb
2173:     PetscScalar, pointer :: xx_v(:)
2174:     ....
2175:     call VecGetArrayReadF90(x,xx_v,ierr)
2176:     a = xx_v(3)
2177:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2178: .ve

2180:     Level: beginner

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

2184: M*/

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

2193:    Logically Collective

2195:    Input Parameter:
2196: +  x - the vector
2197: .  m - first dimension of two dimensional array
2198: .  n - second dimension of two dimensional array
2199: .  mstart - first index you will use in first coordinate direction (often 0)
2200: -  nstart - first index in the second coordinate direction (often 0)

2202:    Output Parameter:
2203: .  a - location to put pointer to the array

2205:    Level: developer

2207:   Notes:
2208:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2209:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2210:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2211:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

2217: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2218:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2219:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2220: @*/
2221: PetscErrorCode  VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2222: {
2224:   PetscInt       i,N;
2225:   PetscScalar    *aa;

2231:   VecGetLocalSize(x,&N);
2232:   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);
2233:   VecGetArray(x,&aa);

2235:   PetscMalloc1(m,a);
2236:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2237:   *a -= mstart;
2238:   return(0);
2239: }

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

2246:    Logically Collective

2248:    Input Parameters:
2249: +  x - the vector
2250: .  m - first dimension of two dimensional array
2251: .  n - second dimension of the two dimensional array
2252: .  mstart - first index you will use in first coordinate direction (often 0)
2253: .  nstart - first index in the second coordinate direction (often 0)
2254: -  a - location of pointer to array obtained from VecGetArray2d()

2256:    Level: developer

2258:    Notes:
2259:    For regular PETSc vectors this routine does not involve any copies. For
2260:    any special vectors that do not store local vector data in a contiguous
2261:    array, this routine will copy the data back into the underlying
2262:    vector data structure from the array obtained with VecGetArray().

2264:    This routine actually zeros out the a pointer.

2266: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2267:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2268:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2269: @*/
2270: PetscErrorCode  VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2271: {
2273:   void           *dummy;

2279:   dummy = (void*)(*a + mstart);
2280:   PetscFree(dummy);
2281:   VecRestoreArray(x,NULL);
2282:   return(0);
2283: }

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

2292:    Logically Collective

2294:    Input Parameter:
2295: +  x - the vector
2296: .  m - first dimension of two dimensional array
2297: -  mstart - first index you will use in first coordinate direction (often 0)

2299:    Output Parameter:
2300: .  a - location to put pointer to the array

2302:    Level: developer

2304:   Notes:
2305:    For a vector obtained from DMCreateLocalVector() mstart are likely
2306:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2307:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2311: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2312:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2313:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2314: @*/
2315: PetscErrorCode  VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2316: {
2318:   PetscInt       N;

2324:   VecGetLocalSize(x,&N);
2325:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2326:   VecGetArray(x,a);
2327:   *a  -= mstart;
2328:   return(0);
2329: }

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

2336:    Logically Collective

2338:    Input Parameters:
2339: +  x - the vector
2340: .  m - first dimension of two dimensional array
2341: .  mstart - first index you will use in first coordinate direction (often 0)
2342: -  a - location of pointer to array obtained from VecGetArray21()

2344:    Level: developer

2346:    Notes:
2347:    For regular PETSc vectors this routine does not involve any copies. For
2348:    any special vectors that do not store local vector data in a contiguous
2349:    array, this routine will copy the data back into the underlying
2350:    vector data structure from the array obtained with VecGetArray1d().

2352:    This routine actually zeros out the a pointer.

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

2356: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2357:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2358:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2359: @*/
2360: PetscErrorCode  VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2361: {

2367:   VecRestoreArray(x,NULL);
2368:   return(0);
2369: }


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

2379:    Logically Collective

2381:    Input Parameter:
2382: +  x - the vector
2383: .  m - first dimension of three dimensional array
2384: .  n - second dimension of three dimensional array
2385: .  p - third dimension of three dimensional array
2386: .  mstart - first index you will use in first coordinate direction (often 0)
2387: .  nstart - first index in the second coordinate direction (often 0)
2388: -  pstart - first index in the third coordinate direction (often 0)

2390:    Output Parameter:
2391: .  a - location to put pointer to the array

2393:    Level: developer

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

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

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

2405: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2406:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2407:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2408: @*/
2409: PetscErrorCode  VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2410: {
2412:   PetscInt       i,N,j;
2413:   PetscScalar    *aa,**b;

2419:   VecGetLocalSize(x,&N);
2420:   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);
2421:   VecGetArray(x,&aa);

2423:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
2424:   b    = (PetscScalar**)((*a) + m);
2425:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2426:   for (i=0; i<m; i++)
2427:     for (j=0; j<n; j++)
2428:       b[i*n+j] = aa + i*n*p + j*p - pstart;

2430:   *a -= mstart;
2431:   return(0);
2432: }

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

2439:    Logically Collective

2441:    Input Parameters:
2442: +  x - the vector
2443: .  m - first dimension of three dimensional array
2444: .  n - second dimension of the three dimensional array
2445: .  p - third dimension of the three dimensional array
2446: .  mstart - first index you will use in first coordinate direction (often 0)
2447: .  nstart - first index in the second coordinate direction (often 0)
2448: .  pstart - first index in the third coordinate direction (often 0)
2449: -  a - location of pointer to array obtained from VecGetArray3d()

2451:    Level: developer

2453:    Notes:
2454:    For regular PETSc vectors this routine does not involve any copies. For
2455:    any special vectors that do not store local vector data in a contiguous
2456:    array, this routine will copy the data back into the underlying
2457:    vector data structure from the array obtained with VecGetArray().

2459:    This routine actually zeros out the a pointer.

2461: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2462:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2463:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2464: @*/
2465: PetscErrorCode  VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2466: {
2468:   void           *dummy;

2474:   dummy = (void*)(*a + mstart);
2475:   PetscFree(dummy);
2476:   VecRestoreArray(x,NULL);
2477:   return(0);
2478: }

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

2487:    Logically Collective

2489:    Input Parameter:
2490: +  x - the vector
2491: .  m - first dimension of four dimensional array
2492: .  n - second dimension of four dimensional array
2493: .  p - third dimension of four dimensional array
2494: .  q - fourth dimension of four dimensional array
2495: .  mstart - first index you will use in first coordinate direction (often 0)
2496: .  nstart - first index in the second coordinate direction (often 0)
2497: .  pstart - first index in the third coordinate direction (often 0)
2498: -  qstart - first index in the fourth coordinate direction (often 0)

2500:    Output Parameter:
2501: .  a - location to put pointer to the array

2503:    Level: beginner

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

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

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

2515: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2516:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2517:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2518: @*/
2519: PetscErrorCode  VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2520: {
2522:   PetscInt       i,N,j,k;
2523:   PetscScalar    *aa,***b,**c;

2529:   VecGetLocalSize(x,&N);
2530:   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);
2531:   VecGetArray(x,&aa);

2533:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
2534:   b    = (PetscScalar***)((*a) + m);
2535:   c    = (PetscScalar**)(b + m*n);
2536:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2537:   for (i=0; i<m; i++)
2538:     for (j=0; j<n; j++)
2539:       b[i*n+j] = c + i*n*p + j*p - pstart;
2540:   for (i=0; i<m; i++)
2541:     for (j=0; j<n; j++)
2542:       for (k=0; k<p; k++)
2543:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
2544:   *a -= mstart;
2545:   return(0);
2546: }

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

2553:    Logically Collective

2555:    Input Parameters:
2556: +  x - the vector
2557: .  m - first dimension of four dimensional array
2558: .  n - second dimension of the four dimensional array
2559: .  p - third dimension of the four dimensional array
2560: .  q - fourth dimension of the four dimensional array
2561: .  mstart - first index you will use in first coordinate direction (often 0)
2562: .  nstart - first index in the second coordinate direction (often 0)
2563: .  pstart - first index in the third coordinate direction (often 0)
2564: .  qstart - first index in the fourth coordinate direction (often 0)
2565: -  a - location of pointer to array obtained from VecGetArray4d()

2567:    Level: beginner

2569:    Notes:
2570:    For regular PETSc vectors this routine does not involve any copies. For
2571:    any special vectors that do not store local vector data in a contiguous
2572:    array, this routine will copy the data back into the underlying
2573:    vector data structure from the array obtained with VecGetArray().

2575:    This routine actually zeros out the a pointer.

2577: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2578:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2579:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2580: @*/
2581: PetscErrorCode  VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2582: {
2584:   void           *dummy;

2590:   dummy = (void*)(*a + mstart);
2591:   PetscFree(dummy);
2592:   VecRestoreArray(x,NULL);
2593:   return(0);
2594: }

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

2603:    Logically Collective

2605:    Input Parameter:
2606: +  x - the vector
2607: .  m - first dimension of two dimensional array
2608: .  n - second dimension of two dimensional array
2609: .  mstart - first index you will use in first coordinate direction (often 0)
2610: -  nstart - first index in the second coordinate direction (often 0)

2612:    Output Parameter:
2613: .  a - location to put pointer to the array

2615:    Level: developer

2617:   Notes:
2618:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2619:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2620:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2621:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

2627: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2628:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2629:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2630: @*/
2631: PetscErrorCode  VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2632: {
2633:   PetscErrorCode    ierr;
2634:   PetscInt          i,N;
2635:   const PetscScalar *aa;

2641:   VecGetLocalSize(x,&N);
2642:   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);
2643:   VecGetArrayRead(x,&aa);

2645:   PetscMalloc1(m,a);
2646:   for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
2647:   *a -= mstart;
2648:   return(0);
2649: }

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

2656:    Logically Collective

2658:    Input Parameters:
2659: +  x - the vector
2660: .  m - first dimension of two dimensional array
2661: .  n - second dimension of the two dimensional array
2662: .  mstart - first index you will use in first coordinate direction (often 0)
2663: .  nstart - first index in the second coordinate direction (often 0)
2664: -  a - location of pointer to array obtained from VecGetArray2d()

2666:    Level: developer

2668:    Notes:
2669:    For regular PETSc vectors this routine does not involve any copies. For
2670:    any special vectors that do not store local vector data in a contiguous
2671:    array, this routine will copy the data back into the underlying
2672:    vector data structure from the array obtained with VecGetArray().

2674:    This routine actually zeros out the a pointer.

2676: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2677:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2678:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2679: @*/
2680: PetscErrorCode  VecRestoreArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2681: {
2683:   void           *dummy;

2689:   dummy = (void*)(*a + mstart);
2690:   PetscFree(dummy);
2691:   VecRestoreArrayRead(x,NULL);
2692:   return(0);
2693: }

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

2702:    Logically Collective

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

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

2712:    Level: developer

2714:   Notes:
2715:    For a vector obtained from DMCreateLocalVector() mstart are likely
2716:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2717:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2721: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2722:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2723:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2724: @*/
2725: PetscErrorCode  VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2726: {
2728:   PetscInt       N;

2734:   VecGetLocalSize(x,&N);
2735:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2736:   VecGetArrayRead(x,(const PetscScalar**)a);
2737:   *a  -= mstart;
2738:   return(0);
2739: }

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

2746:    Logically Collective

2748:    Input Parameters:
2749: +  x - the vector
2750: .  m - first dimension of two dimensional array
2751: .  mstart - first index you will use in first coordinate direction (often 0)
2752: -  a - location of pointer to array obtained from VecGetArray21()

2754:    Level: developer

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

2762:    This routine actually zeros out the a pointer.

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

2766: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2767:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2768:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2769: @*/
2770: PetscErrorCode  VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2771: {

2777:   VecRestoreArrayRead(x,NULL);
2778:   return(0);
2779: }


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

2789:    Logically Collective

2791:    Input Parameter:
2792: +  x - the vector
2793: .  m - first dimension of three dimensional array
2794: .  n - second dimension of three dimensional array
2795: .  p - third dimension of three dimensional array
2796: .  mstart - first index you will use in first coordinate direction (often 0)
2797: .  nstart - first index in the second coordinate direction (often 0)
2798: -  pstart - first index in the third coordinate direction (often 0)

2800:    Output Parameter:
2801: .  a - location to put pointer to the array

2803:    Level: developer

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

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

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

2815: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2816:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2817:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2818: @*/
2819: PetscErrorCode  VecGetArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2820: {
2821:   PetscErrorCode    ierr;
2822:   PetscInt          i,N,j;
2823:   const PetscScalar *aa;
2824:   PetscScalar       **b;

2830:   VecGetLocalSize(x,&N);
2831:   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);
2832:   VecGetArrayRead(x,&aa);

2834:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
2835:   b    = (PetscScalar**)((*a) + m);
2836:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2837:   for (i=0; i<m; i++)
2838:     for (j=0; j<n; j++)
2839:       b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;

2841:   *a -= mstart;
2842:   return(0);
2843: }

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

2850:    Logically Collective

2852:    Input Parameters:
2853: +  x - the vector
2854: .  m - first dimension of three dimensional array
2855: .  n - second dimension of the three dimensional array
2856: .  p - third dimension of the three dimensional array
2857: .  mstart - first index you will use in first coordinate direction (often 0)
2858: .  nstart - first index in the second coordinate direction (often 0)
2859: .  pstart - first index in the third coordinate direction (often 0)
2860: -  a - location of pointer to array obtained from VecGetArray3dRead()

2862:    Level: developer

2864:    Notes:
2865:    For regular PETSc vectors this routine does not involve any copies. For
2866:    any special vectors that do not store local vector data in a contiguous
2867:    array, this routine will copy the data back into the underlying
2868:    vector data structure from the array obtained with VecGetArray().

2870:    This routine actually zeros out the a pointer.

2872: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2873:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2874:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2875: @*/
2876: PetscErrorCode  VecRestoreArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2877: {
2879:   void           *dummy;

2885:   dummy = (void*)(*a + mstart);
2886:   PetscFree(dummy);
2887:   VecRestoreArrayRead(x,NULL);
2888:   return(0);
2889: }

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

2898:    Logically Collective

2900:    Input Parameter:
2901: +  x - the vector
2902: .  m - first dimension of four dimensional array
2903: .  n - second dimension of four dimensional array
2904: .  p - third dimension of four dimensional array
2905: .  q - fourth dimension of four dimensional array
2906: .  mstart - first index you will use in first coordinate direction (often 0)
2907: .  nstart - first index in the second coordinate direction (often 0)
2908: .  pstart - first index in the third coordinate direction (often 0)
2909: -  qstart - first index in the fourth coordinate direction (often 0)

2911:    Output Parameter:
2912: .  a - location to put pointer to the array

2914:    Level: beginner

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

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

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

2926: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2927:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2928:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2929: @*/
2930: PetscErrorCode  VecGetArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2931: {
2932:   PetscErrorCode    ierr;
2933:   PetscInt          i,N,j,k;
2934:   const PetscScalar *aa;
2935:   PetscScalar       ***b,**c;

2941:   VecGetLocalSize(x,&N);
2942:   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);
2943:   VecGetArrayRead(x,&aa);

2945:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
2946:   b    = (PetscScalar***)((*a) + m);
2947:   c    = (PetscScalar**)(b + m*n);
2948:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2949:   for (i=0; i<m; i++)
2950:     for (j=0; j<n; j++)
2951:       b[i*n+j] = c + i*n*p + j*p - pstart;
2952:   for (i=0; i<m; i++)
2953:     for (j=0; j<n; j++)
2954:       for (k=0; k<p; k++)
2955:         c[i*n*p+j*p+k] = (PetscScalar*) aa + i*n*p*q + j*p*q + k*q - qstart;
2956:   *a -= mstart;
2957:   return(0);
2958: }

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

2965:    Logically Collective

2967:    Input Parameters:
2968: +  x - the vector
2969: .  m - first dimension of four dimensional array
2970: .  n - second dimension of the four dimensional array
2971: .  p - third dimension of the four dimensional array
2972: .  q - fourth dimension of the four dimensional array
2973: .  mstart - first index you will use in first coordinate direction (often 0)
2974: .  nstart - first index in the second coordinate direction (often 0)
2975: .  pstart - first index in the third coordinate direction (often 0)
2976: .  qstart - first index in the fourth coordinate direction (often 0)
2977: -  a - location of pointer to array obtained from VecGetArray4dRead()

2979:    Level: beginner

2981:    Notes:
2982:    For regular PETSc vectors this routine does not involve any copies. For
2983:    any special vectors that do not store local vector data in a contiguous
2984:    array, this routine will copy the data back into the underlying
2985:    vector data structure from the array obtained with VecGetArray().

2987:    This routine actually zeros out the a pointer.

2989: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2990:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2991:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2992: @*/
2993: PetscErrorCode  VecRestoreArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2994: {
2996:   void           *dummy;

3002:   dummy = (void*)(*a + mstart);
3003:   PetscFree(dummy);
3004:   VecRestoreArrayRead(x,NULL);
3005:   return(0);
3006: }

3008: #if defined(PETSC_USE_DEBUG)

3012: /*@
3013:    VecLockGet  - Gets the current lock status of a vector

3015:    Logically Collective on Vec

3017:    Input Parameter:
3018: .  x - the vector

3020:    Output Parameter:
3021: .  state - greater than zero indicates the vector is still locked

3023:    Level: beginner

3025:    Concepts: vector^accessing local values

3027: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockPush(), VecLockGet()
3028: @*/
3029: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
3030: {
3033:   *state = x->lock;
3034:   return(0);
3035: }

3039: /*@
3040:    VecLockPush  - Lock a vector from writing

3042:    Logically Collective on Vec

3044:    Input Parameter:
3045: .  x - the vector

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

3049:     Call VecLockPop() to remove the latest lock

3051:    Level: beginner

3053:    Concepts: vector^accessing local values

3055: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockPop(), VecLockGet()
3056: @*/
3057: PetscErrorCode VecLockPush(Vec x)
3058: {
3061:   x->lock++;
3062:   return(0);
3063: }

3067: /*@
3068:    VecLockPop  - Unlock a vector from writing

3070:    Logically Collective on Vec

3072:    Input Parameter:
3073: .  x - the vector

3075:    Level: beginner

3077:    Concepts: vector^accessing local values

3079: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockPush(), VecLockGet()
3080: @*/
3081: PetscErrorCode VecLockPop(Vec x)
3082: {
3085:   x->lock--;
3086:   if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector has been unlocked too many times");
3087:   return(0);
3088: }

3090: #endif