Actual source code: rvector.c

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

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

 50:    Logically Collective on Vec

 52:    Input Parameters:
 53: .  x, y  - the vectors

 55:    Output Parameter:
 56: .  max - the result

 58:    Level: advanced

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

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


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

 84: /*@
 85:    VecDot - Computes the vector dot product.

 87:    Collective on Vec

 89:    Input Parameters:
 90: .  x, y - the vectors

 92:    Output Parameter:
 93: .  val - the dot product

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

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

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

111:    Level: intermediate

113:    Concepts: inner product
114:    Concepts: vector^inner product

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


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

139: /*@
140:    VecDotRealPart - Computes the real part of the vector dot product.

142:    Collective on Vec

144:    Input Parameters:
145: .  x, y - the vectors

147:    Output Parameter:
148: .  val - the real part of the dot product;

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

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

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

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

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

165:    Level: intermediate

167:    Concepts: inner product
168:    Concepts: vector^inner product

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

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

185: /*@
186:    VecNorm  - Computes the vector norm.

188:    Collective on Vec

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

196:    Output Parameter:
197: .  val - the norm

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

204:    Level: intermediate

206:    Performance Issues:
207: $    per-processor memory bandwidth
208: $    interprocessor latency
209: $    work load inbalance that causes certain processes to arrive much earlier than others

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

216:    Concepts: norm
217:    Concepts: vector^norm

219: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNormAvailable(),
220:           VecNormBegin(), VecNormEnd()

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

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

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

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

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

256:    Not Collective

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

264:    Output Parameter:
265: +  available - PETSC_TRUE if the val returned is valid
266: -  val - the norm

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

273:    Level: intermediate

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

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

285:    Concepts: norm
286:    Concepts: vector^norm

288: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNorm()
289:           VecNormBegin(), VecNormEnd()

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


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

310: /*@
311:    VecNormalize - Normalizes a vector by 2-norm.

313:    Collective on Vec

315:    Input Parameters:
316: +  x - the vector

318:    Output Parameter:
319: .  x - the normalized vector
320: -  val - the vector norm before normalization

322:    Level: intermediate

324:    Concepts: vector^normalizing
325:    Concepts: normalizing^vector

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

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

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

354:    Collective on Vec

356:    Input Parameter:
357: .  x - the vector

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

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

366:    Returns the smallest index with the maximum value
367:    Level: intermediate

369:    Concepts: maximum^of vector
370:    Concepts: vector^maximum value

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

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

390: /*@
391:    VecMin - Determines the minimum vector component and its location.

393:    Collective on Vec

395:    Input Parameters:
396: .  x - the vector

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

402:    Level: intermediate

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

407:    This returns the smallest index with the minumum value

409:    Concepts: minimum^of vector
410:    Concepts: vector^minimum entry

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

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

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

434:    Collective on Vec

436:    Input Parameters:
437: .  x, y - the vectors

439:    Output Parameter:
440: .  val - the dot product

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

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

451:    Level: intermediate

453:    Concepts: inner product^non-Hermitian
454:    Concepts: vector^inner product
455:    Concepts: non-Hermitian inner product

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


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

480: /*@
481:    VecScale - Scales a vector.

483:    Not collective on Vec

485:    Input Parameters:
486: +  x - the vector
487: -  alpha - the scalar

489:    Output Parameter:
490: .  x - the scaled vector

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

496:    Level: intermediate

498:    Concepts: vector^scaling
499:    Concepts: scaling^vector

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

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

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

537:    Logically Collective on Vec

539:    Input Parameters:
540: +  x  - the vector
541: -  alpha - the scalar

543:    Output Parameter:
544: .  x  - the vector

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

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

556:    Level: beginner

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

560:    Concepts: vector^setting to constant

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

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

574:   PetscLogEventBegin(VEC_Set,x,0,0,0);
575:   (*x->ops->set)(x,alpha);
576:   PetscLogEventEnd(VEC_Set,x,0,0,0);
577:   PetscObjectStateIncrease((PetscObject)x);

579:   /*  norms can be simply set */
580:   val  = PetscAbsScalar(alpha);
581:   PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],x->map->N * val);
582:   PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
583:   val  = PetscSqrtReal((PetscReal)x->map->N) * val;
584:   PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],val);
585:   PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],val);
586:   return(0);
587: }


