Actual source code: rvector.c

petsc-3.4.4 2014-03-13
  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 (vec->petscnative || vec->ops->getarray) {
 24:     VecGetLocalSize(vec,&n);
 25:     VecGetArrayRead(vec,&x);
 26:     for (i=0; i<n; i++) {
 27:       if (begin) {
 28:         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);
 29:       } else {
 30:         if (PetscIsInfOrNanScalar(x[i])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FP,"Vec entry at local location %D is not-a-number or infinite at end of function: Parameter number %D",i,argnum);
 31:       }
 32:     }
 33:     VecRestoreArrayRead(vec,&x);
 34:   }
 35:   return(0);
 36: #else
 37:   return 0;
 38: #endif
 39: }

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

 46:    Logically Collective on Vec

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

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

 54:    Level: advanced

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

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


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

 80: /*@
 81:    VecDot - Computes the vector dot product.

 83:    Collective on Vec

 85:    Input Parameters:
 86: .  x, y - the vectors

 88:    Output Parameter:
 89: .  val - the dot product

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

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

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

107:    Level: intermediate

109:    Concepts: inner product
110:    Concepts: vector^inner product

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


127:   PetscLogEventBarrierBegin(VEC_DotBarrier,x,y,0,0,PetscObjectComm((PetscObject)x));
128:   (*x->ops->dot)(x,y,val);
129:   PetscLogEventBarrierEnd(VEC_DotBarrier,x,y,0,0,PetscObjectComm((PetscObject)x));
130:   return(0);
131: }

135: /*@
136:    VecDotRealPart - Computes the real part of the vector dot product.

138:    Collective on Vec

140:    Input Parameters:
141: .  x, y - the vectors

143:    Output Parameter:
144: .  val - the real part of the dot product;

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

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

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

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

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

161:    Level: intermediate

163:    Concepts: inner product
164:    Concepts: vector^inner product

166: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDot(), VecDotNorm2()
167: @*/
168: PetscErrorCode  VecDotRealPart(Vec x,Vec y,PetscReal *val)
169: {
171:   PetscScalar    fdot;

174:   VecDot(x,y,&fdot);
175:   *val = PetscRealPart(fdot);
176:   return(0);
177: }

181: /*@
182:    VecNorm  - Computes the vector norm.

184:    Collective on Vec

186:    Input Parameters:
187: +  x - the vector
188: -  type - one of NORM_1, NORM_2, NORM_INFINITY.  Also available
189:           NORM_1_AND_2, which computes both norms and stores them
190:           in a two element array.

192:    Output Parameter:
193: .  val - the norm

195:    Notes:
196: $     NORM_1 denotes sum_i |x_i|
197: $     NORM_2 denotes sqrt(sum_i (x_i)^2)
198: $     NORM_INFINITY denotes max_i |x_i|

200:    Level: intermediate

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

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

212:    Concepts: norm
213:    Concepts: vector^norm

215: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNormAvailable(),
216:           VecNormBegin(), VecNormEnd()

218: @*/
219: PetscErrorCode  VecNorm(Vec x,NormType type,PetscReal *val)
220: {
221:   PetscBool      flg;

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

230:   /*
231:    * Cached data?
232:    */
233:   if (type!=NORM_1_AND_2) {
234:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,flg);
235:     if (flg) return(0);
236:   }
237:   PetscLogEventBarrierBegin(VEC_NormBarrier,x,0,0,0,PetscObjectComm((PetscObject)x));
238:   (*x->ops->norm)(x,type,val);
239:   PetscLogEventBarrierEnd(VEC_NormBarrier,x,0,0,0,PetscObjectComm((PetscObject)x));

241:   if (type!=NORM_1_AND_2) {
242:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[type],*val);
243:   }
244:   return(0);
245: }

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

252:    Not Collective

254:    Input Parameters:
255: +  x - the vector
256: -  type - one of NORM_1, NORM_2, NORM_INFINITY.  Also available
257:           NORM_1_AND_2, which computes both norms and stores them
258:           in a two element array.

260:    Output Parameter:
261: +  available - PETSC_TRUE if the val returned is valid
262: -  val - the norm

264:    Notes:
265: $     NORM_1 denotes sum_i |x_i|
266: $     NORM_2 denotes sqrt(sum_i (x_i)^2)
267: $     NORM_INFINITY denotes max_i |x_i|

269:    Level: intermediate

271:    Performance Issues:
272: $    per-processor memory bandwidth
273: $    interprocessor latency
274: $    work load inbalance that causes certain processes to arrive much earlier than others

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

281:    Concepts: norm
282:    Concepts: vector^norm

284: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNorm()
285:           VecNormBegin(), VecNormEnd()

287: @*/
288: PetscErrorCode  VecNormAvailable(Vec x,NormType type,PetscBool  *available,PetscReal *val)
289: {


297:   *available = PETSC_FALSE;
298:   if (type!=NORM_1_AND_2) {
299:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,*available);
300:   }
301:   return(0);
302: }

