Actual source code: rvector.c

petsc-master 2014-12-20
Report Typos and Errors
  2: /*
  3:      Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
  4:    These are the vector functions the user calls.
  5: */
  6: #include <petsc-private/vecimpl.h>       /*I  "petscvec.h"   I*/
  7: static PetscInt VecGetSubVectorSavedStateId = -1;

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

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

 23: #if defined(PETSC_HAVE_CUSP)
 24:   if ((vec->petscnative || vec->ops->getarray) && (vec->valid_GPU_array == PETSC_CUSP_CPU || vec->valid_GPU_array == PETSC_CUSP_BOTH)) {
 25: #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 (if |alpha|*N not too large) */
580:   val  = PetscAbsScalar(alpha);
581:   if (x->map->N == 0) {
582:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],0.0l);
583:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],0.0);
584:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],0.0);
585:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],0.0);
586:   } else if (val > PETSC_MAX_REAL/x->map->N) {
587:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
588:   } else {
589:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],x->map->N * val);
590:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
591:     val  = PetscSqrtReal((PetscReal)x->map->N) * val;
592:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],val);
593:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],val);
594:   }
595:   return(0);
596: }


601: /*@
602:    VecAXPY - Computes y = alpha x + y.

604:    Logically Collective on Vec

606:    Input Parameters:
607: +  alpha - the scalar
608: -  x, y  - the vectors

610:    Output Parameter:
611: .  y - output vector

613:    Level: intermediate

615:    Notes: x and y MUST be different vectors

617:    Concepts: vector^BLAS
618:    Concepts: BLAS

620: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY()
621: @*/
622: PetscErrorCode  VecAXPY(Vec y,PetscScalar alpha,Vec x)
623: {

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

636:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
637:   (*y->ops->axpy)(y,alpha,x);
638:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
639:   PetscObjectStateIncrease((PetscObject)y);
640:   return(0);
641: }

645: /*@
646:    VecAXPBY - Computes y = alpha x + beta y.

648:    Logically Collective on Vec

650:    Input Parameters:
651: +  alpha,beta - the scalars
652: -  x, y  - the vectors

654:    Output Parameter:
655: .  y - output vector

657:    Level: intermediate

659:    Notes: x and y MUST be different vectors

661:    Concepts: BLAS
662:    Concepts: vector^BLAS

664: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
665: @*/
666: PetscErrorCode  VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
667: {

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

681:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
682:   (*y->ops->axpby)(y,alpha,beta,x);
683:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
684:   PetscObjectStateIncrease((PetscObject)y);
685:   return(0);
686: }