592: /*@
593:    VecAXPY - Computes y = alpha x + y.

595:    Logically Collective on Vec

597:    Input Parameters:
598: +  alpha - the scalar
599: -  x, y  - the vectors

601:    Output Parameter:
602: .  y - output vector

604:    Level: intermediate

606:    Notes: x and y MUST be different vectors

608:    Concepts: vector^BLAS
609:    Concepts: BLAS

611: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY()
612: @*/
613: PetscErrorCode  VecAXPY(Vec y,PetscScalar alpha,Vec x)
614: {

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

627:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
628:   (*y->ops->axpy)(y,alpha,x);
629:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
630:   PetscObjectStateIncrease((PetscObject)y);
631:   return(0);
632: }

636: /*@
637:    VecAXPBY - Computes y = alpha x + beta y.

639:    Logically Collective on Vec

641:    Input Parameters:
642: +  alpha,beta - the scalars
643: -  x, y  - the vectors

645:    Output Parameter:
646: .  y - output vector

648:    Level: intermediate

650:    Notes: x and y MUST be different vectors

652:    Concepts: BLAS
653:    Concepts: vector^BLAS

655: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
656: @*/
657: PetscErrorCode  VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
658: {

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

672:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
673:   (*y->ops->axpby)(y,alpha,beta,x);
674:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
675:   PetscObjectStateIncrease((PetscObject)y);
676:   return(0);
677: }

681: /*@
682:    VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z

684:    Logically Collective on Vec

686:    Input Parameters:
687: +  alpha,beta, gamma - the scalars
688: -  x, y, z  - the vectors

690:    Output Parameter:
691: .  z - output vector

693:    Level: intermediate

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

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

699:    Concepts: BLAS
700:    Concepts: vector^BLAS

702: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
703: @*/
704: PetscErrorCode  VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
705: {

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

725:   PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
726:   (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
727:   PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
728:   PetscObjectStateIncrease((PetscObject)z);
729:   return(0);
730: }

734: /*@
735:    VecAYPX - Computes y = x + alpha y.

737:    Logically Collective on Vec

739:    Input Parameters:
740: +  alpha - the scalar
741: -  x, y  - the vectors

743:    Output Parameter:
744: .  y - output vector

746:    Level: intermediate

748:    Notes: x and y MUST be different vectors

750:    Concepts: vector^BLAS
751:    Concepts: BLAS

753: .seealso: VecAXPY(), VecWAXPY()
754: @*/
755: PetscErrorCode  VecAYPX(Vec y,PetscScalar alpha,Vec x)
756: {

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

767:   PetscLogEventBegin(VEC_AYPX,x,y,0,0);
768:    (*y->ops->aypx)(y,alpha,x);
769:   PetscLogEventEnd(VEC_AYPX,x,y,0,0);
770:   PetscObjectStateIncrease((PetscObject)y);
771:   return(0);
772: }


777: /*@
778:    VecWAXPY - Computes w = alpha x + y.

780:    Logically Collective on Vec

782:    Input Parameters:
783: +  alpha - the scalar
784: -  x, y  - the vectors

786:    Output Parameter:
787: .  w - the result

789:    Level: intermediate

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

793:    Concepts: vector^BLAS
794:    Concepts: BLAS

796: .seealso: VecAXPY(), VecAYPX(), VecAXPBY()
797: @*/
798: PetscErrorCode  VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
799: {

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

817:   PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
818:    (*w->ops->waxpy)(w,alpha,x,y);
819:   PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
820:   PetscObjectStateIncrease((PetscObject)w);
821:   return(0);
822: }


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

830:    Not Collective

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

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

844:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
845:    options cannot be mixed without intervening calls to the assembly
846:    routines.

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

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

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

859:    Level: beginner

861:    Concepts: vector^setting values

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

875:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
876:   (*x->ops->setvalues)(x,ni,ix,y,iora);
877:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
878:   PetscObjectStateIncrease((PetscObject)x);
879:   return(0);
880: }

