Actual source code: rvector.c

petsc-master 2017-04-24
Report Typos and Errors

  2: /*
  3:      Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
  4:    These are the vector functions the user calls.
  5: */
  6:  #include <petsc/private/vecimpl.h>
  7: static PetscInt VecGetSubVectorSavedStateId = -1;

  9: PETSC_EXTERN PetscErrorCode VecValidValues(Vec vec,PetscInt argnum,PetscBool begin)
 10: {
 11: #if defined(PETSC_USE_DEBUG)
 12:   PetscErrorCode    ierr;
 13:   PetscInt          n,i;
 14:   const PetscScalar *x;

 17: #if defined(PETSC_HAVE_CUSP)
 18:   if ((vec->petscnative || vec->ops->getarray) && (vec->valid_GPU_array == PETSC_CUSP_CPU || vec->valid_GPU_array == PETSC_CUSP_BOTH)) {
 19: #elif defined(PETSC_HAVE_VECCUDA)
 20:   if ((vec->petscnative || vec->ops->getarray) && (vec->valid_GPU_array == PETSC_CUDA_CPU || vec->valid_GPU_array == PETSC_CUDA_BOTH)) {
 21: #elif defined(PETSC_HAVE_VIENNACL)
 22:   if ((vec->petscnative || vec->ops->getarray) && (vec->valid_GPU_array == PETSC_VIENNACL_CPU || vec->valid_GPU_array == PETSC_VIENNACL_BOTH)) {
 23: #else
 24:   if (vec->petscnative || vec->ops->getarray) {
 25: #endif
 26:     VecGetLocalSize(vec,&n);
 27:     VecGetArrayRead(vec,&x);
 28:     for (i=0; i<n; i++) {
 29:       if (begin) {
 30:         if (PetscIsInfOrNanScalar(x[i])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FP,"Vec entry at local location %D is not-a-number or infinite at beginning of function: Parameter number %D",i,argnum);
 31:       } else {
 32:         if (PetscIsInfOrNanScalar(x[i])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FP,"Vec entry at local location %D is not-a-number or infinite at end of function: Parameter number %D",i,argnum);
 33:       }
 34:     }
 35:     VecRestoreArrayRead(vec,&x);
 36:   }
 37:   return(0);
 38: #else
 39:   return 0;
 40: #endif
 41: }

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

 46:    Logically Collective on Vec

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

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

 54:    Level: advanced

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

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

 72:   VecCheckSameSize(x,1,y,2);
 73:   (*x->ops->maxpointwisedivide)(x,y,max);
 74:   return(0);
 75: }

 77: /*@
 78:    VecDot - Computes the vector dot product.

 80:    Collective on Vec

 82:    Input Parameters:
 83: .  x, y - the vectors

 85:    Output Parameter:
 86: .  val - the dot product

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

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

100:    Use VecTDot() for the indefinite form
101: $     val = (x,y) = y^T x,
102:    where y^T denotes the transpose of y.

104:    Level: intermediate

106:    Concepts: inner product
107:    Concepts: vector^inner product

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

122:   VecCheckSameSize(x,1,y,2);

124:   PetscLogEventBarrierBegin(VEC_DotBarrier,x,y,0,0,PetscObjectComm((PetscObject)x));
125:   (*x->ops->dot)(x,y,val);
126:   PetscLogEventBarrierEnd(VEC_DotBarrier,x,y,0,0,PetscObjectComm((PetscObject)x));
127:   return(0);
128: }

130: /*@
131:    VecDotRealPart - Computes the real part of the vector dot product.

133:    Collective on Vec

135:    Input Parameters:
136: .  x, y - the vectors

138:    Output Parameter:
139: .  val - the real part of the dot product;

141:    Performance Issues:
142: $    per-processor memory bandwidth
143: $    interprocessor latency
144: $    work load inbalance that causes certain processes to arrive much earlier than others

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

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

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

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

156:    Level: intermediate

158:    Concepts: inner product
159:    Concepts: vector^inner product

161: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDot(), VecDotNorm2()
162: @*/
163: PetscErrorCode  VecDotRealPart(Vec x,Vec y,PetscReal *val)
164: {
166:   PetscScalar    fdot;

169:   VecDot(x,y,&fdot);
170:   *val = PetscRealPart(fdot);
171:   return(0);
172: }

174: /*@
175:    VecNorm  - Computes the vector norm.

177:    Collective on Vec

179:    Input Parameters:
180: +  x - the vector
181: -  type - one of NORM_1, NORM_2, NORM_INFINITY.  Also available
182:           NORM_1_AND_2, which computes both norms and stores them
183:           in a two element array.

185:    Output Parameter:
186: .  val - the norm

188:    Notes:
189: $     NORM_1 denotes sum_i |x_i|
190: $     NORM_2 denotes sqrt(sum_i |x_i|^2)
191: $     NORM_INFINITY denotes max_i |x_i|

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

198:    Level: intermediate

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

205:    Concepts: norm
206:    Concepts: vector^norm

208: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNormAvailable(),
209:           VecNormBegin(), VecNormEnd()

211: @*/
212: PetscErrorCode  VecNorm(Vec x,NormType type,PetscReal *val)
213: {
214:   PetscBool      flg;


222:   /*
223:    * Cached data?
224:    */
225:   if (type!=NORM_1_AND_2) {
226:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,flg);
227:     if (flg) return(0);
228:   }
229:   PetscLogEventBarrierBegin(VEC_NormBarrier,x,0,0,0,PetscObjectComm((PetscObject)x));
230:   (*x->ops->norm)(x,type,val);
231:   PetscLogEventBarrierEnd(VEC_NormBarrier,x,0,0,0,PetscObjectComm((PetscObject)x));

233:   if (type!=NORM_1_AND_2) {
234:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[type],*val);
235:   }
236:   return(0);
237: }

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

242:    Not Collective

244:    Input Parameters:
245: +  x - the vector
246: -  type - one of NORM_1, NORM_2, NORM_INFINITY.  Also available
247:           NORM_1_AND_2, which computes both norms and stores them
248:           in a two element array.

250:    Output Parameter:
251: +  available - PETSC_TRUE if the val returned is valid
252: -  val - the norm

254:    Notes:
255: $     NORM_1 denotes sum_i |x_i|
256: $     NORM_2 denotes sqrt(sum_i (x_i)^2)
257: $     NORM_INFINITY denotes max_i |x_i|

259:    Level: intermediate

261:    Performance Issues:
262: $    per-processor memory bandwidth
263: $    interprocessor latency
264: $    work load inbalance that causes certain processes to arrive much earlier than others

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

271:    Concepts: norm
272:    Concepts: vector^norm

274: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNorm()
275:           VecNormBegin(), VecNormEnd()

277: @*/
278: PetscErrorCode  VecNormAvailable(Vec x,NormType type,PetscBool  *available,PetscReal *val)
279: {


287:   *available = PETSC_FALSE;
288:   if (type!=NORM_1_AND_2) {
289:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,*available);
290:   }
291:   return(0);
292: }

294: /*@
295:    VecNormalize - Normalizes a vector by 2-norm.

297:    Collective on Vec

299:    Input Parameters:
300: +  x - the vector

302:    Output Parameter:
303: .  x - the normalized vector
304: -  val - the vector norm before normalization

306:    Level: intermediate

308:    Concepts: vector^normalizing
309:    Concepts: normalizing^vector

311: @*/
312: PetscErrorCode  VecNormalize(Vec x,PetscReal *val)
313: {
315:   PetscReal      norm;

320:   PetscLogEventBegin(VEC_Normalize,x,0,0,0);
321:   VecNorm(x,NORM_2,&norm);
322:   if (norm == 0.0) {
323:     PetscInfo(x,"Vector of zero norm can not be normalized; Returning only the zero norm\n");
324:   } else if (norm != 1.0) {
325:     PetscScalar tmp = 1.0/norm;
326:     VecScale(x,tmp);
327:   }
328:   if (val) *val = norm;
329:   PetscLogEventEnd(VEC_Normalize,x,0,0,0);
330:   return(0);
331: }

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

336:    Collective on Vec

338:    Input Parameter:
339: .  x - the vector