690: /*@
691:    VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z

693:    Logically Collective on Vec

695:    Input Parameters:
696: +  alpha,beta, gamma - the scalars
697: -  x, y, z  - the vectors

699:    Output Parameter:
700: .  z - output vector

702:    Level: intermediate

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

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

708:    Concepts: BLAS
709:    Concepts: vector^BLAS

711: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
712: @*/
713: PetscErrorCode  VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
714: {

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

734:   PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
735:   (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
736:   PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
737:   PetscObjectStateIncrease((PetscObject)z);
738:   return(0);
739: }

743: /*@
744:    VecAYPX - Computes y = x + alpha y.

746:    Logically Collective on Vec

748:    Input Parameters:
749: +  alpha - the scalar
750: -  x, y  - the vectors

752:    Output Parameter:
753: .  y - output vector

755:    Level: intermediate

757:    Notes: x and y MUST be different vectors

759:    Concepts: vector^BLAS
760:    Concepts: BLAS

762: .seealso: VecAXPY(), VecWAXPY()
763: @*/
764: PetscErrorCode  VecAYPX(Vec y,PetscScalar alpha,Vec x)
765: {

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

776:   PetscLogEventBegin(VEC_AYPX,x,y,0,0);
777:    (*y->ops->aypx)(y,alpha,x);
778:   PetscLogEventEnd(VEC_AYPX,x,y,0,0);
779:   PetscObjectStateIncrease((PetscObject)y);
780:   return(0);
781: }


786: /*@
787:    VecWAXPY - Computes w = alpha x + y.

789:    Logically Collective on Vec

791:    Input Parameters:
792: +  alpha - the scalar
793: -  x, y  - the vectors

795:    Output Parameter:
796: .  w - the result

798:    Level: intermediate

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

802:    Concepts: vector^BLAS
803:    Concepts: BLAS

805: .seealso: VecAXPY(), VecAYPX(), VecAXPBY()
806: @*/
807: PetscErrorCode  VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
808: {

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

826:   PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
827:    (*w->ops->waxpy)(w,alpha,x,y);
828:   PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
829:   PetscObjectStateIncrease((PetscObject)w);
830:   return(0);
831: }


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

839:    Not Collective

841:    Input Parameters:
842: +  x - vector to insert in
843: .  ni - number of elements to add
844: .  ix - indices where to add
845: .  y - array of values
846: -  iora - either INSERT_VALUES or ADD_VALUES, where
847:    ADD_VALUES adds values to any existing entries, and
848:    INSERT_VALUES replaces existing entries with new values

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

853:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
854:    options cannot be mixed without intervening calls to the assembly
855:    routines.

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

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

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

868:    Level: beginner

870:    Concepts: vector^setting values

872: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
873:            VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
874: @*/
875: PetscErrorCode  VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
876: {

881:   if (!ni) return(0);
885:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
886:   (*x->ops->setvalues)(x,ni,ix,y,iora);
887:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
888:   PetscObjectStateIncrease((PetscObject)x);
889:   return(0);
890: }

894: /*@
895:    VecGetValues - Gets values from certain locations of a vector. Currently
896:           can only get values on the same processor

898:     Not Collective

900:    Input Parameters:
901: +  x - vector to get values from
902: .  ni - number of elements to get
903: -  ix - indices where to get them from (in global 1d numbering)

905:    Output Parameter:
906: .   y - array of values

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

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

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

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

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

921:    Level: beginner

923:    Concepts: vector^getting values

925: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecGetValuesLocal(),
926:            VecGetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecSetValues()
927: @*/
928: PetscErrorCode  VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
929: {

934:   if (!ni) return(0);
938:   (*x->ops->getvalues)(x,ni,ix,y);
939:   return(0);
940: }

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

947:    Not Collective

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

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

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

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

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

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

976:    Level: intermediate

978:    Concepts: vector^setting values blocked

980: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
981:            VecSetValues()
982: @*/
983: PetscErrorCode  VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
984: {

992:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
993:   (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
994:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
995:   PetscObjectStateIncrease((PetscObject)x);
996:   return(0);
997: }


1002: /*@
1003:    VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1004:    using a local ordering of the nodes.

1006:    Not Collective

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

1017:    Level: intermediate

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

1022:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
1023:    options cannot be mixed without intervening calls to the assembly
1024:    routines.

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

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

1031:    Concepts: vector^setting values with local numbering

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

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

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

1069: /*@
1070:    VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1071:    using a local ordering of the nodes.

1073:    Not Collective

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

1084:    Level: intermediate

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

1090:    Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1091:    options cannot be mixed without intervening calls to the assembly
1092:    routines.

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

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


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

1102: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
1103:            VecSetLocalToGlobalMapping()
1104: @*/
1105: PetscErrorCode  VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1106: {
1108:   PetscInt       lixp[128],*lix = lixp;

1115:   if (ni > 128) {
1116:     PetscMalloc1(ni,&lix);
1117:   }

1119:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1120:   ISLocalToGlobalMappingApplyBlock(x->map->mapping,ni,(PetscInt*)ix,lix);
1121:   (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1122:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1123:   if (ni > 128) {
1124:     PetscFree(lix);
1125:   }
1126:   PetscObjectStateIncrease((PetscObject)x);
1127:   return(0);
1128: }

1132: /*@
1133:    VecMTDot - Computes indefinite vector multiple dot products.
1134:    That is, it does NOT use the complex conjugate.

1136:    Collective on Vec

1138:    Input Parameters:
1139: +  x - one vector
1140: .  nv - number of vectors
1141: -  y - array of vectors.  Note that vectors are pointers

1143:    Output Parameter:
1144: .  val - array of the dot products

1146:    Notes for Users of Complex Numbers:
1147:    For complex vectors, VecMTDot() computes the indefinite form
1148: $      val = (x,y) = y^T x,
1149:    where y^T denotes the transpose of y.

1151:    Use VecMDot() for the inner product
1152: $      val = (x,y) = y^H x,
1153:    where y^H denotes the conjugate transpose of y.

1155:    Level: intermediate

1157:    Concepts: inner product^multiple
1158:    Concepts: vector^multiple inner products

1160: .seealso: VecMDot(), VecTDot()
1161: @*/
1162: PetscErrorCode  VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1163: {


1176:   PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1177:   (*x->ops->mtdot)(x,nv,y,val);
1178:   PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1179:   return(0);
1180: }

1184: /*@
1185:    VecMDot - Computes vector multiple dot products.

1187:    Collective on Vec

1189:    Input Parameters:
1190: +  x - one vector
1191: .  nv - number of vectors
1192: -  y - array of vectors.

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

1197:    Notes for Users of Complex Numbers:
1198:    For complex vectors, VecMDot() computes
1199: $     val = (x,y) = y^H x,
1200:    where y^H denotes the conjugate transpose of y.

1202:    Use VecMTDot() for the indefinite form
1203: $     val = (x,y) = y^T x,
1204:    where y^T denotes the transpose of y.

1206:    Level: intermediate

1208:    Concepts: inner product^multiple
1209:    Concepts: vector^multiple inner products

1211: .seealso: VecMTDot(), VecDot()
1212: @*/
1213: PetscErrorCode  VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1214: {

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

1229:   PetscLogEventBarrierBegin(VEC_MDotBarrier,x,*y,0,0,PetscObjectComm((PetscObject)x));
1230:   (*x->ops->mdot)(x,nv,y,val);
1231:   PetscLogEventBarrierEnd(VEC_MDotBarrier,x,*y,0,0,PetscObjectComm((PetscObject)x));
1232:   return(0);
1233: }

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

1240:    Logically Collective on Vec

1242:    Input Parameters:
1243: +  nv - number of scalars and x-vectors
1244: .  alpha - array of scalars
1245: .  y - one vector
1246: -  x - array of vectors

1248:    Level: intermediate

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

1252:    Concepts: BLAS

1254: .seealso: VecAXPY(), VecWAXPY(), VecAYPX()
1255: @*/
1256: PetscErrorCode  VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1257: {
1259:   PetscInt       i;

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

1274:   PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1275:   (*y->ops->maxpy)(y,nv,alpha,x);
1276:   PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1277:   PetscObjectStateIncrease((PetscObject)y);
1278:   return(0);
1279: }

1283: /*@
1284:    VecGetSubVector - Gets a vector representing part of another vector

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

1288:    Input Arguments:
1289: + X - vector from which to extract a subvector
1290: - is - index set representing portion of X to extract

1292:    Output Arguments:
1293: . Y - subvector corresponding to is

1295:    Level: advanced

1297:    Notes:
1298:    The subvector Y should be returned with VecRestoreSubVector().

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

1303: .seealso: MatGetSubMatrix()
1304: @*/
1305: PetscErrorCode  VecGetSubVector(Vec X,IS is,Vec *Y)
1306: {
1307:   PetscErrorCode   ierr;
1308:   Vec              Z;
1309:   PetscObjectState state;

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

1363: /*@
1364:    VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()

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

1368:    Input Arguments:
1369: + X - vector from which subvector was obtained
1370: . is - index set representing the subset of X
1371: - Y - subvector being restored

1373:    Level: advanced

1375: .seealso: VecGetSubVector()
1376: @*/
1377: PetscErrorCode  VecRestoreSubVector(Vec X,IS is,Vec *Y)
1378: {

1386:   if (X->ops->restoresubvector) {
1387:     (*X->ops->restoresubvector)(X,is,Y);
1388:   } else {
1389:     PetscObjectState savedstate=0,newstate;
1390:     PetscBool valid;
1391:     PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,savedstate,valid);
1392:     PetscObjectStateGet((PetscObject)*Y,&newstate);
1393:     if (valid && savedstate < newstate) {
1394:       /* We might need to copy entries back, first check whether we have no-copy view */
1395:       PetscInt  gstart,gend,start;
1396:       PetscBool contiguous,gcontiguous;
1397:       VecGetOwnershipRange(X,&gstart,&gend);
1398:       ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1399:       MPI_Allreduce(&contiguous,&gcontiguous,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1400:       if (!gcontiguous) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Unhandled case, values have been changed and need to be copied back into X");
1401:     }
1402:     VecDestroy(Y);
1403:   }
1404:   return(0);
1405: }

1407: /*@C
1408:    VecGetLocalVectorRead - Maps the local portion of a vector into a
1409:    sequential vector.  This function is similar to VecGetArray which
1410:    maps the local portion into a raw pointer.
1411:    */
1414: PetscErrorCode VecGetLocalVectorRead(Vec v,Vec *w)
1415: {
1417:   PetscScalar    *a;
1418:   PetscInt       m1,m2;

1424:   VecGetLocalSize(v,&m1);
1425:   VecGetLocalSize(*w,&m2);
1426:   if (m1 != m2) {
1427:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vectors of different local sizes.");
1428:   }
1429:   if (v->ops->getlocalvectorread) {
1430:     (*v->ops->getlocalvectorread)(v,w);
1431:   } else {
1432:     VecGetArrayRead(v,(const PetscScalar**)&a);
1433:     VecPlaceArray(*w,a);
1434:   }
1435:   return(0);
1436: }

1438: /*@C
1439:    VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1440:    previously mapped into a sequential vector using
1441:    VecGetLocalVectorRead.  This function is similar to VecGetArray which
1442:    maps the local portion into a raw pointer.
1443:    */
1446: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec *w)
1447: {
1449:   PetscScalar    *a;

1455:   if (v->ops->restorelocalvectorread) {
1456:     (*v->ops->restorelocalvectorread)(v,w);
1457:   } else {
1458:     VecGetArrayRead(*w,(const PetscScalar**)&a);
1459:     VecRestoreArrayRead(v,(const PetscScalar**)&a);
1460:     VecResetArray(*w);
1461:   }
1462:   return(0);
1463: }

1465: /*@C
1466:    VecGetLocalVector - Maps the local portion of a vector into a
1467:    sequential vector.  This function is similar to VecGetArray which
1468:    maps the local portion into a raw pointer.
1469:    */
1472: PetscErrorCode VecGetLocalVector(Vec v,Vec *w)
1473: {
1475:   PetscScalar    *a;
1476:   PetscInt       m1,m2;

1482:   VecGetLocalSize(v,&m1);
1483:   VecGetLocalSize(*w,&m2);
1484:   if (m1 != m2) {
1485:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vectors of different local sizes.");
1486:   }
1487:   if (v->ops->getlocalvector) {
1488:     (*v->ops->getlocalvector)(v,w);
1489:   } else {
1490:     VecGetArray(v,&a);
1491:     VecPlaceArray(*w,a);
1492:   }
1493:   return(0);
1494: }

1496: /*@C
1497:    VecRestoreLocalVector - Unmaps the local portion of a vector
1498:    previously mapped into a sequential vector using
1499:    VecGetLocalVectorRead.  This function is similar to VecGetArray which
1500:    maps the local portion into a raw pointer.
1501:    */
1504: PetscErrorCode VecRestoreLocalVector(Vec v,Vec *w)
1505: {
1507:   PetscScalar    *a;

1513:   if (v->ops->restorelocalvector) {
1514:     (*v->ops->restorelocalvector)(v,w);
1515:   } else {
1516:     VecGetArray(*w,&a);
1517:     VecRestoreArray(v,&a);
1518:     VecResetArray(*w);
1519:   }
1520:   return(0);
1521: }

1525: /*@C
1526:    VecGetArray - Returns a pointer to a contiguous array that contains this
1527:    processor's portion of the vector data. For the standard PETSc
1528:    vectors, VecGetArray() returns a pointer to the local data array and
1529:    does not use any copies. If the underlying vector data is not stored
1530:    in a contiquous array this routine will copy the data to a contiquous
1531:    array and return a pointer to that. You MUST call VecRestoreArray()
1532:    when you no longer need access to the array.

1534:    Logically Collective on Vec

1536:    Input Parameter:
1537: .  x - the vector

1539:    Output Parameter:
1540: .  a - location to put pointer to the array

1542:    Fortran Note:
1543:    This routine is used differently from Fortran 77
1544: $    Vec         x
1545: $    PetscScalar x_array(1)
1546: $    PetscOffset i_x
1547: $    PetscErrorCode ierr
1548: $       call VecGetArray(x,x_array,i_x,ierr)
1549: $
1550: $   Access first local entry in vector with
1551: $      value = x_array(i_x + 1)
1552: $
1553: $      ...... other code
1554: $       call VecRestoreArray(x,x_array,i_x,ierr)
1555:    For Fortran 90 see VecGetArrayF90()

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

1560:    Level: beginner

1562:    Concepts: vector^accessing local values

1564: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(), VecGetArray2d()
1565: @*/
1566: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1567: {

1572:   if (x->petscnative) {
1573: #if defined(PETSC_HAVE_CUSP)
1574:     if (x->valid_GPU_array == PETSC_CUSP_GPU) {
1575:       VecCUSPCopyFromGPU(x);
1576:     } else if (x->valid_GPU_array == PETSC_CUSP_UNALLOCATED) {
1577:       VecCUSPAllocateCheckHost(x);
1578:     }
1579: #endif
1580: #if defined(PETSC_HAVE_VIENNACL)
1581:     if (x->valid_GPU_array == PETSC_VIENNACL_GPU) {
1582:       VecViennaCLCopyFromGPU(x);
1583:     }
1584: #endif
1585:     *a = *((PetscScalar**)x->data);
1586:   } else {
1587:     (*x->ops->getarray)(x,a);
1588:   }
1589:   return(0);
1590: }

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

1597:    Not Collective

1599:    Input Parameters:
1600: .  x - the vector

1602:    Output Parameter:
1603: .  a - the array

1605:    Level: beginner

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

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

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

1616: .seealso: VecGetArray(), VecRestoreArray()
1617: @*/
1618: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1619: {

1624:   if (x->petscnative) {
1625: #if defined(PETSC_HAVE_CUSP)
1626:     if (x->valid_GPU_array == PETSC_CUSP_GPU) {
1627:       VecCUSPCopyFromGPU(x);
1628:     }
1629: #endif
1630: #if defined(PETSC_HAVE_VIENNACL)
1631:     if (x->valid_GPU_array == PETSC_VIENNACL_GPU) {
1632:       VecViennaCLCopyFromGPU(x);
1633:     }
1634: #endif
1635:     *a = *((PetscScalar **)x->data);
1636:   } else if (x->ops->getarrayread) {
1637:     (*x->ops->getarrayread)(x,a);
1638:   } else {
1639:     (*x->ops->getarray)(x,(PetscScalar**)a);
1640:   }
1641:   return(0);
1642: }

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

1651:    Logically Collective on Vec

1653:    Input Parameter:
1654: +  x - the vectors
1655: -  n - the number of vectors

1657:    Output Parameter:
1658: .  a - location to put pointer to the array

1660:    Fortran Note:
1661:    This routine is not supported in Fortran.

1663:    Level: intermediate

1665: .seealso: VecGetArray(), VecRestoreArrays()
1666: @*/
1667: PetscErrorCode  VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1668: {
1670:   PetscInt       i;
1671:   PetscScalar    **q;

1677:   if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1678:   PetscMalloc1(n,&q);
1679:   for (i=0; i<n; ++i) {
1680:     VecGetArray(x[i],&q[i]);
1681:   }
1682:   *a = q;
1683:   return(0);
1684: }

1688: /*@C
1689:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1690:    has been called.

1692:    Logically Collective on Vec

1694:    Input Parameters:
1695: +  x - the vector
1696: .  n - the number of vectors
1697: -  a - location of pointer to arrays obtained from VecGetArrays()

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

1705:    Fortran Note:
1706:    This routine is not supported in Fortran.

1708:    Level: intermediate

1710: .seealso: VecGetArrays(), VecRestoreArray()
1711: @*/
1712: PetscErrorCode  VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1713: {
1715:   PetscInt       i;
1716:   PetscScalar    **q = *a;


1723:   for (i=0; i<n; ++i) {
1724:     VecRestoreArray(x[i],&q[i]);
1725:   }
1726:   PetscFree(q);
1727:   return(0);
1728: }

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

1735:    Logically Collective on Vec

1737:    Input Parameters:
1738: +  x - the vector
1739: -  a - location of pointer to array obtained from VecGetArray()

1741:    Level: beginner

1743:    Notes:
1744:    For regular PETSc vectors this routine does not involve any copies. For
1745:    any special vectors that do not store local vector data in a contiguous
1746:    array, this routine will copy the data back into the underlying
1747:    vector data structure from the array obtained with VecGetArray().

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

1753:    Fortran Note:
1754:    This routine is used differently from Fortran 77
1755: $    Vec         x
1756: $    PetscScalar x_array(1)
1757: $    PetscOffset i_x
1758: $    PetscErrorCode ierr
1759: $       call VecGetArray(x,x_array,i_x,ierr)
1760: $
1761: $   Access first local entry in vector with
1762: $      value = x_array(i_x + 1)
1763: $
1764: $      ...... other code
1765: $       call VecRestoreArray(x,x_array,i_x,ierr)

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

1771: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(), VecRestoreArray2d()
1772: @*/
1773: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1774: {

1779:   if (x->petscnative) {
1780: #if defined(PETSC_HAVE_CUSP)
1781:     x->valid_GPU_array = PETSC_CUSP_CPU;
1782: #endif
1783: #if defined(PETSC_HAVE_VIENNACL)
1784:     x->valid_GPU_array = PETSC_VIENNACL_CPU;
1785: #endif
1786:   } else {
1787:     (*x->ops->restorearray)(x,a);
1788:   }
1789:   if (a) *a = NULL;
1790:   PetscObjectStateIncrease((PetscObject)x);
1791:   return(0);
1792: }

1796: /*@C
1797:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

1799:    Not Collective

1801:    Input Parameters:
1802: +  vec - the vector
1803: -  array - the array

1805:    Level: beginner

1807: .seealso: VecGetArray(), VecRestoreArray()
1808: @*/
1809: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1810: {

1815:   if (x->petscnative) {
1816: #if defined(PETSC_HAVE_VIENNACL)
1817:     x->valid_GPU_array = PETSC_VIENNACL_CPU;
1818: #endif
1819:   } else if (x->ops->restorearrayread) {
1820:     (*x->ops->restorearrayread)(x,a);
1821:   } else {
1822:     (*x->ops->restorearray)(x,(PetscScalar**)a);
1823:   }
1824:   if (a) *a = NULL;
1825:   return(0);
1826: }

1830: /*@
1831:    VecPlaceArray - Allows one to replace the array in a vector with an
1832:    array provided by the user. This is useful to avoid copying an array
1833:    into a vector.

1835:    Not Collective

1837:    Input Parameters:
1838: +  vec - the vector
1839: -  array - the array

1841:    Notes:
1842:    You can return to the original array with a call to VecResetArray()

1844:    Level: developer

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

1848: @*/
1849: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
1850: {

1857:   if (vec->ops->placearray) {
1858:     (*vec->ops->placearray)(vec,array);
1859:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
1860:   PetscObjectStateIncrease((PetscObject)vec);
1861:   return(0);
1862: }

1866: /*@C
1867:    VecReplaceArray - Allows one to replace the array in a vector with an
1868:    array provided by the user. This is useful to avoid copying an array
1869:    into a vector.

1871:    Not Collective

1873:    Input Parameters:
1874: +  vec - the vector
1875: -  array - the array

1877:    Notes:
1878:    This permanently replaces the array and frees the memory associated
1879:    with the old array.

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

1884:    Not supported from Fortran

1886:    Level: developer

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

1890: @*/
1891: PetscErrorCode  VecReplaceArray(Vec vec,const PetscScalar array[])
1892: {

1898:   if (vec->ops->replacearray) {
1899:     (*vec->ops->replacearray)(vec,array);
1900:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
1901:   PetscObjectStateIncrease((PetscObject)vec);
1902:   return(0);
1903: }

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

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

1912:     Collective on Vec

1914:     Input Parameters:
1915: +   x - a vector to mimic
1916: -   n - the number of vectors to obtain

1918:     Output Parameters:
1919: +   y - Fortran90 pointer to the array of vectors
1920: -   ierr - error code

1922:     Example of Usage:
1923: .vb
1924:     Vec x
1925:     Vec, pointer :: y(:)
1926:     ....
1927:     call VecDuplicateVecsF90(x,2,y,ierr)
1928:     call VecSet(y(2),alpha,ierr)
1929:     call VecSet(y(2),alpha,ierr)
1930:     ....
1931:     call VecDestroyVecsF90(2,y,ierr)
1932: .ve

1934:     Notes:
1935:     Not yet supported for all F90 compilers

1937:     Use VecDestroyVecsF90() to free the space.

1939:     Level: beginner

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

1943: M*/

1945: /*MC
1946:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
1947:     VecGetArrayF90().

1949:     Synopsis:
1950:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

1952:     Logically Collective on Vec

1954:     Input Parameters:
1955: +   x - vector
1956: -   xx_v - the Fortran90 pointer to the array

1958:     Output Parameter:
1959: .   ierr - error code

1961:     Example of Usage:
1962: .vb
1963:     PetscScalar, pointer :: xx_v(:)
1964:     ....
1965:     call VecGetArrayF90(x,xx_v,ierr)
1966:     a = xx_v(3)
1967:     call VecRestoreArrayF90(x,xx_v,ierr)
1968: .ve

1970:     Level: beginner

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

1974: M*/

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

1979:     Synopsis:
1980:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

1982:     Collective on Vec

1984:     Input Parameters:
1985: +   n - the number of vectors previously obtained
1986: -   x - pointer to array of vector pointers

1988:     Output Parameter:
1989: .   ierr - error code

1991:     Notes:
1992:     Not yet supported for all F90 compilers

1994:     Level: beginner

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

1998: M*/

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

2006:     Synopsis:
2007:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2009:     Logically Collective on Vec

2011:     Input Parameter:
2012: .   x - vector

2014:     Output Parameters:
2015: +   xx_v - the Fortran90 pointer to the array
2016: -   ierr - error code

2018:     Example of Usage:
2019: .vb
2020:     PetscScalar, pointer :: xx_v(:)
2021:     ....
2022:     call VecGetArrayF90(x,xx_v,ierr)
2023:     a = xx_v(3)
2024:     call VecRestoreArrayF90(x,xx_v,ierr)
2025: .ve

2027:     Level: beginner

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

2031: M*/


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

2041:    Logically Collective

2043:    Input Parameter:
2044: +  x - the vector
2045: .  m - first dimension of two dimensional array
2046: .  n - second dimension of two dimensional array
2047: .  mstart - first index you will use in first coordinate direction (often 0)
2048: -  nstart - first index in the second coordinate direction (often 0)

2050:    Output Parameter:
2051: .  a - location to put pointer to the array

2053:    Level: developer

2055:   Notes:
2056:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2057:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2058:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2059:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

2065: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2066:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2067:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2068: @*/
2069: PetscErrorCode  VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2070: {
2072:   PetscInt       i,N;
2073:   PetscScalar    *aa;

2079:   VecGetLocalSize(x,&N);
2080:   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);
2081:   VecGetArray(x,&aa);

2083:   PetscMalloc1(m,a);
2084:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2085:   *a -= mstart;
2086:   return(0);
2087: }

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

2094:    Logically Collective

2096:    Input Parameters:
2097: +  x - the vector
2098: .  m - first dimension of two dimensional array
2099: .  n - second dimension of the two dimensional array
2100: .  mstart - first index you will use in first coordinate direction (often 0)
2101: .  nstart - first index in the second coordinate direction (often 0)
2102: -  a - location of pointer to array obtained from VecGetArray2d()

2104:    Level: developer

2106:    Notes:
2107:    For regular PETSc vectors this routine does not involve any copies. For
2108:    any special vectors that do not store local vector data in a contiguous
2109:    array, this routine will copy the data back into the underlying
2110:    vector data structure from the array obtained with VecGetArray().

2112:    This routine actually zeros out the a pointer.

2114: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2115:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2116:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2117: @*/
2118: PetscErrorCode  VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2119: {
2121:   void           *dummy;

2127:   dummy = (void*)(*a + mstart);
2128:   PetscFree(dummy);
2129:   VecRestoreArray(x,NULL);
2130:   return(0);
2131: }

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

2140:    Logically Collective

2142:    Input Parameter:
2143: +  x - the vector
2144: .  m - first dimension of two dimensional array
2145: -  mstart - first index you will use in first coordinate direction (often 0)

2147:    Output Parameter:
2148: .  a - location to put pointer to the array

2150:    Level: developer

2152:   Notes:
2153:    For a vector obtained from DMCreateLocalVector() mstart are likely
2154:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2155:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2159: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2160:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2161:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2162: @*/
2163: PetscErrorCode  VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2164: {
2166:   PetscInt       N;

2172:   VecGetLocalSize(x,&N);
2173:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2174:   VecGetArray(x,a);
2175:   *a  -= mstart;
2176:   return(0);
2177: }

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

2184:    Logically Collective

2186:    Input Parameters:
2187: +  x - the vector
2188: .  m - first dimension of two dimensional array
2189: .  mstart - first index you will use in first coordinate direction (often 0)
2190: -  a - location of pointer to array obtained from VecGetArray21()

2192:    Level: developer

2194:    Notes:
2195:    For regular PETSc vectors this routine does not involve any copies. For
2196:    any special vectors that do not store local vector data in a contiguous
2197:    array, this routine will copy the data back into the underlying
2198:    vector data structure from the array obtained with VecGetArray1d().

2200:    This routine actually zeros out the a pointer.

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

2204: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2205:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2206:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2207: @*/
2208: PetscErrorCode  VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2209: {

2215:   VecRestoreArray(x,NULL);
2216:   return(0);
2217: }


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

2227:    Logically Collective

2229:    Input Parameter:
2230: +  x - the vector
2231: .  m - first dimension of three dimensional array
2232: .  n - second dimension of three dimensional array
2233: .  p - third dimension of three dimensional array
2234: .  mstart - first index you will use in first coordinate direction (often 0)
2235: .  nstart - first index in the second coordinate direction (often 0)
2236: -  pstart - first index in the third coordinate direction (often 0)

2238:    Output Parameter:
2239: .  a - location to put pointer to the array

2241:    Level: developer

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

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

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

2253: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2254:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2255:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2256: @*/
2257: PetscErrorCode  VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2258: {
2260:   PetscInt       i,N,j;
2261:   PetscScalar    *aa,**b;

2267:   VecGetLocalSize(x,&N);
2268:   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);
2269:   VecGetArray(x,&aa);

2271:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
2272:   b    = (PetscScalar**)((*a) + m);
2273:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2274:   for (i=0; i<m; i++)
2275:     for (j=0; j<n; j++)
2276:       b[i*n+j] = aa + i*n*p + j*p - pstart;

2278:   *a -= mstart;
2279:   return(0);
2280: }

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

2287:    Logically Collective

2289:    Input Parameters:
2290: +  x - the vector
2291: .  m - first dimension of three dimensional array
2292: .  n - second dimension of the three dimensional array
2293: .  p - third dimension of the three dimensional array
2294: .  mstart - first index you will use in first coordinate direction (often 0)
2295: .  nstart - first index in the second coordinate direction (often 0)
2296: .  pstart - first index in the third coordinate direction (often 0)
2297: -  a - location of pointer to array obtained from VecGetArray3d()

2299:    Level: developer

2301:    Notes:
2302:    For regular PETSc vectors this routine does not involve any copies. For
2303:    any special vectors that do not store local vector data in a contiguous
2304:    array, this routine will copy the data back into the underlying
2305:    vector data structure from the array obtained with VecGetArray().

2307:    This routine actually zeros out the a pointer.

2309: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2310:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2311:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2312: @*/
2313: PetscErrorCode  VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2314: {
2316:   void           *dummy;

2322:   dummy = (void*)(*a + mstart);
2323:   PetscFree(dummy);
2324:   VecRestoreArray(x,NULL);
2325:   return(0);
2326: }

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

2335:    Logically Collective

2337:    Input Parameter:
2338: +  x - the vector
2339: .  m - first dimension of four dimensional array
2340: .  n - second dimension of four dimensional array
2341: .  p - third dimension of four dimensional array
2342: .  q - fourth dimension of four dimensional array
2343: .  mstart - first index you will use in first coordinate direction (often 0)
2344: .  nstart - first index in the second coordinate direction (often 0)
2345: .  pstart - first index in the third coordinate direction (often 0)
2346: -  qstart - first index in the fourth coordinate direction (often 0)

2348:    Output Parameter:
2349: .  a - location to put pointer to the array

2351:    Level: beginner

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

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

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

2363: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2364:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2365:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2366: @*/
2367: PetscErrorCode  VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2368: {
2370:   PetscInt       i,N,j,k;
2371:   PetscScalar    *aa,***b,**c;

2377:   VecGetLocalSize(x,&N);
2378:   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);
2379:   VecGetArray(x,&aa);

2381:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
2382:   b    = (PetscScalar***)((*a) + m);
2383:   c    = (PetscScalar**)(b + m*n);
2384:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2385:   for (i=0; i<m; i++)
2386:     for (j=0; j<n; j++)
2387:       b[i*n+j] = c + i*n*p + j*p - pstart;
2388:   for (i=0; i<m; i++)
2389:     for (j=0; j<n; j++)
2390:       for (k=0; k<p; k++)
2391:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
2392:   *a -= mstart;
2393:   return(0);
2394: }

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

2401:    Logically Collective

2403:    Input Parameters:
2404: +  x - the vector
2405: .  m - first dimension of four dimensional array
2406: .  n - second dimension of the four dimensional array
2407: .  p - third dimension of the four dimensional array
2408: .  q - fourth dimension of the four dimensional array
2409: .  mstart - first index you will use in first coordinate direction (often 0)
2410: .  nstart - first index in the second coordinate direction (often 0)
2411: .  pstart - first index in the third coordinate direction (often 0)
2412: .  qstart - first index in the fourth coordinate direction (often 0)
2413: -  a - location of pointer to array obtained from VecGetArray4d()

2415:    Level: beginner

2417:    Notes:
2418:    For regular PETSc vectors this routine does not involve any copies. For
2419:    any special vectors that do not store local vector data in a contiguous
2420:    array, this routine will copy the data back into the underlying
2421:    vector data structure from the array obtained with VecGetArray().

2423:    This routine actually zeros out the a pointer.

2425: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2426:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2427:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2428: @*/
2429: PetscErrorCode  VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2430: {
2432:   void           *dummy;

2438:   dummy = (void*)(*a + mstart);
2439:   PetscFree(dummy);
2440:   VecRestoreArray(x,NULL);
2441:   return(0);
2442: }