884: /*@
885:    VecGetValues - Gets values from certain locations of a vector. Currently
886:           can only get values on the same processor

888:     Not Collective

890:    Input Parameters:
891: +  x - vector to get values from
892: .  ni - number of elements to get
893: -  ix - indices where to get them from (in global 1d numbering)

895:    Output Parameter:
896: .   y - array of values

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

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

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

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

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

911:    Level: beginner

913:    Concepts: vector^getting values

915: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecGetValuesLocal(),
916:            VecGetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecSetValues()
917: @*/
918: PetscErrorCode  VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
919: {

927:   (*x->ops->getvalues)(x,ni,ix,y);
928:   return(0);
929: }

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

936:    Not Collective

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

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

951:    Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
952:    options cannot be mixed without intervening calls to the assembly
953:    routines.

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

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

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

965:    Level: intermediate

967:    Concepts: vector^setting values blocked

969: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
970:            VecSetValues()
971: @*/
972: PetscErrorCode  VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
973: {

981:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
982:   (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
983:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
984:   PetscObjectStateIncrease((PetscObject)x);
985:   return(0);
986: }


991: /*@
992:    VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
993:    using a local ordering of the nodes.

995:    Not Collective

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

1006:    Level: intermediate

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

1011:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
1012:    options cannot be mixed without intervening calls to the assembly
1013:    routines.

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

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

1020:    Concepts: vector^setting values with local numbering

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


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

1057: /*@
1058:    VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1059:    using a local ordering of the nodes.

1061:    Not Collective

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

1072:    Level: intermediate

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

1078:    Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1079:    options cannot be mixed without intervening calls to the assembly
1080:    routines.

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

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


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

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

1103:   if (!x->map->bmapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMappingBlock()");
1104:   if (ni > 128) {
1105:     PetscMalloc1(ni,&lix);
1106:   }

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

1121: /*@
1122:    VecMTDot - Computes indefinite vector multiple dot products.
1123:    That is, it does NOT use the complex conjugate.

1125:    Collective on Vec

1127:    Input Parameters:
1128: +  x - one vector
1129: .  nv - number of vectors
1130: -  y - array of vectors.  Note that vectors are pointers

1132:    Output Parameter:
1133: .  val - array of the dot products

1135:    Notes for Users of Complex Numbers:
1136:    For complex vectors, VecMTDot() computes the indefinite form
1137: $      val = (x,y) = y^T x,
1138:    where y^T denotes the transpose of y.

1140:    Use VecMDot() for the inner product
1141: $      val = (x,y) = y^H x,
1142:    where y^H denotes the conjugate transpose of y.

1144:    Level: intermediate

1146:    Concepts: inner product^multiple
1147:    Concepts: vector^multiple inner products

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


1165:   PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1166:   (*x->ops->mtdot)(x,nv,y,val);
1167:   PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1168:   return(0);
1169: }

1173: /*@
1174:    VecMDot - Computes vector multiple dot products.

1176:    Collective on Vec

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

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

1186:    Notes for Users of Complex Numbers:
1187:    For complex vectors, VecMDot() computes
1188: $     val = (x,y) = y^H x,
1189:    where y^H denotes the conjugate transpose of y.

1191:    Use VecMTDot() for the indefinite form
1192: $     val = (x,y) = y^T x,
1193:    where y^T denotes the transpose of y.

1195:    Level: intermediate

1197:    Concepts: inner product^multiple
1198:    Concepts: vector^multiple inner products

1200: .seealso: VecMTDot(), VecDot()
1201: @*/
1202: PetscErrorCode  VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1203: {

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

1218:   PetscLogEventBarrierBegin(VEC_MDotBarrier,x,*y,0,0,PetscObjectComm((PetscObject)x));
1219:   (*x->ops->mdot)(x,nv,y,val);
1220:   PetscLogEventBarrierEnd(VEC_MDotBarrier,x,*y,0,0,PetscObjectComm((PetscObject)x));
1221:   return(0);
1222: }

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

1229:    Logically Collective on Vec

1231:    Input Parameters:
1232: +  nv - number of scalars and x-vectors
1233: .  alpha - array of scalars
1234: .  y - one vector
1235: -  x - array of vectors

1237:    Level: intermediate

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

1241:    Concepts: BLAS

1243: .seealso: VecAXPY(), VecWAXPY(), VecAYPX()
1244: @*/
1245: PetscErrorCode  VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1246: {
1248:   PetscInt       i;

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

1263:   PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1264:   (*y->ops->maxpy)(y,nv,alpha,x);
1265:   PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1266:   PetscObjectStateIncrease((PetscObject)y);
1267:   return(0);
1268: }

1272: /*@
1273:    VecGetSubVector - Gets a vector representing part of another vector

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

1277:    Input Arguments:
1278: + X - vector from which to extract a subvector
1279: - is - index set representing portion of X to extract

1281:    Output Arguments:
1282: . Y - subvector corresponding to is

1284:    Level: advanced

1286:    Notes:
1287:    The subvector Y should be returned with VecRestoreSubVector().

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

1292: .seealso: MatGetSubMatrix()
1293: @*/
1294: PetscErrorCode  VecGetSubVector(Vec X,IS is,Vec *Y)
1295: {
1296:   PetscErrorCode   ierr;
1297:   Vec              Z;
1298:   PetscObjectState state;

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

1352: /*@
1353:    VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()

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

1357:    Input Arguments:
1358: + X - vector from which subvector was obtained
1359: . is - index set representing the subset of X
1360: - Y - subvector being restored

1362:    Level: advanced

1364: .seealso: VecGetSubVector()
1365: @*/
1366: PetscErrorCode  VecRestoreSubVector(Vec X,IS is,Vec *Y)
1367: {

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

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

1407:    Logically Collective on Vec

1409:    Input Parameter:
1410: .  x - the vector

1412:    Output Parameter:
1413: .  a - location to put pointer to the array

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

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

1433:    Level: beginner

1435:    Concepts: vector^accessing local values

1437: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(), VecGetArray2d()
1438: @*/
1439: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1440: {

1445:   if (x->petscnative) {
1446: #if defined(PETSC_HAVE_CUSP)
1447:     if (x->valid_GPU_array == PETSC_CUSP_GPU) {
1448:       VecCUSPCopyFromGPU(x);
1449:     }
1450: #endif
1451: #if defined(PETSC_HAVE_VIENNACL)
1452:     if (x->valid_GPU_array == PETSC_VIENNACL_GPU) {
1453:       VecViennaCLCopyFromGPU(x);
1454:     }
1455: #endif
1456:     *a = *((PetscScalar **)x->data);
1457:   } else {
1458:     (*x->ops->getarray)(x,a);
1459:   }
1460:   return(0);
1461: }

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

1468:    Not Collective

1470:    Input Parameters:
1471: .  x - the vector

1473:    Output Parameter:
1474: .  a - the array

1476:    Level: beginner

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

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

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

1487: .seealso: VecGetArray(), VecRestoreArray()
1488: @*/
1489: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1490: {

1495:   if (x->petscnative) {
1496: #if defined(PETSC_HAVE_CUSP)
1497:     if (x->valid_GPU_array == PETSC_CUSP_GPU) {
1498:       VecCUSPCopyFromGPU(x);
1499:     }
1500: #endif
1501: #if defined(PETSC_HAVE_VIENNACL)
1502:     if (x->valid_GPU_array == PETSC_VIENNACL_GPU) {
1503:       VecViennaCLCopyFromGPU(x);
1504:     }
1505: #endif
1506:     *a = *((PetscScalar **)x->data);
1507:   } else if (x->ops->getarrayread) {
1508:     (*x->ops->getarrayread)(x,a);
1509:   } else {
1510:     (*x->ops->getarray)(x,(PetscScalar**)a);
1511:   }
1512:   return(0);
1513: }

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

1522:    Logically Collective on Vec

1524:    Input Parameter:
1525: +  x - the vectors
1526: -  n - the number of vectors

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

1531:    Fortran Note:
1532:    This routine is not supported in Fortran.

1534:    Level: intermediate

1536: .seealso: VecGetArray(), VecRestoreArrays()
1537: @*/
1538: PetscErrorCode  VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1539: {
1541:   PetscInt       i;
1542:   PetscScalar    **q;

1548:   if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1549:   PetscMalloc1(n,&q);
1550:   for (i=0; i<n; ++i) {
1551:     VecGetArray(x[i],&q[i]);
1552:   }
1553:   *a = q;
1554:   return(0);
1555: }

1559: /*@C
1560:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1561:    has been called.

1563:    Logically Collective on Vec

1565:    Input Parameters:
1566: +  x - the vector
1567: .  n - the number of vectors
1568: -  a - location of pointer to arrays obtained from VecGetArrays()

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

1576:    Fortran Note:
1577:    This routine is not supported in Fortran.

1579:    Level: intermediate

1581: .seealso: VecGetArrays(), VecRestoreArray()
1582: @*/
1583: PetscErrorCode  VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1584: {
1586:   PetscInt       i;
1587:   PetscScalar    **q = *a;


1594:   for (i=0; i<n; ++i) {
1595:     VecRestoreArray(x[i],&q[i]);
1596:   }
1597:   PetscFree(q);
1598:   return(0);
1599: }

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

1606:    Logically Collective on Vec

1608:    Input Parameters:
1609: +  x - the vector
1610: -  a - location of pointer to array obtained from VecGetArray()

1612:    Level: beginner

1614:    Notes:
1615:    For regular PETSc vectors this routine does not involve any copies. For
1616:    any special vectors that do not store local vector data in a contiguous
1617:    array, this routine will copy the data back into the underlying
1618:    vector data structure from the array obtained with VecGetArray().

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

1624:    Fortran Note:
1625:    This routine is used differently from Fortran 77
1626: $    Vec         x
1627: $    PetscScalar x_array(1)
1628: $    PetscOffset i_x
1629: $    PetscErrorCode ierr
1630: $       call VecGetArray(x,x_array,i_x,ierr)
1631: $
1632: $   Access first local entry in vector with
1633: $      value = x_array(i_x + 1)
1634: $
1635: $      ...... other code
1636: $       call VecRestoreArray(x,x_array,i_x,ierr)

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

1642: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(), VecRestoreArray2d()
1643: @*/
1644: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1645: {

1650:   if (x->petscnative) {
1651: #if defined(PETSC_HAVE_CUSP)
1652:     x->valid_GPU_array = PETSC_CUSP_CPU;
1653: #endif
1654: #if defined(PETSC_HAVE_VIENNACL)
1655:     x->valid_GPU_array = PETSC_VIENNACL_CPU;
1656: #endif
1657:   } else {
1658:     (*x->ops->restorearray)(x,a);
1659:   }
1660:   if (a) *a = NULL;
1661:   PetscObjectStateIncrease((PetscObject)x);
1662:   return(0);
1663: }

1667: /*@C
1668:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

1670:    Not Collective

1672:    Input Parameters:
1673: +  vec - the vector
1674: -  array - the array

1676:    Level: beginner

1678: .seealso: VecGetArray(), VecRestoreArray()
1679: @*/
1680: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1681: {

1686:   if (x->petscnative) {
1687: #if defined(PETSC_HAVE_VIENNACL)
1688:     x->valid_GPU_array = PETSC_VIENNACL_CPU;
1689: #endif
1690:   } else if (x->ops->restorearrayread) {
1691:     (*x->ops->restorearrayread)(x,a);
1692:   } else {
1693:     (*x->ops->restorearray)(x,(PetscScalar**)a);
1694:   }
1695:   if (a) *a = NULL;
1696:   return(0);
1697: }

1701: /*@
1702:    VecPlaceArray - Allows one to replace the array in a vector with an
1703:    array provided by the user. This is useful to avoid copying an array
1704:    into a vector.

1706:    Not Collective

1708:    Input Parameters:
1709: +  vec - the vector
1710: -  array - the array

1712:    Notes:
1713:    You can return to the original array with a call to VecResetArray()

1715:    Level: developer

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

1719: @*/
1720: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
1721: {

1728:   if (vec->ops->placearray) {
1729:     (*vec->ops->placearray)(vec,array);
1730:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
1731:   PetscObjectStateIncrease((PetscObject)vec);
1732:   return(0);
1733: }

1737: /*@C
1738:    VecReplaceArray - Allows one to replace the array in a vector with an
1739:    array provided by the user. This is useful to avoid copying an array
1740:    into a vector.

1742:    Not Collective

1744:    Input Parameters:
1745: +  vec - the vector
1746: -  array - the array

1748:    Notes:
1749:    This permanently replaces the array and frees the memory associated
1750:    with the old array.

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

1755:    Not supported from Fortran

1757:    Level: developer

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

1761: @*/
1762: PetscErrorCode  VecReplaceArray(Vec vec,const PetscScalar array[])
1763: {

1769:   if (vec->ops->replacearray) {
1770:     (*vec->ops->replacearray)(vec,array);
1771:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
1772:   PetscObjectStateIncrease((PetscObject)vec);
1773:   return(0);
1774: }

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

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

1783:     Collective on Vec

1785:     Input Parameters:
1786: +   x - a vector to mimic
1787: -   n - the number of vectors to obtain

1789:     Output Parameters:
1790: +   y - Fortran90 pointer to the array of vectors
1791: -   ierr - error code

1793:     Example of Usage:
1794: .vb
1795:     Vec x
1796:     Vec, pointer :: y(:)
1797:     ....
1798:     call VecDuplicateVecsF90(x,2,y,ierr)
1799:     call VecSet(y(2),alpha,ierr)
1800:     call VecSet(y(2),alpha,ierr)
1801:     ....
1802:     call VecDestroyVecsF90(2,y,ierr)
1803: .ve

1805:     Notes:
1806:     Not yet supported for all F90 compilers

1808:     Use VecDestroyVecsF90() to free the space.

1810:     Level: beginner

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

1814: M*/

1816: /*MC
1817:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
1818:     VecGetArrayF90().

1820:     Synopsis:
1821:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

1823:     Logically Collective on Vec

1825:     Input Parameters:
1826: +   x - vector
1827: -   xx_v - the Fortran90 pointer to the array

1829:     Output Parameter:
1830: .   ierr - error code

1832:     Example of Usage:
1833: .vb
1834:     PetscScalar, pointer :: xx_v(:)
1835:     ....
1836:     call VecGetArrayF90(x,xx_v,ierr)
1837:     a = xx_v(3)
1838:     call VecRestoreArrayF90(x,xx_v,ierr)
1839: .ve

1841:     Level: beginner

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

1845: M*/

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

1850:     Synopsis:
1851:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

1853:     Collective on Vec

1855:     Input Parameters:
1856: +   n - the number of vectors previously obtained
1857: -   x - pointer to array of vector pointers

1859:     Output Parameter:
1860: .   ierr - error code

1862:     Notes:
1863:     Not yet supported for all F90 compilers

1865:     Level: beginner

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

1869: M*/

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

1877:     Synopsis:
1878:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

1880:     Logically Collective on Vec

1882:     Input Parameter:
1883: .   x - vector

1885:     Output Parameters:
1886: +   xx_v - the Fortran90 pointer to the array
1887: -   ierr - error code

1889:     Example of Usage:
1890: .vb
1891:     PetscScalar, pointer :: xx_v(:)
1892:     ....
1893:     call VecGetArrayF90(x,xx_v,ierr)
1894:     a = xx_v(3)
1895:     call VecRestoreArrayF90(x,xx_v,ierr)
1896: .ve

1898:     Level: beginner

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

1902: M*/


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

1912:    Logically Collective

1914:    Input Parameter:
1915: +  x - the vector
1916: .  m - first dimension of two dimensional array
1917: .  n - second dimension of two dimensional array
1918: .  mstart - first index you will use in first coordinate direction (often 0)
1919: -  nstart - first index in the second coordinate direction (often 0)

1921:    Output Parameter:
1922: .  a - location to put pointer to the array

1924:    Level: developer

1926:   Notes:
1927:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
1928:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
1929:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
1930:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

1936: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
1937:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
1938:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
1939: @*/
1940: PetscErrorCode  VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
1941: {
1943:   PetscInt       i,N;
1944:   PetscScalar    *aa;

1950:   VecGetLocalSize(x,&N);
1951:   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);
1952:   VecGetArray(x,&aa);

1954:   PetscMalloc1(m,a);
1955:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
1956:   *a -= mstart;
1957:   return(0);
1958: }

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

1965:    Logically Collective

1967:    Input Parameters:
1968: +  x - the vector
1969: .  m - first dimension of two dimensional array
1970: .  n - second dimension of the two dimensional array
1971: .  mstart - first index you will use in first coordinate direction (often 0)
1972: .  nstart - first index in the second coordinate direction (often 0)
1973: -  a - location of pointer to array obtained from VecGetArray2d()

1975:    Level: developer

1977:    Notes:
1978:    For regular PETSc vectors this routine does not involve any copies. For
1979:    any special vectors that do not store local vector data in a contiguous
1980:    array, this routine will copy the data back into the underlying
1981:    vector data structure from the array obtained with VecGetArray().

1983:    This routine actually zeros out the a pointer.

1985: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
1986:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
1987:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
1988: @*/
1989: PetscErrorCode  VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
1990: {
1992:   void           *dummy;

1998:   dummy = (void*)(*a + mstart);
1999:   PetscFree(dummy);
2000:   VecRestoreArray(x,NULL);
2001:   return(0);
2002: }

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

2011:    Logically Collective

2013:    Input Parameter:
2014: +  x - the vector
2015: .  m - first dimension of two dimensional array
2016: -  mstart - first index you will use in first coordinate direction (often 0)

2018:    Output Parameter:
2019: .  a - location to put pointer to the array

2021:    Level: developer

2023:   Notes:
2024:    For a vector obtained from DMCreateLocalVector() mstart are likely
2025:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2026:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2030: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2031:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2032:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2033: @*/
2034: PetscErrorCode  VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2035: {
2037:   PetscInt       N;

2043:   VecGetLocalSize(x,&N);
2044:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2045:   VecGetArray(x,a);
2046:   *a  -= mstart;
2047:   return(0);
2048: }

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

2055:    Logically Collective

2057:    Input Parameters:
2058: +  x - the vector
2059: .  m - first dimension of two dimensional array
2060: .  mstart - first index you will use in first coordinate direction (often 0)
2061: -  a - location of pointer to array obtained from VecGetArray21()

2063:    Level: developer

2065:    Notes:
2066:    For regular PETSc vectors this routine does not involve any copies. For
2067:    any special vectors that do not store local vector data in a contiguous
2068:    array, this routine will copy the data back into the underlying
2069:    vector data structure from the array obtained with VecGetArray1d().

2071:    This routine actually zeros out the a pointer.

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

2075: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2076:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2077:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2078: @*/
2079: PetscErrorCode  VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2080: {

2086:   VecRestoreArray(x,NULL);
2087:   return(0);
2088: }


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

2098:    Logically Collective

2100:    Input Parameter:
2101: +  x - the vector
2102: .  m - first dimension of three dimensional array
2103: .  n - second dimension of three dimensional array
2104: .  p - third dimension of three dimensional array
2105: .  mstart - first index you will use in first coordinate direction (often 0)
2106: .  nstart - first index in the second coordinate direction (often 0)
2107: -  pstart - first index in the third coordinate direction (often 0)

2109:    Output Parameter:
2110: .  a - location to put pointer to the array

2112:    Level: developer

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

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

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

2124: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2125:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2126:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2127: @*/
2128: PetscErrorCode  VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2129: {
2131:   PetscInt       i,N,j;
2132:   PetscScalar    *aa,**b;

2138:   VecGetLocalSize(x,&N);
2139:   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);
2140:   VecGetArray(x,&aa);

2142:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
2143:   b    = (PetscScalar**)((*a) + m);
2144:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2145:   for (i=0; i<m; i++)
2146:     for (j=0; j<n; j++)
2147:       b[i*n+j] = aa + i*n*p + j*p - pstart;

2149:   *a -= mstart;
2150:   return(0);
2151: }

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

2158:    Logically Collective

2160:    Input Parameters:
2161: +  x - the vector
2162: .  m - first dimension of three dimensional array
2163: .  n - second dimension of the three dimensional array
2164: .  p - third dimension of the three dimensional array
2165: .  mstart - first index you will use in first coordinate direction (often 0)
2166: .  nstart - first index in the second coordinate direction (often 0)
2167: .  pstart - first index in the third coordinate direction (often 0)
2168: -  a - location of pointer to array obtained from VecGetArray3d()

2170:    Level: developer

2172:    Notes:
2173:    For regular PETSc vectors this routine does not involve any copies. For
2174:    any special vectors that do not store local vector data in a contiguous
2175:    array, this routine will copy the data back into the underlying
2176:    vector data structure from the array obtained with VecGetArray().

2178:    This routine actually zeros out the a pointer.

2180: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2181:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2182:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2183: @*/
2184: PetscErrorCode  VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2185: {
2187:   void           *dummy;

2193:   dummy = (void*)(*a + mstart);
2194:   PetscFree(dummy);
2195:   VecRestoreArray(x,NULL);
2196:   return(0);
2197: }

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

2206:    Logically Collective

2208:    Input Parameter:
2209: +  x - the vector
2210: .  m - first dimension of four dimensional array
2211: .  n - second dimension of four dimensional array
2212: .  p - third dimension of four dimensional array
2213: .  q - fourth dimension of four dimensional array
2214: .  mstart - first index you will use in first coordinate direction (often 0)
2215: .  nstart - first index in the second coordinate direction (often 0)
2216: .  pstart - first index in the third coordinate direction (often 0)
2217: -  qstart - first index in the fourth coordinate direction (often 0)

2219:    Output Parameter:
2220: .  a - location to put pointer to the array

2222:    Level: beginner

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

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

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

2234: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2235:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2236:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2237: @*/
2238: PetscErrorCode  VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2239: {
2241:   PetscInt       i,N,j,k;
2242:   PetscScalar    *aa,***b,**c;

2248:   VecGetLocalSize(x,&N);
2249:   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);
2250:   VecGetArray(x,&aa);

2252:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
2253:   b    = (PetscScalar***)((*a) + m);
2254:   c    = (PetscScalar**)(b + m*n);
2255:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2256:   for (i=0; i<m; i++)
2257:     for (j=0; j<n; j++)
2258:       b[i*n+j] = c + i*n*p + j*p - pstart;
2259:   for (i=0; i<m; i++)
2260:     for (j=0; j<n; j++)
2261:       for (k=0; k<p; k++)
2262:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
2263:   *a -= mstart;
2264:   return(0);
2265: }

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

2272:    Logically Collective

2274:    Input Parameters:
2275: +  x - the vector
2276: .  m - first dimension of four dimensional array
2277: .  n - second dimension of the four dimensional array
2278: .  p - third dimension of the four dimensional array
2279: .  q - fourth dimension of the four dimensional array
2280: .  mstart - first index you will use in first coordinate direction (often 0)
2281: .  nstart - first index in the second coordinate direction (often 0)
2282: .  pstart - first index in the third coordinate direction (often 0)
2283: .  qstart - first index in the fourth coordinate direction (often 0)
2284: -  a - location of pointer to array obtained from VecGetArray4d()

2286:    Level: beginner

2288:    Notes:
2289:    For regular PETSc vectors this routine does not involve any copies. For
2290:    any special vectors that do not store local vector data in a contiguous
2291:    array, this routine will copy the data back into the underlying
2292:    vector data structure from the array obtained with VecGetArray().

2294:    This routine actually zeros out the a pointer.

2296: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2297:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2298:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2299: @*/
2300: PetscErrorCode  VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2301: {
2303:   void           *dummy;

2309:   dummy = (void*)(*a + mstart);
2310:   PetscFree(dummy);
2311:   VecRestoreArray(x,NULL);
2312:   return(0);
2313: }