341:    Output Parameters:
342: +  val - the maximum component
343: -  p - the location of val (pass NULL if you don't want this)

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

348:    Returns the smallest index with the maximum value
349:    Level: intermediate

351:    Concepts: maximum^of vector
352:    Concepts: vector^maximum value

354: .seealso: VecNorm(), VecMin()
355: @*/
356: PetscErrorCode  VecMax(Vec x,PetscInt *p,PetscReal *val)
357: {

364:   PetscLogEventBegin(VEC_Max,x,0,0,0);
365:   (*x->ops->max)(x,p,val);
366:   PetscLogEventEnd(VEC_Max,x,0,0,0);
367:   return(0);
368: }

370: /*@C
371:    VecMin - Determines the minimum vector component and its location.

373:    Collective on Vec

375:    Input Parameters:
376: .  x - the vector

378:    Output Parameter:
379: +  val - the minimum component
380: -  p - the location of val (pass NULL if you don't want this location)

382:    Level: intermediate

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

387:    This returns the smallest index with the minumum value

389:    Concepts: minimum^of vector
390:    Concepts: vector^minimum entry

392: .seealso: VecMax()
393: @*/
394: PetscErrorCode  VecMin(Vec x,PetscInt *p,PetscReal *val)
395: {

402:   PetscLogEventBegin(VEC_Min,x,0,0,0);
403:   (*x->ops->min)(x,p,val);
404:   PetscLogEventEnd(VEC_Min,x,0,0,0);
405:   return(0);
406: }

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

412:    Collective on Vec

414:    Input Parameters:
415: .  x, y - the vectors

417:    Output Parameter:
418: .  val - the dot product

420:    Notes for Users of Complex Numbers:
421:    For complex vectors, VecTDot() computes the indefinite form
422: $     val = (x,y) = y^T x,
423:    where y^T denotes the transpose of y.

425:    Use VecDot() for the inner product
426: $     val = (x,y) = y^H x,
427:    where y^H denotes the conjugate transpose of y.

429:    Level: intermediate

431:    Concepts: inner product^non-Hermitian
432:    Concepts: vector^inner product
433:    Concepts: non-Hermitian inner product

435: .seealso: VecDot(), VecMTDot()
436: @*/
437: PetscErrorCode  VecTDot(Vec x,Vec y,PetscScalar *val)
438: {

448:   VecCheckSameSize(x,1,y,2);

450:   PetscLogEventBegin(VEC_TDot,x,y,0,0);
451:   (*x->ops->tdot)(x,y,val);
452:   PetscLogEventEnd(VEC_TDot,x,y,0,0);
453:   return(0);
454: }

456: /*@
457:    VecScale - Scales a vector.

459:    Not collective on Vec

461:    Input Parameters:
462: +  x - the vector
463: -  alpha - the scalar

465:    Output Parameter:
466: .  x - the scaled vector

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

472:    Level: intermediate

474:    Concepts: vector^scaling
475:    Concepts: scaling^vector

477: @*/
478: PetscErrorCode  VecScale(Vec x, PetscScalar alpha)
479: {
480:   PetscReal      norms[4] = {0.0,0.0,0.0, 0.0};
481:   PetscBool      flgs[4];
483:   PetscInt       i;

488:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
489:   PetscLogEventBegin(VEC_Scale,x,0,0,0);
490:   if (alpha != (PetscScalar)1.0) {
491:     /* get current stashed norms */
492:     for (i=0; i<4; i++) {
493:       PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
494:     }
495:     (*x->ops->scale)(x,alpha);
496:     PetscObjectStateIncrease((PetscObject)x);
497:     /* put the scaled stashed norms back into the Vec */
498:     for (i=0; i<4; i++) {
499:       if (flgs[i]) {
500:         PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],PetscAbsScalar(alpha)*norms[i]);
501:       }
502:     }
503:   }
504:   PetscLogEventEnd(VEC_Scale,x,0,0,0);
505:   return(0);
506: }

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

511:    Logically Collective on Vec

513:    Input Parameters:
514: +  x  - the vector
515: -  alpha - the scalar

517:    Output Parameter:
518: .  x  - the vector

520:    Note:
521:    For a vector of dimension n, VecSet() computes
522: $     x[i] = alpha, for i=1,...,n,
523:    so that all vector entries then equal the identical
524:    scalar value, alpha.  Use the more general routine
525:    VecSetValues() to set different vector entries.

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

530:    Level: beginner

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

534:    Concepts: vector^setting to constant

536: @*/
537: PetscErrorCode  VecSet(Vec x,PetscScalar alpha)
538: {
539:   PetscReal      val;

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

549:   PetscLogEventBegin(VEC_Set,x,0,0,0);
550:   (*x->ops->set)(x,alpha);
551:   PetscLogEventEnd(VEC_Set,x,0,0,0);
552:   PetscObjectStateIncrease((PetscObject)x);

554:   /*  norms can be simply set (if |alpha|*N not too large) */
555:   val  = PetscAbsScalar(alpha);
556:   if (x->map->N == 0) {
557:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],0.0l);
558:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],0.0);
559:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],0.0);
560:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],0.0);
561:   } else if (val > PETSC_MAX_REAL/x->map->N) {
562:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
563:   } else {
564:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],x->map->N * val);
565:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
566:     val  = PetscSqrtReal((PetscReal)x->map->N) * val;
567:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],val);
568:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],val);
569:   }
570:   return(0);
571: }


574: /*@
575:    VecAXPY - Computes y = alpha x + y.

577:    Logically Collective on Vec

579:    Input Parameters:
580: +  alpha - the scalar
581: -  x, y  - the vectors

583:    Output Parameter:
584: .  y - output vector

586:    Level: intermediate

588:    Notes: x and y MUST be different vectors

590:    Concepts: vector^BLAS
591:    Concepts: BLAS

593: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY()
594: @*/
595: PetscErrorCode  VecAXPY(Vec y,PetscScalar alpha,Vec x)
596: {

605:   VecCheckSameSize(x,1,y,3);
606:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
608:   VecLocked(y,1);

610:   VecLockPush(x);
611:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
612:   (*y->ops->axpy)(y,alpha,x);
613:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
614:   VecLockPop(x);
615:   PetscObjectStateIncrease((PetscObject)y);
616:   return(0);
617: }

619: /*@
620:    VecAXPBY - Computes y = alpha x + beta y.

622:    Logically Collective on Vec

624:    Input Parameters:
625: +  alpha,beta - the scalars
626: -  x, y  - the vectors

628:    Output Parameter:
629: .  y - output vector

631:    Level: intermediate

633:    Notes: x and y MUST be different vectors

635:    Concepts: BLAS
636:    Concepts: vector^BLAS

638: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
639: @*/
640: PetscErrorCode  VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
641: {

650:   VecCheckSameSize(x,1,y,4);
651:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");

655:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
656:   (*y->ops->axpby)(y,alpha,beta,x);
657:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
658:   PetscObjectStateIncrease((PetscObject)y);
659:   return(0);
660: }