306: /*@
307:    VecNormalize - Normalizes a vector by 2-norm.

309:    Collective on Vec

311:    Input Parameters:
312: +  x - the vector

314:    Output Parameter:
315: .  x - the normalized vector
316: -  val - the vector norm before normalization

318:    Level: intermediate

320:    Concepts: vector^normalizing
321:    Concepts: normalizing^vector

323: @*/
324: PetscErrorCode  VecNormalize(Vec x,PetscReal *val)
325: {
327:   PetscReal      norm;

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

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

350:    Collective on Vec

352:    Input Parameter:
353: .  x - the vector

355:    Output Parameters:
356: +  val - the maximum component
357: -  p - the location of val (pass NULL if you don't want this)

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

362:    Returns the smallest index with the maximum value
363:    Level: intermediate

365:    Concepts: maximum^of vector
366:    Concepts: vector^maximum value

368: .seealso: VecNorm(), VecMin()
369: @*/
370: PetscErrorCode  VecMax(Vec x,PetscInt *p,PetscReal *val)
371: {

378:   PetscLogEventBegin(VEC_Max,x,0,0,0);
379:   (*x->ops->max)(x,p,val);
380:   PetscLogEventEnd(VEC_Max,x,0,0,0);
381:   return(0);
382: }

386: /*@
387:    VecMin - Determines the minimum vector component and its location.

389:    Collective on Vec

391:    Input Parameters:
392: .  x - the vector

394:    Output Parameter:
395: +  val - the minimum component
396: -  p - the location of val (pass NULL if you don't want this location)

398:    Level: intermediate

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

403:    This returns the smallest index with the minumum value

405:    Concepts: minimum^of vector
406:    Concepts: vector^minimum entry

408: .seealso: VecMax()
409: @*/
410: PetscErrorCode  VecMin(Vec x,PetscInt *p,PetscReal *val)
411: {

418:   PetscLogEventBegin(VEC_Min,x,0,0,0);
419:   (*x->ops->min)(x,p,val);
420:   PetscLogEventEnd(VEC_Min,x,0,0,0);
421:   return(0);
422: }

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

430:    Collective on Vec

432:    Input Parameters:
433: .  x, y - the vectors

435:    Output Parameter:
436: .  val - the dot product

438:    Notes for Users of Complex Numbers:
439:    For complex vectors, VecTDot() computes the indefinite form
440: $     val = (x,y) = y^T x,
441:    where y^T denotes the transpose of y.

443:    Use VecDot() for the inner product
444: $     val = (x,y) = y^H x,
445:    where y^H denotes the conjugate transpose of y.

447:    Level: intermediate

449:    Concepts: inner product^non-Hermitian
450:    Concepts: vector^inner product
451:    Concepts: non-Hermitian inner product

453: .seealso: VecDot(), VecMTDot()
454: @*/
455: PetscErrorCode  VecTDot(Vec x,Vec y,PetscScalar *val)
456: {


468:   PetscLogEventBegin(VEC_TDot,x,y,0,0);
469:   (*x->ops->tdot)(x,y,val);
470:   PetscLogEventEnd(VEC_TDot,x,y,0,0);
471:   return(0);
472: }

476: /*@
477:    VecScale - Scales a vector.

479:    Not collective on Vec

481:    Input Parameters:
482: +  x - the vector
483: -  alpha - the scalar

485:    Output Parameter:
486: .  x - the scaled vector

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

492:    Level: intermediate

494:    Concepts: vector^scaling
495:    Concepts: scaling^vector

497: @*/
498: PetscErrorCode  VecScale(Vec x, PetscScalar alpha)
499: {
500:   PetscReal      norms[4] = {0.0,0.0,0.0, 0.0};
501:   PetscBool      flgs[4];
503:   PetscInt       i;

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

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

533:    Logically Collective on Vec

535:    Input Parameters:
536: +  x  - the vector
537: -  alpha - the scalar

539:    Output Parameter:
540: .  x  - the vector

542:    Note:
543:    For a vector of dimension n, VecSet() computes
544: $     x[i] = alpha, for i=1,...,n,
545:    so that all vector entries then equal the identical
546:    scalar value, alpha.  Use the more general routine
547:    VecSetValues() to set different vector entries.

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

552:    Level: beginner

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

556:    Concepts: vector^setting to constant

558: @*/
559: PetscErrorCode  VecSet(Vec x,PetscScalar alpha)
560: {
561:   PetscReal      val;

567:   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()");

570:   PetscLogEventBegin(VEC_Set,x,0,0,0);
571:   (*x->ops->set)(x,alpha);
572:   PetscLogEventEnd(VEC_Set,x,0,0,0);
573:   PetscObjectStateIncrease((PetscObject)x);

575:   /*  norms can be simply set */
576:   val  = PetscAbsScalar(alpha);
577:   PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],x->map->N * val);
578:   PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
579:   val  = PetscSqrtReal((double)x->map->N) * val;
580:   PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],val);
581:   PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],val);
582:   return(0);
583: }


588: /*@
589:    VecAXPY - Computes y = alpha x + y.

591:    Logically Collective on Vec

593:    Input Parameters:
594: +  alpha - the scalar
595: -  x, y  - the vectors

597:    Output Parameter:
598: .  y - output vector

600:    Level: intermediate

602:    Notes: x and y MUST be different vectors

604:    Concepts: vector^BLAS
605:    Concepts: BLAS

607: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY()
608: @*/
609: PetscErrorCode  VecAXPY(Vec y,PetscScalar alpha,Vec x)
610: {

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

623:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
624:   (*y->ops->axpy)(y,alpha,x);
625:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
626:   PetscObjectStateIncrease((PetscObject)y);
627:   return(0);
628: }

632: /*@
633:    VecAXPBY - Computes y = alpha x + beta y.

635:    Logically Collective on Vec

637:    Input Parameters:
638: +  alpha,beta - the scalars
639: -  x, y  - the vectors

641:    Output Parameter:
642: .  y - output vector

644:    Level: intermediate

646:    Notes: x and y MUST be different vectors

648:    Concepts: BLAS
649:    Concepts: vector^BLAS

651: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
652: @*/
653: PetscErrorCode  VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
654: {

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

668:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
669:   (*y->ops->axpby)(y,alpha,beta,x);
670:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
671:   PetscObjectStateIncrease((PetscObject)y);
672:   return(0);
673: }