662: /*@
663:    VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z

665:    Logically Collective on Vec

667:    Input Parameters:
668: +  alpha,beta, gamma - the scalars
669: -  x, y, z  - the vectors

671:    Output Parameter:
672: .  z - output vector

674:    Level: intermediate

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

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

680:    Concepts: BLAS
681:    Concepts: vector^BLAS

683: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
684: @*/
685: PetscErrorCode  VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
686: {

698:   VecCheckSameSize(x,1,y,5);
699:   VecCheckSameSize(x,1,z,6);
700:   if (x == y || x == z) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
701:   if (y == z) SETERRQ(PetscObjectComm((PetscObject)y),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");

706:   PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
707:   (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
708:   PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
709:   PetscObjectStateIncrease((PetscObject)z);
710:   return(0);
711: }

713: /*@
714:    VecAYPX - Computes y = x + alpha y.

716:    Logically Collective on Vec

718:    Input Parameters:
719: +  alpha - the scalar
720: -  x, y  - the vectors

722:    Output Parameter:
723: .  y - output vector

725:    Level: intermediate

727:    Notes: x and y MUST be different vectors

729:    Concepts: vector^BLAS
730:    Concepts: BLAS

732: .seealso: VecAXPY(), VecWAXPY()
733: @*/
734: PetscErrorCode  VecAYPX(Vec y,PetscScalar alpha,Vec x)
735: {

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

748:   PetscLogEventBegin(VEC_AYPX,x,y,0,0);
749:    (*y->ops->aypx)(y,alpha,x);
750:   PetscLogEventEnd(VEC_AYPX,x,y,0,0);
751:   PetscObjectStateIncrease((PetscObject)y);
752:   return(0);
753: }


756: /*@
757:    VecWAXPY - Computes w = alpha x + y.

759:    Logically Collective on Vec

761:    Input Parameters:
762: +  alpha - the scalar
763: -  x, y  - the vectors

765:    Output Parameter:
766: .  w - the result

768:    Level: intermediate

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

772:    Concepts: vector^BLAS
773:    Concepts: BLAS

775: .seealso: VecAXPY(), VecAYPX(), VecAXPBY()
776: @*/
777: PetscErrorCode  VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
778: {

790:   VecCheckSameSize(x,3,y,4);
791:   VecCheckSameSize(x,3,w,1);
792:   if (w == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector y, suggest VecAXPY()");
793:   if (w == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector x, suggest VecAYPX()");

796:   PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
797:    (*w->ops->waxpy)(w,alpha,x,y);
798:   PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
799:   PetscObjectStateIncrease((PetscObject)w);
800:   return(0);
801: }


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

807:    Not Collective

809:    Input Parameters:
810: +  x - vector to insert in
811: .  ni - number of elements to add
812: .  ix - indices where to add
813: .  y - array of values
814: -  iora - either INSERT_VALUES or ADD_VALUES, where
815:    ADD_VALUES adds values to any existing entries, and
816:    INSERT_VALUES replaces existing entries with new values

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

821:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
822:    options cannot be mixed without intervening calls to the assembly
823:    routines.

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

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

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

836:    Level: beginner

838:    Concepts: vector^setting values

840: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
841:            VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
842: @*/
843: PetscErrorCode  VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
844: {

849:   if (!ni) return(0);
853:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
854:   (*x->ops->setvalues)(x,ni,ix,y,iora);
855:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
856:   PetscObjectStateIncrease((PetscObject)x);
857:   return(0);
858: }

860: /*@
861:    VecGetValues - Gets values from certain locations of a vector. Currently
862:           can only get values on the same processor

864:     Not Collective

866:    Input Parameters:
867: +  x - vector to get values from
868: .  ni - number of elements to get
869: -  ix - indices where to get them from (in global 1d numbering)

871:    Output Parameter:
872: .   y - array of values

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

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

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

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

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

887:    Level: beginner

889:    Concepts: vector^getting values

891: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues()
892: @*/
893: PetscErrorCode  VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
894: {

899:   if (!ni) return(0);
903:   (*x->ops->getvalues)(x,ni,ix,y);
904:   return(0);
905: }

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

910:    Not Collective

912:    Input Parameters:
913: +  x - vector to insert in
914: .  ni - number of blocks to add
915: .  ix - indices where to add in block count, rather than element count
916: .  y - array of values
917: -  iora - either INSERT_VALUES or ADD_VALUES, where
918:    ADD_VALUES adds values to any existing entries, and
919:    INSERT_VALUES replaces existing entries with new values

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

925:    Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
926:    options cannot be mixed without intervening calls to the assembly
927:    routines.

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

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

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

939:    Level: intermediate

941:    Concepts: vector^setting values blocked

943: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
944:            VecSetValues()
945: @*/
946: PetscErrorCode  VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
947: {

955:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
956:   (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
957:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
958:   PetscObjectStateIncrease((PetscObject)x);
959:   return(0);
960: }


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

967:    Not Collective

969:    Input Parameters:
970: +  x - vector to insert in
971: .  ni - number of elements to add
972: .  ix - indices where to add
973: .  y - array of values
974: -  iora - either INSERT_VALUES or ADD_VALUES, where
975:    ADD_VALUES adds values to any existing entries, and
976:    INSERT_VALUES replaces existing entries with new values

978:    Level: intermediate

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

983:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
984:    options cannot be mixed without intervening calls to the assembly
985:    routines.

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

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

992:    Concepts: vector^setting values with local numbering

994: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetLocalToGlobalMapping(),
995:            VecSetValuesBlockedLocal()
996: @*/
997: PetscErrorCode  VecSetValuesLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
998: {
1000:   PetscInt       lixp[128],*lix = lixp;

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

1009:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1010:   if (!x->ops->setvalueslocal) {
1011:     if (!x->map->mapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMapping()");
1012:     if (ni > 128) {
1013:       PetscMalloc1(ni,&lix);
1014:     }
1015:     ISLocalToGlobalMappingApply(x->map->mapping,ni,(PetscInt*)ix,lix);
1016:     (*x->ops->setvalues)(x,ni,lix,y,iora);
1017:     if (ni > 128) {
1018:       PetscFree(lix);
1019:     }
1020:   } else {
1021:     (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
1022:   }
1023:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1024:   PetscObjectStateIncrease((PetscObject)x);
1025:   return(0);
1026: }

1028: /*@
1029:    VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1030:    using a local ordering of the nodes.

1032:    Not Collective

1034:    Input Parameters:
1035: +  x - vector to insert in
1036: .  ni - number of blocks to add
1037: .  ix - indices where to add in block count, not element count
1038: .  y - array of values
1039: -  iora - either INSERT_VALUES or ADD_VALUES, where
1040:    ADD_VALUES adds values to any existing entries, and
1041:    INSERT_VALUES replaces existing entries with new values

1043:    Level: intermediate

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

1049:    Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1050:    options cannot be mixed without intervening calls to the assembly
1051:    routines.

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

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


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

1061: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
1062:            VecSetLocalToGlobalMapping()
1063: @*/
1064: PetscErrorCode  VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1065: {
1067:   PetscInt       lixp[128],*lix = lixp;

1074:   if (ni > 128) {
1075:     PetscMalloc1(ni,&lix);
1076:   }

1078:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1079:   ISLocalToGlobalMappingApplyBlock(x->map->mapping,ni,(PetscInt*)ix,lix);
1080:   (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1081:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1082:   if (ni > 128) {
1083:     PetscFree(lix);
1084:   }
1085:   PetscObjectStateIncrease((PetscObject)x);
1086:   return(0);
1087: }

1089: /*@
1090:    VecMTDot - Computes indefinite vector multiple dot products.
1091:    That is, it does NOT use the complex conjugate.

1093:    Collective on Vec

1095:    Input Parameters:
1096: +  x - one vector
1097: .  nv - number of vectors
1098: -  y - array of vectors.  Note that vectors are pointers

1100:    Output Parameter:
1101: .  val - array of the dot products

1103:    Notes for Users of Complex Numbers:
1104:    For complex vectors, VecMTDot() computes the indefinite form
1105: $      val = (x,y) = y^T x,
1106:    where y^T denotes the transpose of y.

1108:    Use VecMDot() for the inner product
1109: $      val = (x,y) = y^H x,
1110:    where y^H denotes the conjugate transpose of y.

1112:    Level: intermediate

1114:    Concepts: inner product^multiple
1115:    Concepts: vector^multiple inner products

1117: .seealso: VecMDot(), VecTDot()
1118: @*/
1119: PetscErrorCode  VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1120: {

1131:   VecCheckSameSize(x,1,*y,3);

1133:   PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1134:   (*x->ops->mtdot)(x,nv,y,val);
1135:   PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1136:   return(0);
1137: }

1139: /*@
1140:    VecMDot - Computes vector multiple dot products.

1142:    Collective on Vec

1144:    Input Parameters:
1145: +  x - one vector
1146: .  nv - number of vectors
1147: -  y - array of vectors.

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

1152:    Notes for Users of Complex Numbers:
1153:    For complex vectors, VecMDot() computes
1154: $     val = (x,y) = y^H x,
1155:    where y^H denotes the conjugate transpose of y.

1157:    Use VecMTDot() for the indefinite form
1158: $     val = (x,y) = y^T x,
1159:    where y^T denotes the transpose of y.

1161:    Level: intermediate

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

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

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

1184:   PetscLogEventBarrierBegin(VEC_MDotBarrier,x,*y,0,0,PetscObjectComm((PetscObject)x));
1185:   (*x->ops->mdot)(x,nv,y,val);
1186:   PetscLogEventBarrierEnd(VEC_MDotBarrier,x,*y,0,0,PetscObjectComm((PetscObject)x));
1187:   return(0);
1188: }

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

1193:    Logically Collective on Vec

1195:    Input Parameters:
1196: +  nv - number of scalars and x-vectors
1197: .  alpha - array of scalars
1198: .  y - one vector
1199: -  x - array of vectors

1201:    Level: intermediate

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

1205:    Concepts: BLAS

1207: .seealso: VecAXPY(), VecWAXPY(), VecAYPX()
1208: @*/
1209: PetscErrorCode  VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1210: {
1212:   PetscInt       i;

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

1227:   PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1228:   (*y->ops->maxpy)(y,nv,alpha,x);
1229:   PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1230:   PetscObjectStateIncrease((PetscObject)y);
1231:   return(0);
1232: }

1234: /*@
1235:    VecGetSubVector - Gets a vector representing part of another vector

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

1239:    Input Arguments:
1240: + X - vector from which to extract a subvector
1241: - is - index set representing portion of X to extract

1243:    Output Arguments:
1244: . Y - subvector corresponding to is

1246:    Level: advanced

1248:    Notes:
1249:    The subvector Y should be returned with VecRestoreSubVector().

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

1254: .seealso: MatCreateSubMatrix()
1255: @*/
1256: PetscErrorCode  VecGetSubVector(Vec X,IS is,Vec *Y)
1257: {
1258:   PetscErrorCode   ierr;
1259:   Vec              Z;

1265:   if (X->ops->getsubvector) {
1266:     (*X->ops->getsubvector)(X,is,&Z);
1267:   } else {                      /* Default implementation currently does no caching */
1268:     PetscInt  gstart,gend,start;
1269:     PetscBool contiguous,gcontiguous;
1270:     VecGetOwnershipRange(X,&gstart,&gend);
1271:     ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1272:     MPIU_Allreduce(&contiguous,&gcontiguous,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1273:     if (gcontiguous) {          /* We can do a no-copy implementation */
1274:       PetscInt    n,N,bs;
1275:       PetscMPIInt size;
1276:       PetscInt    state;

1278:       ISGetLocalSize(is,&n);
1279:       VecGetBlockSize(X,&bs);
1280:       if (n%bs || bs == 1) bs = -1; /* Do not decide block size if we do not have to */
1281:       MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1282:       VecLockGet(X,&state);
1283:       if (state) {
1284:         const PetscScalar *x;
1285:         VecGetArrayRead(X,&x);
1286:         if (size == 1) {
1287:           VecCreateSeqWithArray(PetscObjectComm((PetscObject)X),bs,n,(PetscScalar*)x+start,&Z);
1288:         } else {
1289:           ISGetSize(is,&N);
1290:           VecCreateMPIWithArray(PetscObjectComm((PetscObject)X),bs,n,N,(PetscScalar*)x+start,&Z);
1291:         }
1292:         VecRestoreArrayRead(X,&x);
1293:         VecLockPush(Z);
1294:       } else {
1295:         PetscScalar *x;
1296:         VecGetArray(X,&x);
1297:         if (size == 1) {
1298:           VecCreateSeqWithArray(PetscObjectComm((PetscObject)X),bs,n,x+start,&Z);
1299:         } else {
1300:           ISGetSize(is,&N);
1301:           VecCreateMPIWithArray(PetscObjectComm((PetscObject)X),bs,n,N,x+start,&Z);
1302:         }
1303:         VecRestoreArray(X,&x);
1304:       }
1305:     } else {                    /* Have to create a scatter and do a copy */
1306:       VecScatter scatter;
1307:       PetscInt   n,N;
1308:       ISGetLocalSize(is,&n);
1309:       ISGetSize(is,&N);
1310:       VecCreate(PetscObjectComm((PetscObject)is),&Z);
1311:       VecSetSizes(Z,n,N);
1312:       VecSetType(Z,((PetscObject)X)->type_name);
1313:       VecScatterCreate(X,is,Z,NULL,&scatter);
1314:       VecScatterBegin(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1315:       VecScatterEnd(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1316:       PetscObjectCompose((PetscObject)Z,"VecGetSubVector_Scatter",(PetscObject)scatter);
1317:       VecScatterDestroy(&scatter);
1318:     }
1319:   }
1320:   /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1321:   if (VecGetSubVectorSavedStateId < 0) {PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);}
1322:   PetscObjectComposedDataSetInt((PetscObject)Z,VecGetSubVectorSavedStateId,1);
1323:   *Y   = Z;
1324:   return(0);
1325: }

1327: /*@
1328:    VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()

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

1332:    Input Arguments:
1333: + X - vector from which subvector was obtained
1334: . is - index set representing the subset of X
1335: - Y - subvector being restored

1337:    Level: advanced

1339: .seealso: VecGetSubVector()
1340: @*/
1341: PetscErrorCode  VecRestoreSubVector(Vec X,IS is,Vec *Y)
1342: {

1350:   if (X->ops->restoresubvector) {
1351:     (*X->ops->restoresubvector)(X,is,Y);
1352:   } else {
1353:     PETSC_UNUSED PetscObjectState dummystate = 0;
1354:     PetscBool valid;
1355:     PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,dummystate,valid);
1356:     if (!valid) {
1357:       VecScatter scatter;

1359:       PetscObjectQuery((PetscObject)*Y,"VecGetSubVector_Scatter",(PetscObject*)&scatter);
1360:       if (scatter) {
1361:         VecScatterBegin(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1362:         VecScatterEnd(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1363:       }
1364:     }
1365:     VecDestroy(Y);
1366:   }
1367:   return(0);
1368: }

1370: /*@
1371:    VecGetLocalVectorRead - Maps the local portion of a vector into a
1372:    vector.  You must call VecRestoreLocalVectorRead() when the local
1373:    vector is no longer needed.

1375:    Not collective.

1377:    Input parameter:
1378: .  v - The vector for which the local vector is desired.

1380:    Output parameter:
1381: .  w - Upon exit this contains the local vector.

1383:    Level: beginner
1384:    
1385:    Notes:
1386:    This function is similar to VecGetArrayRead() which maps the local
1387:    portion into a raw pointer.  VecGetLocalVectorRead() is usually
1388:    almost as efficient as VecGetArrayRead() but in certain circumstances
1389:    VecGetLocalVectorRead() can be much more efficient than
1390:    VecGetArrayRead().  This is because the construction of a contiguous
1391:    array representing the vector data required by VecGetArrayRead() can
1392:    be an expensive operation for certain vector types.  For example, for
1393:    GPU vectors VecGetArrayRead() requires that the data between device
1394:    and host is synchronized.  

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

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

1409:   VecCheckSameLocalSize(v,1,w,2);
1410:   if (v->ops->getlocalvectorread) {
1411:     (*v->ops->getlocalvectorread)(v,w);
1412:   } else {
1413:     VecGetArrayRead(v,(const PetscScalar**)&a);
1414:     VecPlaceArray(w,a);
1415:   }
1416:   return(0);
1417: }

1419: /*@
1420:    VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1421:    previously mapped into a vector using VecGetLocalVectorRead().

1423:    Not collective.

1425:    Input parameter:
1426: .  v - The local portion of this vector was previously mapped into w using VecGetLocalVectorRead().
1427: .  w - The vector into which the local portion of v was mapped.

1429:    Level: beginner

1431: .seealso: VecGetLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1432: @*/
1433: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec w)
1434: {
1436:   PetscScalar    *a;

1441:   if (v->ops->restorelocalvectorread) {
1442:     (*v->ops->restorelocalvectorread)(v,w);
1443:   } else {
1444:     VecGetArrayRead(w,(const PetscScalar**)&a);
1445:     VecRestoreArrayRead(v,(const PetscScalar**)&a);
1446:     VecResetArray(w);
1447:   }
1448:   return(0);
1449: }

1451: /*@
1452:    VecGetLocalVector - Maps the local portion of a vector into a
1453:    vector.

1455:    Collective on v, not collective on w.

1457:    Input parameter:
1458: .  v - The vector for which the local vector is desired.

1460:    Output parameter:
1461: .  w - Upon exit this contains the local vector.

1463:    Level: beginner

1465:    Notes:
1466:    This function is similar to VecGetArray() which maps the local
1467:    portion into a raw pointer.  VecGetLocalVector() is usually about as
1468:    efficient as VecGetArray() but in certain circumstances
1469:    VecGetLocalVector() can be much more efficient than VecGetArray().
1470:    This is because the construction of a contiguous array representing
1471:    the vector data required by VecGetArray() can be an expensive
1472:    operation for certain vector types.  For example, for GPU vectors
1473:    VecGetArray() requires that the data between device and host is
1474:    synchronized.

1476: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1477: @*/
1478: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1479: {
1481:   PetscScalar    *a;

1486:   VecCheckSameLocalSize(v,1,w,2);
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: /*@
1497:    VecRestoreLocalVector - Unmaps the local portion of a vector
1498:    previously mapped into a vector using VecGetLocalVector().

1500:    Logically collective.

1502:    Input parameter:
1503: .  v - The local portion of this vector was previously mapped into w using VecGetLocalVector().
1504: .  w - The vector into which the local portion of v was mapped.

1506:    Level: beginner

1508: .seealso: VecGetLocalVector(), VecGetLocalVectorRead(), VecRestoreLocalVectorRead(), LocalVectorRead(), VecGetArrayRead(), VecGetArray()
1509: @*/
1510: PetscErrorCode VecRestoreLocalVector(Vec v,Vec w)
1511: {
1513:   PetscScalar    *a;

1518:   if (v->ops->restorelocalvector) {
1519:     (*v->ops->restorelocalvector)(v,w);
1520:   } else {
1521:     VecGetArray(w,&a);
1522:     VecRestoreArray(v,&a);
1523:     VecResetArray(w);
1524:   }
1525:   return(0);
1526: }

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

1537:    Logically Collective on Vec

1539:    Input Parameter:
1540: .  x - the vector

1542:    Output Parameter:
1543: .  a - location to put pointer to the array

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

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

1563:    Level: beginner

1565:    Concepts: vector^accessing local values

1567: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1568:           VecGetArrayPair(), VecRestoreArrayPair()
1569: @*/
1570: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1571: {

1576:   VecLocked(x,1);
1577:   if (x->petscnative) {
1578: #if defined(PETSC_HAVE_CUSP)
1579:     if (x->valid_GPU_array == PETSC_CUSP_GPU) {
1580:       VecCUSPCopyFromGPU(x);
1581:     } else if (x->valid_GPU_array == PETSC_CUSP_UNALLOCATED) {
1582:       VecCUSPAllocateCheckHost(x);
1583:     }
1584: #elif defined(PETSC_HAVE_VIENNACL)
1585:     if (x->valid_GPU_array == PETSC_VIENNACL_GPU) {
1586:       VecViennaCLCopyFromGPU(x);
1587:     } else if (x->valid_GPU_array == PETSC_VIENNACL_UNALLOCATED) {
1588:       VecViennaCLAllocateCheckHost(x);
1589:     }
1590: #elif defined(PETSC_HAVE_VECCUDA)
1591:     if (x->valid_GPU_array == PETSC_CUDA_GPU) {
1592:       VecCUDACopyFromGPU(x);
1593:     } else if (x->valid_GPU_array == PETSC_CUDA_UNALLOCATED) {
1594:       VecCUDAAllocateCheckHost(x);
1595:     }
1596: #endif
1597:     *a = *((PetscScalar**)x->data);
1598:   } else {
1599:     (*x->ops->getarray)(x,a);
1600:   }
1601:   return(0);
1602: }

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

1607:    Not Collective

1609:    Input Parameters:
1610: .  x - the vector

1612:    Output Parameter:
1613: .  a - the array

1615:    Level: beginner

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

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

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

1626: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1627: @*/
1628: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1629: {

1634:   if (x->petscnative) {
1635: #if defined(PETSC_HAVE_CUSP)
1636:     if (x->valid_GPU_array == PETSC_CUSP_GPU) {
1637:       VecCUSPCopyFromGPU(x);
1638:     }
1639: #elif defined(PETSC_HAVE_VIENNACL)
1640:     if (x->valid_GPU_array == PETSC_VIENNACL_GPU) {
1641:       VecViennaCLCopyFromGPU(x);
1642:     }
1643: #elif defined(PETSC_HAVE_VECCUDA)
1644:     if (x->valid_GPU_array == PETSC_CUDA_GPU) {
1645:       VecCUDACopyFromGPU(x);
1646:     }
1647: #endif
1648:     *a = *((PetscScalar **)x->data);
1649:   } else if (x->ops->getarrayread) {
1650:     (*x->ops->getarrayread)(x,a);
1651:   } else {
1652:     (*x->ops->getarray)(x,(PetscScalar**)a);
1653:   }
1654:   return(0);
1655: }

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

1662:    Logically Collective on Vec

1664:    Input Parameter:
1665: +  x - the vectors
1666: -  n - the number of vectors

1668:    Output Parameter:
1669: .  a - location to put pointer to the array

1671:    Fortran Note:
1672:    This routine is not supported in Fortran.

1674:    Level: intermediate

1676: .seealso: VecGetArray(), VecRestoreArrays()
1677: @*/
1678: PetscErrorCode  VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1679: {
1681:   PetscInt       i;
1682:   PetscScalar    **q;

1688:   if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1689:   PetscMalloc1(n,&q);
1690:   for (i=0; i<n; ++i) {
1691:     VecGetArray(x[i],&q[i]);
1692:   }
1693:   *a = q;
1694:   return(0);
1695: }

1697: /*@C
1698:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1699:    has been called.

1701:    Logically Collective on Vec

1703:    Input Parameters:
1704: +  x - the vector
1705: .  n - the number of vectors
1706: -  a - location of pointer to arrays obtained from VecGetArrays()

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

1714:    Fortran Note:
1715:    This routine is not supported in Fortran.

1717:    Level: intermediate

1719: .seealso: VecGetArrays(), VecRestoreArray()
1720: @*/
1721: PetscErrorCode  VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1722: {
1724:   PetscInt       i;
1725:   PetscScalar    **q = *a;


1732:   for (i=0; i<n; ++i) {
1733:     VecRestoreArray(x[i],&q[i]);
1734:   }
1735:   PetscFree(q);
1736:   return(0);
1737: }

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

1742:    Logically Collective on Vec

1744:    Input Parameters:
1745: +  x - the vector
1746: -  a - location of pointer to array obtained from VecGetArray()

1748:    Level: beginner

1750:    Notes:
1751:    For regular PETSc vectors this routine does not involve any copies. For
1752:    any special vectors that do not store local vector data in a contiguous
1753:    array, this routine will copy the data back into the underlying
1754:    vector data structure from the array obtained with VecGetArray().

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

1760:    Fortran Note:
1761:    This routine is used differently from Fortran 77
1762: $    Vec         x
1763: $    PetscScalar x_array(1)
1764: $    PetscOffset i_x
1765: $    PetscErrorCode ierr
1766: $       call VecGetArray(x,x_array,i_x,ierr)
1767: $
1768: $   Access first local entry in vector with
1769: $      value = x_array(i_x + 1)
1770: $
1771: $      ...... other code
1772: $       call VecRestoreArray(x,x_array,i_x,ierr)

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

1778: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1779:           VecGetArrayPair(), VecRestoreArrayPair()
1780: @*/
1781: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1782: {

1787:   if (x->petscnative) {
1788: #if defined(PETSC_HAVE_CUSP)
1789:     x->valid_GPU_array = PETSC_CUSP_CPU;
1790: #elif defined(PETSC_HAVE_VIENNACL)
1791:     x->valid_GPU_array = PETSC_VIENNACL_CPU;
1792: #elif defined(PETSC_HAVE_VECCUDA)
1793:     x->valid_GPU_array = PETSC_CUDA_CPU;
1794: #endif
1795:   } else {
1796:     (*x->ops->restorearray)(x,a);
1797:   }
1798:   if (a) *a = NULL;
1799:   PetscObjectStateIncrease((PetscObject)x);
1800:   return(0);
1801: }

1803: /*@C
1804:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

1806:    Not Collective

1808:    Input Parameters:
1809: +  vec - the vector
1810: -  array - the array

1812:    Level: beginner

1814: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1815: @*/
1816: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1817: {

1822:   if (x->petscnative) {
1823: #if defined(PETSC_HAVE_VIENNACL)
1824:     x->valid_GPU_array = PETSC_VIENNACL_CPU;
1825: #endif
1826:   } else if (x->ops->restorearrayread) {
1827:     (*x->ops->restorearrayread)(x,a);
1828:   } else {
1829:     (*x->ops->restorearray)(x,(PetscScalar**)a);
1830:   }
1831:   if (a) *a = NULL;
1832:   return(0);
1833: }

1835: /*@
1836:    VecPlaceArray - Allows one to replace the array in a vector with an
1837:    array provided by the user. This is useful to avoid copying an array
1838:    into a vector.

1840:    Not Collective

1842:    Input Parameters:
1843: +  vec - the vector
1844: -  array - the array

1846:    Notes:
1847:    You can return to the original array with a call to VecResetArray()

1849:    Level: developer

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

1853: @*/
1854: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
1855: {

1862:   if (vec->ops->placearray) {
1863:     (*vec->ops->placearray)(vec,array);
1864:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
1865:   PetscObjectStateIncrease((PetscObject)vec);
1866:   return(0);
1867: }

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

1874:    Not Collective

1876:    Input Parameters:
1877: +  vec - the vector
1878: -  array - the array

1880:    Notes:
1881:    This permanently replaces the array and frees the memory associated
1882:    with the old array.

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

1887:    Not supported from Fortran

1889:    Level: developer

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

1893: @*/
1894: PetscErrorCode  VecReplaceArray(Vec vec,const PetscScalar array[])
1895: {

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

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

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

1915:     Collective on Vec

1917:     Input Parameters:
1918: +   x - a vector to mimic
1919: -   n - the number of vectors to obtain

1921:     Output Parameters:
1922: +   y - Fortran90 pointer to the array of vectors
1923: -   ierr - error code

1925:     Example of Usage:
1926: .vb
1927: #include <petsc/finclude/petscvec.h90>

1929:     Vec x
1930:     Vec, pointer :: y(:)
1931:     ....
1932:     call VecDuplicateVecsF90(x,2,y,ierr)
1933:     call VecSet(y(2),alpha,ierr)
1934:     call VecSet(y(2),alpha,ierr)
1935:     ....
1936:     call VecDestroyVecsF90(2,y,ierr)
1937: .ve

1939:     Notes:
1940:     Not yet supported for all F90 compilers

1942:     Use VecDestroyVecsF90() to free the space.

1944:     Level: beginner

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

1948: M*/

1950: /*MC
1951:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
1952:     VecGetArrayF90().

1954:     Synopsis:
1955:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

1957:     Logically Collective on Vec

1959:     Input Parameters:
1960: +   x - vector
1961: -   xx_v - the Fortran90 pointer to the array

1963:     Output Parameter:
1964: .   ierr - error code

1966:     Example of Usage:
1967: .vb
1968: #include <petsc/finclude/petscvec.h90>

1970:     PetscScalar, pointer :: xx_v(:)
1971:     ....
1972:     call VecGetArrayF90(x,xx_v,ierr)
1973:     xx_v(3) = a
1974:     call VecRestoreArrayF90(x,xx_v,ierr)
1975: .ve

1977:     Level: beginner

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

1981: M*/

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

1986:     Synopsis:
1987:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

1989:     Collective on Vec

1991:     Input Parameters:
1992: +   n - the number of vectors previously obtained
1993: -   x - pointer to array of vector pointers

1995:     Output Parameter:
1996: .   ierr - error code

1998:     Notes:
1999:     Not yet supported for all F90 compilers

2001:     Level: beginner

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

2005: M*/

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

2013:     Synopsis:
2014:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2016:     Logically Collective on Vec

2018:     Input Parameter:
2019: .   x - vector

2021:     Output Parameters:
2022: +   xx_v - the Fortran90 pointer to the array
2023: -   ierr - error code

2025:     Example of Usage:
2026: .vb
2027: #include <petsc/finclude/petscvec.h90>

2029:     PetscScalar, pointer :: xx_v(:)
2030:     ....
2031:     call VecGetArrayF90(x,xx_v,ierr)
2032:     xx_v(3) = a
2033:     call VecRestoreArrayF90(x,xx_v,ierr)
2034: .ve

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

2038:     Level: beginner

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

2042: M*/

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

2050:     Synopsis:
2051:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2053:     Logically Collective on Vec

2055:     Input Parameter:
2056: .   x - vector

2058:     Output Parameters:
2059: +   xx_v - the Fortran90 pointer to the array
2060: -   ierr - error code

2062:     Example of Usage:
2063: .vb
2064: #include <petsc/finclude/petscvec.h90>

2066:     PetscScalar, pointer :: xx_v(:)
2067:     ....
2068:     call VecGetArrayReadF90(x,xx_v,ierr)
2069:     a = xx_v(3)
2070:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2071: .ve

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

2075:     Level: beginner

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

2079: M*/

2081: /*MC
2082:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2083:     VecGetArrayReadF90().

2085:     Synopsis:
2086:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2088:     Logically Collective on Vec

2090:     Input Parameters:
2091: +   x - vector
2092: -   xx_v - the Fortran90 pointer to the array

2094:     Output Parameter:
2095: .   ierr - error code

2097:     Example of Usage:
2098: .vb
2099: #include <petsc/finclude/petscvec.h90>

2101:     PetscScalar, pointer :: xx_v(:)
2102:     ....
2103:     call VecGetArrayReadF90(x,xx_v,ierr)
2104:     a = xx_v(3)
2105:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2106: .ve

2108:     Level: beginner

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

2112: M*/

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

2119:    Logically Collective

2121:    Input Parameter:
2122: +  x - the vector
2123: .  m - first dimension of two dimensional array
2124: .  n - second dimension of two dimensional array
2125: .  mstart - first index you will use in first coordinate direction (often 0)
2126: -  nstart - first index in the second coordinate direction (often 0)

2128:    Output Parameter:
2129: .  a - location to put pointer to the array

2131:    Level: developer

2133:   Notes:
2134:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2135:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2136:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2137:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

2143: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2144:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2145:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2146: @*/
2147: PetscErrorCode  VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2148: {
2150:   PetscInt       i,N;
2151:   PetscScalar    *aa;

2157:   VecGetLocalSize(x,&N);
2158:   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);
2159:   VecGetArray(x,&aa);

2161:   PetscMalloc1(m,a);
2162:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2163:   *a -= mstart;
2164:   return(0);
2165: }

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

2170:    Logically Collective

2172:    Input Parameters:
2173: +  x - the vector
2174: .  m - first dimension of two dimensional array
2175: .  n - second dimension of the two dimensional array
2176: .  mstart - first index you will use in first coordinate direction (often 0)
2177: .  nstart - first index in the second coordinate direction (often 0)
2178: -  a - location of pointer to array obtained from VecGetArray2d()

2180:    Level: developer

2182:    Notes:
2183:    For regular PETSc vectors this routine does not involve any copies. For
2184:    any special vectors that do not store local vector data in a contiguous
2185:    array, this routine will copy the data back into the underlying
2186:    vector data structure from the array obtained with VecGetArray().

2188:    This routine actually zeros out the a pointer.

2190: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2191:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2192:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2193: @*/
2194: PetscErrorCode  VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2195: {
2197:   void           *dummy;

2203:   dummy = (void*)(*a + mstart);
2204:   PetscFree(dummy);
2205:   VecRestoreArray(x,NULL);
2206:   return(0);
2207: }

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

2214:    Logically Collective

2216:    Input Parameter:
2217: +  x - the vector
2218: .  m - first dimension of two dimensional array
2219: -  mstart - first index you will use in first coordinate direction (often 0)

2221:    Output Parameter:
2222: .  a - location to put pointer to the array

2224:    Level: developer

2226:   Notes:
2227:    For a vector obtained from DMCreateLocalVector() mstart are likely
2228:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2229:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2233: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2234:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2235:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2236: @*/
2237: PetscErrorCode  VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2238: {
2240:   PetscInt       N;

2246:   VecGetLocalSize(x,&N);
2247:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2248:   VecGetArray(x,a);
2249:   *a  -= mstart;
2250:   return(0);
2251: }

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

2256:    Logically Collective

2258:    Input Parameters:
2259: +  x - the vector
2260: .  m - first dimension of two dimensional array
2261: .  mstart - first index you will use in first coordinate direction (often 0)
2262: -  a - location of pointer to array obtained from VecGetArray21()

2264:    Level: developer

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

2272:    This routine actually zeros out the a pointer.

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

2276: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2277:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2278:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2279: @*/
2280: PetscErrorCode  VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2281: {

2287:   VecRestoreArray(x,NULL);
2288:   return(0);
2289: }


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

2297:    Logically Collective

2299:    Input Parameter:
2300: +  x - the vector
2301: .  m - first dimension of three dimensional array
2302: .  n - second dimension of three dimensional array
2303: .  p - third dimension of three dimensional array
2304: .  mstart - first index you will use in first coordinate direction (often 0)
2305: .  nstart - first index in the second coordinate direction (often 0)
2306: -  pstart - first index in the third coordinate direction (often 0)

2308:    Output Parameter:
2309: .  a - location to put pointer to the array

2311:    Level: developer

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

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

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

2323: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2324:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2325:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2326: @*/
2327: PetscErrorCode  VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2328: {
2330:   PetscInt       i,N,j;
2331:   PetscScalar    *aa,**b;

2337:   VecGetLocalSize(x,&N);
2338:   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);
2339:   VecGetArray(x,&aa);

2341:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
2342:   b    = (PetscScalar**)((*a) + m);
2343:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2344:   for (i=0; i<m; i++)
2345:     for (j=0; j<n; j++)
2346:       b[i*n+j] = aa + i*n*p + j*p - pstart;

2348:   *a -= mstart;
2349:   return(0);
2350: }

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

2355:    Logically Collective

2357:    Input Parameters:
2358: +  x - the vector
2359: .  m - first dimension of three dimensional array
2360: .  n - second dimension of the three dimensional array
2361: .  p - third dimension of the three dimensional array
2362: .  mstart - first index you will use in first coordinate direction (often 0)
2363: .  nstart - first index in the second coordinate direction (often 0)
2364: .  pstart - first index in the third coordinate direction (often 0)
2365: -  a - location of pointer to array obtained from VecGetArray3d()

2367:    Level: developer

2369:    Notes:
2370:    For regular PETSc vectors this routine does not involve any copies. For
2371:    any special vectors that do not store local vector data in a contiguous
2372:    array, this routine will copy the data back into the underlying
2373:    vector data structure from the array obtained with VecGetArray().

2375:    This routine actually zeros out the a pointer.

2377: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2378:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2379:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2380: @*/
2381: PetscErrorCode  VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2382: {
2384:   void           *dummy;

2390:   dummy = (void*)(*a + mstart);
2391:   PetscFree(dummy);
2392:   VecRestoreArray(x,NULL);
2393:   return(0);
2394: }

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

2401:    Logically Collective

2403:    Input Parameter:
2404: +  x - the vector
2405: .  m - first dimension of four dimensional array
2406: .  n - second dimension of four dimensional array
2407: .  p - third dimension of four dimensional array
2408: .  q - fourth dimension of 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)

2414:    Output Parameter:
2415: .  a - location to put pointer to the array

2417:    Level: beginner

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

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

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

2429: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2430:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2431:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2432: @*/
2433: PetscErrorCode  VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2434: {
2436:   PetscInt       i,N,j,k;
2437:   PetscScalar    *aa,***b,**c;

2443:   VecGetLocalSize(x,&N);
2444:   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);
2445:   VecGetArray(x,&aa);

2447:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
2448:   b    = (PetscScalar***)((*a) + m);
2449:   c    = (PetscScalar**)(b + m*n);
2450:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2451:   for (i=0; i<m; i++)
2452:     for (j=0; j<n; j++)
2453:       b[i*n+j] = c + i*n*p + j*p - pstart;
2454:   for (i=0; i<m; i++)
2455:     for (j=0; j<n; j++)
2456:       for (k=0; k<p; k++)
2457:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
2458:   *a -= mstart;
2459:   return(0);
2460: }

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

2465:    Logically Collective

2467:    Input Parameters:
2468: +  x - the vector
2469: .  m - first dimension of four dimensional array
2470: .  n - second dimension of the four dimensional array
2471: .  p - third dimension of the four dimensional array
2472: .  q - fourth dimension of the four dimensional array
2473: .  mstart - first index you will use in first coordinate direction (often 0)
2474: .  nstart - first index in the second coordinate direction (often 0)
2475: .  pstart - first index in the third coordinate direction (often 0)
2476: .  qstart - first index in the fourth coordinate direction (often 0)
2477: -  a - location of pointer to array obtained from VecGetArray4d()

2479:    Level: beginner

2481:    Notes:
2482:    For regular PETSc vectors this routine does not involve any copies. For
2483:    any special vectors that do not store local vector data in a contiguous
2484:    array, this routine will copy the data back into the underlying
2485:    vector data structure from the array obtained with VecGetArray().

2487:    This routine actually zeros out the a pointer.

2489: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2490:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2491:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2492: @*/
2493: PetscErrorCode  VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2494: {
2496:   void           *dummy;

2502:   dummy = (void*)(*a + mstart);
2503:   PetscFree(dummy);
2504:   VecRestoreArray(x,NULL);
2505:   return(0);
2506: }

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

2513:    Logically Collective

2515:    Input Parameter:
2516: +  x - the vector
2517: .  m - first dimension of two dimensional array
2518: .  n - second dimension of two dimensional array
2519: .  mstart - first index you will use in first coordinate direction (often 0)
2520: -  nstart - first index in the second coordinate direction (often 0)

2522:    Output Parameter:
2523: .  a - location to put pointer to the array

2525:    Level: developer

2527:   Notes:
2528:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2529:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2530:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2531:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

2537: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2538:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2539:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2540: @*/
2541: PetscErrorCode  VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2542: {
2543:   PetscErrorCode    ierr;
2544:   PetscInt          i,N;
2545:   const PetscScalar *aa;

2551:   VecGetLocalSize(x,&N);
2552:   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);
2553:   VecGetArrayRead(x,&aa);

2555:   PetscMalloc1(m,a);
2556:   for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
2557:   *a -= mstart;
2558:   return(0);
2559: }

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

2564:    Logically Collective

2566:    Input Parameters:
2567: +  x - the vector
2568: .  m - first dimension of two dimensional array
2569: .  n - second dimension of the two dimensional array
2570: .  mstart - first index you will use in first coordinate direction (often 0)
2571: .  nstart - first index in the second coordinate direction (often 0)
2572: -  a - location of pointer to array obtained from VecGetArray2d()

2574:    Level: developer

2576:    Notes:
2577:    For regular PETSc vectors this routine does not involve any copies. For
2578:    any special vectors that do not store local vector data in a contiguous
2579:    array, this routine will copy the data back into the underlying
2580:    vector data structure from the array obtained with VecGetArray().

2582:    This routine actually zeros out the a pointer.

2584: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2585:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2586:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2587: @*/
2588: PetscErrorCode  VecRestoreArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2589: {
2591:   void           *dummy;

2597:   dummy = (void*)(*a + mstart);
2598:   PetscFree(dummy);
2599:   VecRestoreArrayRead(x,NULL);
2600:   return(0);
2601: }

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

2608:    Logically Collective

2610:    Input Parameter:
2611: +  x - the vector
2612: .  m - first dimension of two dimensional array
2613: -  mstart - first index you will use in first coordinate direction (often 0)

2615:    Output Parameter:
2616: .  a - location to put pointer to the array

2618:    Level: developer

2620:   Notes:
2621:    For a vector obtained from DMCreateLocalVector() mstart are likely
2622:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2623:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2627: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2628:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2629:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2630: @*/
2631: PetscErrorCode  VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2632: {
2634:   PetscInt       N;

2640:   VecGetLocalSize(x,&N);
2641:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2642:   VecGetArrayRead(x,(const PetscScalar**)a);
2643:   *a  -= mstart;
2644:   return(0);
2645: }

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

2650:    Logically Collective

2652:    Input Parameters:
2653: +  x - the vector
2654: .  m - first dimension of two dimensional array
2655: .  mstart - first index you will use in first coordinate direction (often 0)
2656: -  a - location of pointer to array obtained from VecGetArray21()

2658:    Level: developer

2660:    Notes:
2661:    For regular PETSc vectors this routine does not involve any copies. For
2662:    any special vectors that do not store local vector data in a contiguous
2663:    array, this routine will copy the data back into the underlying
2664:    vector data structure from the array obtained with VecGetArray1dRead().

2666:    This routine actually zeros out the a pointer.

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

2670: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2671:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2672:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2673: @*/
2674: PetscErrorCode  VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2675: {

2681:   VecRestoreArrayRead(x,NULL);
2682:   return(0);
2683: }


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

2691:    Logically Collective

2693:    Input Parameter:
2694: +  x - the vector
2695: .  m - first dimension of three dimensional array
2696: .  n - second dimension of three dimensional array
2697: .  p - third dimension of three dimensional array
2698: .  mstart - first index you will use in first coordinate direction (often 0)
2699: .  nstart - first index in the second coordinate direction (often 0)
2700: -  pstart - first index in the third coordinate direction (often 0)

2702:    Output Parameter:
2703: .  a - location to put pointer to the array

2705:    Level: developer

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

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

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

2717: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2718:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2719:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2720: @*/
2721: PetscErrorCode  VecGetArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2722: {
2723:   PetscErrorCode    ierr;
2724:   PetscInt          i,N,j;
2725:   const PetscScalar *aa;
2726:   PetscScalar       **b;

2732:   VecGetLocalSize(x,&N);
2733:   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);
2734:   VecGetArrayRead(x,&aa);

2736:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
2737:   b    = (PetscScalar**)((*a) + m);
2738:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2739:   for (i=0; i<m; i++)
2740:     for (j=0; j<n; j++)
2741:       b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;

2743:   *a -= mstart;
2744:   return(0);
2745: }

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

2750:    Logically Collective

2752:    Input Parameters:
2753: +  x - the vector
2754: .  m - first dimension of three dimensional array
2755: .  n - second dimension of the three dimensional array
2756: .  p - third dimension of the three dimensional array
2757: .  mstart - first index you will use in first coordinate direction (often 0)
2758: .  nstart - first index in the second coordinate direction (often 0)
2759: .  pstart - first index in the third coordinate direction (often 0)
2760: -  a - location of pointer to array obtained from VecGetArray3dRead()

2762:    Level: developer

2764:    Notes:
2765:    For regular PETSc vectors this routine does not involve any copies. For
2766:    any special vectors that do not store local vector data in a contiguous
2767:    array, this routine will copy the data back into the underlying
2768:    vector data structure from the array obtained with VecGetArray().

2770:    This routine actually zeros out the a pointer.

2772: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2773:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2774:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2775: @*/
2776: PetscErrorCode  VecRestoreArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2777: {
2779:   void           *dummy;

2785:   dummy = (void*)(*a + mstart);
2786:   PetscFree(dummy);
2787:   VecRestoreArrayRead(x,NULL);
2788:   return(0);
2789: }

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

2796:    Logically Collective

2798:    Input Parameter:
2799: +  x - the vector
2800: .  m - first dimension of four dimensional array
2801: .  n - second dimension of four dimensional array
2802: .  p - third dimension of four dimensional array
2803: .  q - fourth dimension of four dimensional array
2804: .  mstart - first index you will use in first coordinate direction (often 0)
2805: .  nstart - first index in the second coordinate direction (often 0)
2806: .  pstart - first index in the third coordinate direction (often 0)
2807: -  qstart - first index in the fourth coordinate direction (often 0)

2809:    Output Parameter:
2810: .  a - location to put pointer to the array

2812:    Level: beginner

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

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

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

2824: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2825:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2826:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2827: @*/
2828: PetscErrorCode  VecGetArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2829: {
2830:   PetscErrorCode    ierr;
2831:   PetscInt          i,N,j,k;
2832:   const PetscScalar *aa;
2833:   PetscScalar       ***b,**c;

2839:   VecGetLocalSize(x,&N);
2840:   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);
2841:   VecGetArrayRead(x,&aa);

2843:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
2844:   b    = (PetscScalar***)((*a) + m);
2845:   c    = (PetscScalar**)(b + m*n);
2846:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2847:   for (i=0; i<m; i++)
2848:     for (j=0; j<n; j++)
2849:       b[i*n+j] = c + i*n*p + j*p - pstart;
2850:   for (i=0; i<m; i++)
2851:     for (j=0; j<n; j++)
2852:       for (k=0; k<p; k++)
2853:         c[i*n*p+j*p+k] = (PetscScalar*) aa + i*n*p*q + j*p*q + k*q - qstart;
2854:   *a -= mstart;
2855:   return(0);
2856: }

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

2861:    Logically Collective

2863:    Input Parameters:
2864: +  x - the vector
2865: .  m - first dimension of four dimensional array
2866: .  n - second dimension of the four dimensional array
2867: .  p - third dimension of the four dimensional array
2868: .  q - fourth dimension of the four dimensional array
2869: .  mstart - first index you will use in first coordinate direction (often 0)
2870: .  nstart - first index in the second coordinate direction (often 0)
2871: .  pstart - first index in the third coordinate direction (often 0)
2872: .  qstart - first index in the fourth coordinate direction (often 0)
2873: -  a - location of pointer to array obtained from VecGetArray4dRead()

2875:    Level: beginner

2877:    Notes:
2878:    For regular PETSc vectors this routine does not involve any copies. For
2879:    any special vectors that do not store local vector data in a contiguous
2880:    array, this routine will copy the data back into the underlying
2881:    vector data structure from the array obtained with VecGetArray().

2883:    This routine actually zeros out the a pointer.

2885: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2886:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2887:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2888: @*/
2889: PetscErrorCode  VecRestoreArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2890: {
2892:   void           *dummy;

2898:   dummy = (void*)(*a + mstart);
2899:   PetscFree(dummy);
2900:   VecRestoreArrayRead(x,NULL);
2901:   return(0);
2902: }

2904: #if defined(PETSC_USE_DEBUG)

2906: /*@
2907:    VecLockGet  - Gets the current lock status of a vector

2909:    Logically Collective on Vec

2911:    Input Parameter:
2912: .  x - the vector

2914:    Output Parameter:
2915: .  state - greater than zero indicates the vector is still locked

2917:    Level: beginner

2919:    Concepts: vector^accessing local values

2921: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockPush(), VecLockGet()
2922: @*/
2923: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
2924: {
2927:   *state = x->lock;
2928:   return(0);
2929: }

2931: /*@
2932:    VecLockPush  - Lock a vector from writing

2934:    Logically Collective on Vec

2936:    Input Parameter:
2937: .  x - the vector

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

2941:     Call VecLockPop() to remove the latest lock

2943:    Level: beginner

2945:    Concepts: vector^accessing local values

2947: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockPop(), VecLockGet()
2948: @*/
2949: PetscErrorCode VecLockPush(Vec x)
2950: {
2953:   x->lock++;
2954:   return(0);
2955: }

2957: /*@
2958:    VecLockPop  - Unlock a vector from writing

2960:    Logically Collective on Vec

2962:    Input Parameter:
2963: .  x - the vector

2965:    Level: beginner

2967:    Concepts: vector^accessing local values

2969: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockPush(), VecLockGet()
2970: @*/
2971: PetscErrorCode VecLockPop(Vec x)
2972: {
2975:   x->lock--;
2976:   if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector has been unlocked too many times");
2977:   return(0);
2978: }

2980: #endif