677: /*@
678:    VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z

680:    Logically Collective on Vec

682:    Input Parameters:
683: +  alpha,beta, gamma - the scalars
684: -  x, y, z  - the vectors

686:    Output Parameter:
687: .  z - output vector

689:    Level: intermediate

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

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

695:    Concepts: BLAS
696:    Concepts: vector^BLAS

698: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
699: @*/
700: PetscErrorCode  VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
701: {

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

721:   PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
722:   (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
723:   PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
724:   PetscObjectStateIncrease((PetscObject)z);
725:   return(0);
726: }

730: /*@
731:    VecAYPX - Computes y = x + alpha y.

733:    Logically Collective on Vec

735:    Input Parameters:
736: +  alpha - the scalar
737: -  x, y  - the vectors

739:    Output Parameter:
740: .  y - output vector

742:    Level: intermediate

744:    Notes: x and y MUST be different vectors

746:    Concepts: vector^BLAS
747:    Concepts: BLAS

749: .seealso: VecAXPY(), VecWAXPY()
750: @*/
751: PetscErrorCode  VecAYPX(Vec y,PetscScalar alpha,Vec x)
752: {

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

763:   PetscLogEventBegin(VEC_AYPX,x,y,0,0);
764:    (*y->ops->aypx)(y,alpha,x);
765:   PetscLogEventEnd(VEC_AYPX,x,y,0,0);
766:   PetscObjectStateIncrease((PetscObject)y);
767:   return(0);
768: }


773: /*@
774:    VecWAXPY - Computes w = alpha x + y.

776:    Logically Collective on Vec

778:    Input Parameters:
779: +  alpha - the scalar
780: -  x, y  - the vectors

782:    Output Parameter:
783: .  w - the result

785:    Level: intermediate

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

789:    Concepts: vector^BLAS
790:    Concepts: BLAS

792: .seealso: VecAXPY(), VecAYPX(), VecAXPBY()
793: @*/
794: PetscErrorCode  VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
795: {

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

813:   PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
814:    (*w->ops->waxpy)(w,alpha,x,y);
815:   PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
816:   PetscObjectStateIncrease((PetscObject)w);
817:   return(0);
818: }


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

826:    Not Collective

828:    Input Parameters:
829: +  x - vector to insert in
830: .  ni - number of elements to add
831: .  ix - indices where to add
832: .  y - array of values
833: -  iora - either INSERT_VALUES or ADD_VALUES, where
834:    ADD_VALUES adds values to any existing entries, and
835:    INSERT_VALUES replaces existing entries with new values

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

840:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
841:    options cannot be mixed without intervening calls to the assembly
842:    routines.

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

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

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

855:    Level: beginner

857:    Concepts: vector^setting values

859: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
860:            VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
861: @*/
862: PetscErrorCode  VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
863: {

871:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
872:   (*x->ops->setvalues)(x,ni,ix,y,iora);
873:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
874:   PetscObjectStateIncrease((PetscObject)x);
875:   return(0);
876: }

880: /*@
881:    VecGetValues - Gets values from certain locations of a vector. Currently
882:           can only get values on the same processor

884:     Not Collective

886:    Input Parameters:
887: +  x - vector to get values from
888: .  ni - number of elements to get
889: -  ix - indices where to get them from (in global 1d numbering)

891:    Output Parameter:
892: .   y - array of values

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

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

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

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

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

907:    Level: beginner

909:    Concepts: vector^getting values

911: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecGetValuesLocal(),
912:            VecGetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecSetValues()
913: @*/
914: PetscErrorCode  VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
915: {

923:   (*x->ops->getvalues)(x,ni,ix,y);
924:   return(0);
925: }

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

932:    Not Collective

934:    Input Parameters:
935: +  x - vector to insert in
936: .  ni - number of blocks to add
937: .  ix - indices where to add in block count, rather than element count
938: .  y - array of values
939: -  iora - either INSERT_VALUES or ADD_VALUES, where
940:    ADD_VALUES adds values to any existing entries, and
941:    INSERT_VALUES replaces existing entries with new values

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

947:    Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
948:    options cannot be mixed without intervening calls to the assembly
949:    routines.

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

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

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

961:    Level: intermediate

963:    Concepts: vector^setting values blocked

965: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
966:            VecSetValues()
967: @*/
968: PetscErrorCode  VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
969: {

977:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
978:   (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
979:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
980:   PetscObjectStateIncrease((PetscObject)x);
981:   return(0);
982: }


987: /*@
988:    VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
989:    using a local ordering of the nodes.

991:    Not Collective

993:    Input Parameters:
994: +  x - vector to insert in
995: .  ni - number of elements to add
996: .  ix - indices where to add
997: .  y - array of values
998: -  iora - either INSERT_VALUES or ADD_VALUES, where
999:    ADD_VALUES adds values to any existing entries, and
1000:    INSERT_VALUES replaces existing entries with new values

1002:    Level: intermediate

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

1007:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
1008:    options cannot be mixed without intervening calls to the assembly
1009:    routines.

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

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

1016:    Concepts: vector^setting values with local numbering

1018: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetLocalToGlobalMapping(),
1019:            VecSetValuesBlockedLocal()
1020: @*/
1021: PetscErrorCode  VecSetValuesLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1022: {
1024:   PetscInt       lixp[128],*lix = lixp;


1032:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1033:   if (!x->ops->setvalueslocal) {
1034:     if (!x->map->mapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMapping()");
1035:     if (ni > 128) {
1036:       PetscMalloc(ni*sizeof(PetscInt),&lix);
1037:     }
1038:     ISLocalToGlobalMappingApply(x->map->mapping,ni,(PetscInt*)ix,lix);
1039:     (*x->ops->setvalues)(x,ni,lix,y,iora);
1040:     if (ni > 128) {
1041:       PetscFree(lix);
1042:     }
1043:   } else {
1044:     (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
1045:   }
1046:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1047:   PetscObjectStateIncrease((PetscObject)x);
1048:   return(0);
1049: }

1053: /*@
1054:    VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1055:    using a local ordering of the nodes.

1057:    Not Collective

1059:    Input Parameters:
1060: +  x - vector to insert in
1061: .  ni - number of blocks to add
1062: .  ix - indices where to add in block count, not element count
1063: .  y - array of values
1064: -  iora - either INSERT_VALUES or ADD_VALUES, where
1065:    ADD_VALUES adds values to any existing entries, and
1066:    INSERT_VALUES replaces existing entries with new values

1068:    Level: intermediate

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

1074:    Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1075:    options cannot be mixed without intervening calls to the assembly
1076:    routines.

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

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


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

1086: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
1087:            VecSetLocalToGlobalMappingBlock()
1088: @*/
1089: PetscErrorCode  VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1090: {
1092:   PetscInt       lixp[128],*lix = lixp;

1099:   if (!x->map->bmapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMappingBlock()");
1100:   if (ni > 128) {
1101:     PetscMalloc(ni*sizeof(PetscInt),&lix);
1102:   }

1104:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1105:   ISLocalToGlobalMappingApply(x->map->bmapping,ni,(PetscInt*)ix,lix);
1106:   (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1107:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1108:   if (ni > 128) {
1109:     PetscFree(lix);
1110:   }
1111:   PetscObjectStateIncrease((PetscObject)x);
1112:   return(0);
1113: }

1117: /*@
1118:    VecMTDot - Computes indefinite vector multiple dot products.
1119:    That is, it does NOT use the complex conjugate.

1121:    Collective on Vec

1123:    Input Parameters:
1124: +  x - one vector
1125: .  nv - number of vectors
1126: -  y - array of vectors.  Note that vectors are pointers

1128:    Output Parameter:
1129: .  val - array of the dot products

1131:    Notes for Users of Complex Numbers:
1132:    For complex vectors, VecMTDot() computes the indefinite form
1133: $      val = (x,y) = y^T x,
1134:    where y^T denotes the transpose of y.

1136:    Use VecMDot() for the inner product
1137: $      val = (x,y) = y^H x,
1138:    where y^H denotes the conjugate transpose of y.

1140:    Level: intermediate

1142:    Concepts: inner product^multiple
1143:    Concepts: vector^multiple inner products

1145: .seealso: VecMDot(), VecTDot()
1146: @*/
1147: PetscErrorCode  VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1148: {


1161:   PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1162:   (*x->ops->mtdot)(x,nv,y,val);
1163:   PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1164:   return(0);
1165: }

1169: /*@
1170:    VecMDot - Computes vector multiple dot products.

1172:    Collective on Vec

1174:    Input Parameters:
1175: +  x - one vector
1176: .  nv - number of vectors
1177: -  y - array of vectors.

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

1182:    Notes for Users of Complex Numbers:
1183:    For complex vectors, VecMDot() computes
1184: $     val = (x,y) = y^H x,
1185:    where y^H denotes the conjugate transpose of y.

1187:    Use VecMTDot() for the indefinite form
1188: $     val = (x,y) = y^T x,
1189:    where y^T denotes the transpose of y.

1191:    Level: intermediate

1193:    Concepts: inner product^multiple
1194:    Concepts: vector^multiple inner products

1196: .seealso: VecMTDot(), VecDot()
1197: @*/
1198: PetscErrorCode  VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1199: {

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

1214:   PetscLogEventBarrierBegin(VEC_MDotBarrier,x,*y,0,0,PetscObjectComm((PetscObject)x));
1215:   (*x->ops->mdot)(x,nv,y,val);
1216:   PetscLogEventBarrierEnd(VEC_MDotBarrier,x,*y,0,0,PetscObjectComm((PetscObject)x));
1217:   return(0);
1218: }

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

1225:    Logically Collective on Vec

1227:    Input Parameters:
1228: +  nv - number of scalars and x-vectors
1229: .  alpha - array of scalars
1230: .  y - one vector
1231: -  x - array of vectors

1233:    Level: intermediate

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

1237:    Concepts: BLAS

1239: .seealso: VecAXPY(), VecWAXPY(), VecAYPX()
1240: @*/
1241: PetscErrorCode  VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1242: {
1244:   PetscInt       i;

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

1259:   PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1260:   (*y->ops->maxpy)(y,nv,alpha,x);
1261:   PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1262:   PetscObjectStateIncrease((PetscObject)y);
1263:   return(0);
1264: }

1268: /*@
1269:    VecGetSubVector - Gets a vector representing part of another vector

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

1273:    Input Arguments:
1274: + X - vector from which to extract a subvector
1275: - is - index set representing portion of X to extract

1277:    Output Arguments:
1278: . Y - subvector corresponding to is

1280:    Level: advanced

1282:    Notes:
1283:    The subvector Y should be returned with VecRestoreSubVector().

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

1288: .seealso: MatGetSubMatrix()
1289: @*/
1290: PetscErrorCode  VecGetSubVector(Vec X,IS is,Vec *Y)
1291: {
1293:   Vec            Z;
1294:   PetscInt       state;

1300:   if (X->ops->getsubvector) {
1301:     (*X->ops->getsubvector)(X,is,&Z);
1302:   } else {                      /* Default implementation currently does no caching */
1303:     PetscInt  gstart,gend,start;
1304:     PetscBool contiguous,gcontiguous;
1305:     VecGetOwnershipRange(X,&gstart,&gend);
1306:     ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1307:     MPI_Allreduce(&contiguous,&gcontiguous,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1308:     if (gcontiguous) {          /* We can do a no-copy implementation */
1309:       PetscInt    n,N;
1310:       PetscScalar *x;
1311:       PetscMPIInt size;
1312:       ISGetLocalSize(is,&n);
1313:       VecGetArray(X,&x);
1314:       MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1315:       if (size == 1) {
1316:         VecCreateSeqWithArray(PetscObjectComm((PetscObject)X),1,n,x+start,&Z);
1317:       } else {
1318:         ISGetSize(is,&N);
1319:         VecCreateMPIWithArray(PetscObjectComm((PetscObject)X),1,n,N,x+start,&Z);
1320:       }
1321:       VecRestoreArray(X,&x);
1322:     } else {                    /* Have to create a scatter and do a copy */
1323:       VecScatter scatter;
1324:       PetscInt   n,N;
1325:       ISGetLocalSize(is,&n);
1326:       ISGetSize(is,&N);
1327:       VecCreate(PetscObjectComm((PetscObject)is),&Z);
1328:       VecSetSizes(Z,n,N);
1329:       VecSetType(Z,((PetscObject)X)->type_name);
1330:       VecScatterCreate(X,is,Z,NULL,&scatter);
1331:       VecScatterBegin(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1332:       VecScatterEnd(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1333:       VecScatterDestroy(&scatter);
1334:     }
1335:   }
1336:   /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1337:   if (VecGetSubVectorSavedStateId < 0) {PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);}
1338:   PetscObjectStateQuery((PetscObject)Z,&state);
1339:   PetscObjectComposedDataSetInt((PetscObject)Z,VecGetSubVectorSavedStateId,state);
1340:   *Y   = Z;
1341:   return(0);
1342: }

1346: /*@
1347:    VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()

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

1351:    Input Arguments:
1352: + X - vector from which subvector was obtained
1353: . is - index set representing the subset of X
1354: - Y - subvector being restored

1356:    Level: advanced

1358: .seealso: VecGetSubVector()
1359: @*/
1360: PetscErrorCode  VecRestoreSubVector(Vec X,IS is,Vec *Y)
1361: {

1369:   if (X->ops->restoresubvector) {
1370:     (*X->ops->restoresubvector)(X,is,Y);
1371:   } else {
1372:     PetscInt  savedstate=0,newstate;
1373:     PetscBool valid;
1374:     PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,savedstate,valid);
1375:     PetscObjectStateQuery((PetscObject)*Y,&newstate);
1376:     if (valid && savedstate < newstate) {
1377:       /* We might need to copy entries back, first check whether we have no-copy view */
1378:       PetscInt  gstart,gend,start;
1379:       PetscBool contiguous,gcontiguous;
1380:       VecGetOwnershipRange(X,&gstart,&gend);
1381:       ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1382:       MPI_Allreduce(&contiguous,&gcontiguous,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1383:       if (!gcontiguous) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Unhandled case, values have been changed and need to be copied back into X");
1384:     }
1385:     VecDestroy(Y);
1386:   }
1387:   return(0);
1388: }

1392: /*@C
1393:    VecGetArray - Returns a pointer to a contiguous array that contains this
1394:    processor's portion of the vector data. For the standard PETSc
1395:    vectors, VecGetArray() returns a pointer to the local data array and
1396:    does not use any copies. If the underlying vector data is not stored
1397:    in a contiquous array this routine will copy the data to a contiquous
1398:    array and return a pointer to that. You MUST call VecRestoreArray()
1399:    when you no longer need access to the array.

1401:    Logically Collective on Vec

1403:    Input Parameter:
1404: .  x - the vector

1406:    Output Parameter:
1407: .  a - location to put pointer to the array

1409:    Fortran Note:
1410:    This routine is used differently from Fortran 77
1411: $    Vec         x
1412: $    PetscScalar x_array(1)
1413: $    PetscOffset i_x
1414: $    PetscErrorCode ierr
1415: $       call VecGetArray(x,x_array,i_x,ierr)
1416: $
1417: $   Access first local entry in vector with
1418: $      value = x_array(i_x + 1)
1419: $
1420: $      ...... other code
1421: $       call VecRestoreArray(x,x_array,i_x,ierr)
1422:    For Fortran 90 see VecGetArrayF90()

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

1427:    Level: beginner

1429:    Concepts: vector^accessing local values

1431: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(), VecGetArray2d()
1432: @*/
1433: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1434: {

1439:   if (x->petscnative) {
1440: #if defined(PETSC_HAVE_CUSP)
1441:     if (x->valid_GPU_array == PETSC_CUSP_GPU) {
1442:       VecCUSPCopyFromGPU(x);
1443:     }
1444: #endif
1445:     *a = *((PetscScalar **)x->data);
1446:   } else {
1447:     (*x->ops->getarray)(x,a);
1448:   }
1449:   return(0);
1450: }

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

1457:    Not Collective

1459:    Input Parameters:
1460: .  x - the vector

1462:    Output Parameter:
1463: .  a - the array

1465:    Level: beginner

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

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

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

1476: .seealso: VecGetArray(), VecRestoreArray()
1477: @*/
1478: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1479: {

1484:   if (x->petscnative) {
1485: #if defined(PETSC_HAVE_CUSP)
1486:     if (x->valid_GPU_array == PETSC_CUSP_GPU) {
1487:       VecCUSPCopyFromGPU(x);
1488:     }
1489: #endif
1490:     *a = *((PetscScalar **)x->data);
1491:   } else if (x->ops->getarrayread) {
1492:     (*x->ops->getarrayread)(x,a);
1493:   } else {
1494:     (*x->ops->getarray)(x,(PetscScalar**)a);
1495:   }
1496:   return(0);
1497: }

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

1506:    Logically Collective on Vec

1508:    Input Parameter:
1509: +  x - the vectors
1510: -  n - the number of vectors

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

1515:    Fortran Note:
1516:    This routine is not supported in Fortran.

1518:    Level: intermediate

1520: .seealso: VecGetArray(), VecRestoreArrays()
1521: @*/
1522: PetscErrorCode  VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1523: {
1525:   PetscInt       i;
1526:   PetscScalar    **q;

1532:   if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1533:   PetscMalloc(n*sizeof(PetscScalar*),&q);
1534:   for (i=0; i<n; ++i) {
1535:     VecGetArray(x[i],&q[i]);
1536:   }
1537:   *a = q;
1538:   return(0);
1539: }

1543: /*@C
1544:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1545:    has been called.

1547:    Logically Collective on Vec

1549:    Input Parameters:
1550: +  x - the vector
1551: .  n - the number of vectors
1552: -  a - location of pointer to arrays obtained from VecGetArrays()

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

1560:    Fortran Note:
1561:    This routine is not supported in Fortran.

1563:    Level: intermediate

1565: .seealso: VecGetArrays(), VecRestoreArray()
1566: @*/
1567: PetscErrorCode  VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1568: {
1570:   PetscInt       i;
1571:   PetscScalar    **q = *a;


1578:   for (i=0; i<n; ++i) {
1579:     VecRestoreArray(x[i],&q[i]);
1580:   }
1581:   PetscFree(q);
1582:   return(0);
1583: }

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

1590:    Logically Collective on Vec

1592:    Input Parameters:
1593: +  x - the vector
1594: -  a - location of pointer to array obtained from VecGetArray()

1596:    Level: beginner

1598:    Notes:
1599:    For regular PETSc vectors this routine does not involve any copies. For
1600:    any special vectors that do not store local vector data in a contiguous
1601:    array, this routine will copy the data back into the underlying
1602:    vector data structure from the array obtained with VecGetArray().

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

1608:    Fortran Note:
1609:    This routine is used differently from Fortran 77
1610: $    Vec         x
1611: $    PetscScalar x_array(1)
1612: $    PetscOffset i_x
1613: $    PetscErrorCode ierr
1614: $       call VecGetArray(x,x_array,i_x,ierr)
1615: $
1616: $   Access first local entry in vector with
1617: $      value = x_array(i_x + 1)
1618: $
1619: $      ...... other code
1620: $       call VecRestoreArray(x,x_array,i_x,ierr)

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

1626: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(), VecRestoreArray2d()
1627: @*/
1628: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1629: {

1634:   if (x->petscnative) {
1635: #if defined(PETSC_HAVE_CUSP)
1636:     x->valid_GPU_array = PETSC_CUSP_CPU;
1637: #endif
1638:   } else {
1639:     (*x->ops->restorearray)(x,a);
1640:   }
1641:   if (a) *a = NULL;
1642:   PetscObjectStateIncrease((PetscObject)x);
1643:   return(0);
1644: }

1648: /*@C
1649:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

1651:    Not Collective

1653:    Input Parameters:
1654: +  vec - the vector
1655: -  array - the array

1657:    Level: beginner

1659: .seealso: VecGetArray(), VecRestoreArray()
1660: @*/
1661: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1662: {

1667:   if (x->petscnative) {
1668: #if defined(PETSC_HAVE_CUSP)
1669:     x->valid_GPU_array = PETSC_CUSP_CPU;
1670: #endif
1671:   } else if (x->ops->restorearrayread) {
1672:     (*x->ops->restorearrayread)(x,a);
1673:   } else {
1674:     (*x->ops->restorearray)(x,(PetscScalar**)a);
1675:   }
1676:   if (a) *a = NULL;
1677:   return(0);
1678: }

1682: /*@
1683:    VecPlaceArray - Allows one to replace the array in a vector with an
1684:    array provided by the user. This is useful to avoid copying an array
1685:    into a vector.

1687:    Not Collective

1689:    Input Parameters:
1690: +  vec - the vector
1691: -  array - the array

1693:    Notes:
1694:    You can return to the original array with a call to VecResetArray()

1696:    Level: developer

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

1700: @*/
1701: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
1702: {

1709:   if (vec->ops->placearray) {
1710:     (*vec->ops->placearray)(vec,array);
1711:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
1712:   PetscObjectStateIncrease((PetscObject)vec);
1713:   return(0);
1714: }

1718: /*@C
1719:    VecReplaceArray - Allows one to replace the array in a vector with an
1720:    array provided by the user. This is useful to avoid copying an array
1721:    into a vector.

1723:    Not Collective

1725:    Input Parameters:
1726: +  vec - the vector
1727: -  array - the array

1729:    Notes:
1730:    This permanently replaces the array and frees the memory associated
1731:    with the old array.

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

1736:    Not supported from Fortran

1738:    Level: developer

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

1742: @*/
1743: PetscErrorCode  VecReplaceArray(Vec vec,const PetscScalar array[])
1744: {

1750:   if (vec->ops->replacearray) {
1751:     (*vec->ops->replacearray)(vec,array);
1752:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
1753:   PetscObjectStateIncrease((PetscObject)vec);
1754:   return(0);
1755: }

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

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

1764:     Collective on Vec

1766:     Input Parameters:
1767: +   x - a vector to mimic
1768: -   n - the number of vectors to obtain

1770:     Output Parameters:
1771: +   y - Fortran90 pointer to the array of vectors
1772: -   ierr - error code

1774:     Example of Usage:
1775: .vb
1776:     Vec x
1777:     Vec, pointer :: y(:)
1778:     ....
1779:     call VecDuplicateVecsF90(x,2,y,ierr)
1780:     call VecSet(y(2),alpha,ierr)
1781:     call VecSet(y(2),alpha,ierr)
1782:     ....
1783:     call VecDestroyVecsF90(2,y,ierr)
1784: .ve

1786:     Notes:
1787:     Not yet supported for all F90 compilers

1789:     Use VecDestroyVecsF90() to free the space.

1791:     Level: beginner

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

1795: M*/

1797: /*MC
1798:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
1799:     VecGetArrayF90().

1801:     Synopsis:
1802:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

1804:     Logically Collective on Vec

1806:     Input Parameters:
1807: +   x - vector
1808: -   xx_v - the Fortran90 pointer to the array

1810:     Output Parameter:
1811: .   ierr - error code

1813:     Example of Usage:
1814: .vb
1815:     PetscScalar, pointer :: xx_v(:)
1816:     ....
1817:     call VecGetArrayF90(x,xx_v,ierr)
1818:     a = xx_v(3)
1819:     call VecRestoreArrayF90(x,xx_v,ierr)
1820: .ve

1822:     Level: beginner

1824: .seealso:  VecGetArrayF90(), VecGetArray(), VecRestoreArray(), UsingFortran

1826: M*/

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

1831:     Synopsis:
1832:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

1834:     Collective on Vec

1836:     Input Parameters:
1837: +   n - the number of vectors previously obtained
1838: -   x - pointer to array of vector pointers

1840:     Output Parameter:
1841: .   ierr - error code

1843:     Notes:
1844:     Not yet supported for all F90 compilers

1846:     Level: beginner

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

1850: M*/

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

1858:     Synopsis:
1859:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

1861:     Logically Collective on Vec

1863:     Input Parameter:
1864: .   x - vector

1866:     Output Parameters:
1867: +   xx_v - the Fortran90 pointer to the array
1868: -   ierr - error code

1870:     Example of Usage:
1871: .vb
1872:     PetscScalar, pointer :: xx_v(:)
1873:     ....
1874:     call VecGetArrayF90(x,xx_v,ierr)
1875:     a = xx_v(3)
1876:     call VecRestoreArrayF90(x,xx_v,ierr)
1877: .ve

1879:     Level: beginner

1881: .seealso:  VecRestoreArrayF90(), VecGetArray(), VecRestoreArray(), UsingFortran

1883: M*/


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

1893:    Logically Collective

1895:    Input Parameter:
1896: +  x - the vector
1897: .  m - first dimension of two dimensional array
1898: .  n - second dimension of two dimensional array
1899: .  mstart - first index you will use in first coordinate direction (often 0)
1900: -  nstart - first index in the second coordinate direction (often 0)

1902:    Output Parameter:
1903: .  a - location to put pointer to the array

1905:    Level: developer

1907:   Notes:
1908:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
1909:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
1910:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
1911:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

1917: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
1918:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
1919:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
1920: @*/
1921: PetscErrorCode  VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
1922: {
1924:   PetscInt       i,N;
1925:   PetscScalar    *aa;

1931:   VecGetLocalSize(x,&N);
1932:   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);
1933:   VecGetArray(x,&aa);

1935:   PetscMalloc(m*sizeof(PetscScalar*),a);
1936:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
1937:   *a -= mstart;
1938:   return(0);
1939: }

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

1946:    Logically Collective

1948:    Input Parameters:
1949: +  x - the vector
1950: .  m - first dimension of two dimensional array
1951: .  n - second dimension of the two dimensional array
1952: .  mstart - first index you will use in first coordinate direction (often 0)
1953: .  nstart - first index in the second coordinate direction (often 0)
1954: -  a - location of pointer to array obtained from VecGetArray2d()

1956:    Level: developer

1958:    Notes:
1959:    For regular PETSc vectors this routine does not involve any copies. For
1960:    any special vectors that do not store local vector data in a contiguous
1961:    array, this routine will copy the data back into the underlying
1962:    vector data structure from the array obtained with VecGetArray().

1964:    This routine actually zeros out the a pointer.

1966: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
1967:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
1968:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
1969: @*/
1970: PetscErrorCode  VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
1971: {
1973:   void           *dummy;

1979:   dummy = (void*)(*a + mstart);
1980:   PetscFree(dummy);
1981:   VecRestoreArray(x,NULL);
1982:   return(0);
1983: }

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

1992:    Logically Collective

1994:    Input Parameter:
1995: +  x - the vector
1996: .  m - first dimension of two dimensional array
1997: -  mstart - first index you will use in first coordinate direction (often 0)

1999:    Output Parameter:
2000: .  a - location to put pointer to the array

2002:    Level: developer

2004:   Notes:
2005:    For a vector obtained from DMCreateLocalVector() mstart are likely
2006:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2007:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2011: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2012:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2013:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2014: @*/
2015: PetscErrorCode  VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2016: {
2018:   PetscInt       N;

2024:   VecGetLocalSize(x,&N);
2025:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2026:   VecGetArray(x,a);
2027:   *a  -= mstart;
2028:   return(0);
2029: }

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

2036:    Logically Collective

2038:    Input Parameters:
2039: +  x - the vector
2040: .  m - first dimension of two dimensional array
2041: .  mstart - first index you will use in first coordinate direction (often 0)
2042: -  a - location of pointer to array obtained from VecGetArray21()

2044:    Level: developer

2046:    Notes:
2047:    For regular PETSc vectors this routine does not involve any copies. For
2048:    any special vectors that do not store local vector data in a contiguous
2049:    array, this routine will copy the data back into the underlying
2050:    vector data structure from the array obtained with VecGetArray1d().

2052:    This routine actually zeros out the a pointer.

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

2056: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2057:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2058:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2059: @*/
2060: PetscErrorCode  VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2061: {

2067:   VecRestoreArray(x,NULL);
2068:   return(0);
2069: }


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

2079:    Logically Collective

2081:    Input Parameter:
2082: +  x - the vector
2083: .  m - first dimension of three dimensional array
2084: .  n - second dimension of three dimensional array
2085: .  p - third dimension of three dimensional array
2086: .  mstart - first index you will use in first coordinate direction (often 0)
2087: .  nstart - first index in the second coordinate direction (often 0)
2088: -  pstart - first index in the third coordinate direction (often 0)

2090:    Output Parameter:
2091: .  a - location to put pointer to the array

2093:    Level: developer

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

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

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

2105: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2106:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2107:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2108: @*/
2109: PetscErrorCode  VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2110: {
2112:   PetscInt       i,N,j;
2113:   PetscScalar    *aa,**b;

2119:   VecGetLocalSize(x,&N);
2120:   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);
2121:   VecGetArray(x,&aa);

2123:   PetscMalloc(m*sizeof(PetscScalar**)+m*n*sizeof(PetscScalar*),a);
2124:   b    = (PetscScalar**)((*a) + m);
2125:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2126:   for (i=0; i<m; i++)
2127:     for (j=0; j<n; j++)
2128:       b[i*n+j] = aa + i*n*p + j*p - pstart;

2130:   *a -= mstart;
2131:   return(0);
2132: }

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

2139:    Logically Collective

2141:    Input Parameters:
2142: +  x - the vector
2143: .  m - first dimension of three dimensional array
2144: .  n - second dimension of the three dimensional array
2145: .  p - third dimension of the three dimensional array
2146: .  mstart - first index you will use in first coordinate direction (often 0)
2147: .  nstart - first index in the second coordinate direction (often 0)
2148: .  pstart - first index in the third coordinate direction (often 0)
2149: -  a - location of pointer to array obtained from VecGetArray3d()

2151:    Level: developer

2153:    Notes:
2154:    For regular PETSc vectors this routine does not involve any copies. For
2155:    any special vectors that do not store local vector data in a contiguous
2156:    array, this routine will copy the data back into the underlying
2157:    vector data structure from the array obtained with VecGetArray().

2159:    This routine actually zeros out the a pointer.

2161: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2162:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2163:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2164: @*/
2165: PetscErrorCode  VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2166: {
2168:   void           *dummy;

2174:   dummy = (void*)(*a + mstart);
2175:   PetscFree(dummy);
2176:   VecRestoreArray(x,NULL);
2177:   return(0);
2178: }

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

2187:    Logically Collective

2189:    Input Parameter:
2190: +  x - the vector
2191: .  m - first dimension of four dimensional array
2192: .  n - second dimension of four dimensional array
2193: .  p - third dimension of four dimensional array
2194: .  q - fourth dimension of four dimensional array
2195: .  mstart - first index you will use in first coordinate direction (often 0)
2196: .  nstart - first index in the second coordinate direction (often 0)
2197: .  pstart - first index in the third coordinate direction (often 0)
2198: -  qstart - first index in the fourth coordinate direction (often 0)

2200:    Output Parameter:
2201: .  a - location to put pointer to the array

2203:    Level: beginner

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

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

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

2215: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2216:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2217:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2218: @*/
2219: PetscErrorCode  VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2220: {
2222:   PetscInt       i,N,j,k;
2223:   PetscScalar    *aa,***b,**c;

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

2233:   PetscMalloc(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p*sizeof(PetscScalar*),a);
2234:   b    = (PetscScalar***)((*a) + m);
2235:   c    = (PetscScalar**)(b + m*n);
2236:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2237:   for (i=0; i<m; i++)
2238:     for (j=0; j<n; j++)
2239:       b[i*n+j] = c + i*n*p + j*p - pstart;
2240:   for (i=0; i<m; i++)
2241:     for (j=0; j<n; j++)
2242:       for (k=0; k<p; k++)
2243:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
2244:   *a -= mstart;
2245:   return(0);
2246: }

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

2253:    Logically Collective

2255:    Input Parameters:
2256: +  x - the vector
2257: .  m - first dimension of four dimensional array
2258: .  n - second dimension of the four dimensional array
2259: .  p - third dimension of the four dimensional array
2260: .  q - fourth dimension of the four dimensional array
2261: .  mstart - first index you will use in first coordinate direction (often 0)
2262: .  nstart - first index in the second coordinate direction (often 0)
2263: .  pstart - first index in the third coordinate direction (often 0)
2264: .  qstart - first index in the fourth coordinate direction (often 0)
2265: -  a - location of pointer to array obtained from VecGetArray4d()

2267:    Level: beginner

2269:    Notes:
2270:    For regular PETSc vectors this routine does not involve any copies. For
2271:    any special vectors that do not store local vector data in a contiguous
2272:    array, this routine will copy the data back into the underlying
2273:    vector data structure from the array obtained with VecGetArray().

2275:    This routine actually zeros out the a pointer.

2277: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2278:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2279:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2280: @*/
2281: PetscErrorCode  VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2282: {
2284:   void           *dummy;

2290:   dummy = (void*)(*a + mstart);
2291:   PetscFree(dummy);
2292:   VecRestoreArray(x,NULL);
2293:   return(0);
2294: }