Actual source code: rvector.c

petsc-master 2016-12-03
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;

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

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

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

 50:    Logically Collective on Vec

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

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

 58:    Level: advanced

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

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

 76:   VecCheckSameSize(x,1,y,2);
 77:   (*x->ops->maxpointwisedivide)(x,y,max);
 78:   return(0);
 79: }

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

 86:    Collective on Vec

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

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

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

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

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

110:    Level: intermediate

112:    Concepts: inner product
113:    Concepts: vector^inner product

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

128:   VecCheckSameSize(x,1,y,1);

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

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

141:    Collective on Vec

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

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

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

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

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

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

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

164:    Level: intermediate

166:    Concepts: inner product
167:    Concepts: vector^inner product

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

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

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

187:    Collective on Vec

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

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

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

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

208:    Level: intermediate

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

215:    Concepts: norm
216:    Concepts: vector^norm

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

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

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

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

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

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

255:    Not Collective

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

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

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

272:    Level: intermediate

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

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

284:    Concepts: norm
285:    Concepts: vector^norm

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

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


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

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

312:    Collective on Vec

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

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

321:    Level: intermediate

323:    Concepts: vector^normalizing
324:    Concepts: normalizing^vector

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

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

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

353:    Collective on Vec

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

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

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

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

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

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

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

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

392:    Collective on Vec

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

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

401:    Level: intermediate

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

406:    This returns the smallest index with the minumum value

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

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

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

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

433:    Collective on Vec

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

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

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

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

450:    Level: intermediate

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

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

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

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

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

482:    Not collective on Vec

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

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

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

495:    Level: intermediate

497:    Concepts: vector^scaling
498:    Concepts: scaling^vector

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

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

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

536:    Logically Collective on Vec

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

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

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

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

555:    Level: beginner

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

559:    Concepts: vector^setting to constant

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

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

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

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


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

604:    Logically Collective on Vec

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

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

613:    Level: intermediate

615:    Notes: x and y MUST be different vectors

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

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

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

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

648: /*@
649:    VecAXPBY - Computes y = alpha x + beta y.

651:    Logically Collective on Vec

653:    Input Parameters:
654: +  alpha,beta - the scalars
655: -  x, y  - the vectors

657:    Output Parameter:
658: .  y - output vector

660:    Level: intermediate

662:    Notes: x and y MUST be different vectors

664:    Concepts: BLAS
665:    Concepts: vector^BLAS

667: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
668: @*/
669: PetscErrorCode  VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
670: {

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

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

693: /*@
694:    VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z

696:    Logically Collective on Vec

698:    Input Parameters:
699: +  alpha,beta, gamma - the scalars
700: -  x, y, z  - the vectors

702:    Output Parameter:
703: .  z - output vector

705:    Level: intermediate

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

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

711:    Concepts: BLAS
712:    Concepts: vector^BLAS

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

729:   VecCheckSameSize(x,1,y,5);
730:   VecCheckSameSize(x,1,z,6);
731:   if (x == y || x == z) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
732:   if (y == z) SETERRQ(PetscObjectComm((PetscObject)y),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");

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

746: /*@
747:    VecAYPX - Computes y = x + alpha y.

749:    Logically Collective on Vec

751:    Input Parameters:
752: +  alpha - the scalar
753: -  x, y  - the vectors

755:    Output Parameter:
756: .  y - output vector

758:    Level: intermediate

760:    Notes: x and y MUST be different vectors

762:    Concepts: vector^BLAS
763:    Concepts: BLAS

765: .seealso: VecAXPY(), VecWAXPY()
766: @*/
767: PetscErrorCode  VecAYPX(Vec y,PetscScalar alpha,Vec x)
768: {

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

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


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

794:    Logically Collective on Vec

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

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

803:    Level: intermediate

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

807:    Concepts: vector^BLAS
808:    Concepts: BLAS

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

825:   VecCheckSameSize(x,3,y,4);
826:   VecCheckSameSize(x,3,w,1);
827:   if (w == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector y, suggest VecAXPY()");
828:   if (w == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector x, suggest VecAYPX()");

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


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

844:    Not Collective

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

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

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

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

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

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

873:    Level: beginner

875:    Concepts: vector^setting values

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

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

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

903:     Not Collective

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

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

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

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

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

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

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

926:    Level: beginner

928:    Concepts: vector^getting values

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

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

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

951:    Not Collective

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

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

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

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

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

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

980:    Level: intermediate

982:    Concepts: vector^setting values blocked

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

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


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

1010:    Not Collective

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

1021:    Level: intermediate

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

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

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

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

1035:    Concepts: vector^setting values with local numbering

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

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

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

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

1077:    Not Collective

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

1088:    Level: intermediate

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

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

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

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


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

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

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

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

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

1140:    Collective on Vec

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

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

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

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

1159:    Level: intermediate

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

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

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

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

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

1191:    Collective on Vec

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

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

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

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

1210:    Level: intermediate

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

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

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

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

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

1244:    Logically Collective on Vec

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

1252:    Level: intermediate

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

1256:    Concepts: BLAS

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

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

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

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

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

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

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

1299:    Level: advanced

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

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

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

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

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

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

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

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

1392:    Level: advanced

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

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

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

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

1430:    Not collective.

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

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

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

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

1454: .seealso: VecRestoreLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1455: @*/
1458: PetscErrorCode VecGetLocalVectorRead(Vec v,Vec w)
1459: {
1461:   PetscScalar    *a;

1466:   VecCheckSameLocalSize(v,1,w,2);
1467:   if (v->ops->getlocalvectorread) {
1468:     (*v->ops->getlocalvectorread)(v,w);
1469:   } else {
1470:     VecGetArrayRead(v,(const PetscScalar**)&a);
1471:     VecPlaceArray(w,a);
1472:   }
1473:   return(0);
1474: }

1476: /*@
1477:    VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1478:    previously mapped into a vector using VecGetLocalVectorRead().

1480:    Not collective.

1482:    Input parameter:
1483: .  v - The local portion of this vector was previously mapped into w using VecGetLocalVectorRead().
1484: .  w - The vector into which the local portion of v was mapped.

1486:    Level: beginner

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

1500:   if (v->ops->restorelocalvectorread) {
1501:     (*v->ops->restorelocalvectorread)(v,w);
1502:   } else {
1503:     VecGetArrayRead(w,(const PetscScalar**)&a);
1504:     VecRestoreArrayRead(v,(const PetscScalar**)&a);
1505:     VecResetArray(w);
1506:   }
1507:   return(0);
1508: }

1510: /*@
1511:    VecGetLocalVector - Maps the local portion of a vector into a
1512:    vector.

1514:    Collective on v, not collective on w.

1516:    Input parameter:
1517: .  v - The vector for which the local vector is desired.

1519:    Output parameter:
1520: .  w - Upon exit this contains the local vector.

1522:    Level: beginner

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

1535: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1536: @*/
1539: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1540: {
1542:   PetscScalar    *a;

1547:   VecCheckSameLocalSize(v,1,w,2);
1548:   if (v->ops->getlocalvector) {
1549:     (*v->ops->getlocalvector)(v,w);
1550:   } else {
1551:     VecGetArray(v,&a);
1552:     VecPlaceArray(w,a);
1553:   }
1554:   return(0);
1555: }

1557: /*@
1558:    VecRestoreLocalVector - Unmaps the local portion of a vector
1559:    previously mapped into a vector using VecGetLocalVector().

1561:    Logically collective.

1563:    Input parameter:
1564: .  v - The local portion of this vector was previously mapped into w using VecGetLocalVector().
1565: .  w - The vector into which the local portion of v was mapped.

1567:    Level: beginner

1569: .seealso: VecGetLocalVector(), VecGetLocalVectorRead(), VecRestoreLocalVectorRead(), LocalVectorRead(), VecGetArrayRead(), VecGetArray()
1570: @*/
1573: PetscErrorCode VecRestoreLocalVector(Vec v,Vec w)
1574: {
1576:   PetscScalar    *a;

1581:   if (v->ops->restorelocalvector) {
1582:     (*v->ops->restorelocalvector)(v,w);
1583:   } else {
1584:     VecGetArray(w,&a);
1585:     VecRestoreArray(v,&a);
1586:     VecResetArray(w);
1587:   }
1588:   return(0);
1589: }

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

1602:    Logically Collective on Vec

1604:    Input Parameter:
1605: .  x - the vector

1607:    Output Parameter:
1608: .  a - location to put pointer to the array

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

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

1628:    Level: beginner

1630:    Concepts: vector^accessing local values

1632: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d()
1633: @*/
1634: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1635: {

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

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

1673:    Not Collective

1675:    Input Parameters:
1676: .  x - the vector

1678:    Output Parameter:
1679: .  a - the array

1681:    Level: beginner

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

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

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

1692: .seealso: VecGetArray(), VecRestoreArray()
1693: @*/
1694: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1695: {

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

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

1730:    Logically Collective on Vec

1732:    Input Parameter:
1733: +  x - the vectors
1734: -  n - the number of vectors

1736:    Output Parameter:
1737: .  a - location to put pointer to the array

1739:    Fortran Note:
1740:    This routine is not supported in Fortran.

1742:    Level: intermediate

1744: .seealso: VecGetArray(), VecRestoreArrays()
1745: @*/
1746: PetscErrorCode  VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1747: {
1749:   PetscInt       i;
1750:   PetscScalar    **q;

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

1767: /*@C
1768:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1769:    has been called.

1771:    Logically Collective on Vec

1773:    Input Parameters:
1774: +  x - the vector
1775: .  n - the number of vectors
1776: -  a - location of pointer to arrays obtained from VecGetArrays()

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

1784:    Fortran Note:
1785:    This routine is not supported in Fortran.

1787:    Level: intermediate

1789: .seealso: VecGetArrays(), VecRestoreArray()
1790: @*/
1791: PetscErrorCode  VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1792: {
1794:   PetscInt       i;
1795:   PetscScalar    **q = *a;


1802:   for (i=0; i<n; ++i) {
1803:     VecRestoreArray(x[i],&q[i]);
1804:   }
1805:   PetscFree(q);
1806:   return(0);
1807: }

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

1814:    Logically Collective on Vec

1816:    Input Parameters:
1817: +  x - the vector
1818: -  a - location of pointer to array obtained from VecGetArray()

1820:    Level: beginner

1822:    Notes:
1823:    For regular PETSc vectors this routine does not involve any copies. For
1824:    any special vectors that do not store local vector data in a contiguous
1825:    array, this routine will copy the data back into the underlying
1826:    vector data structure from the array obtained with VecGetArray().

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

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

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

1850: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d()
1851: @*/
1852: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1853: {

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

1876: /*@C
1877:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

1879:    Not Collective

1881:    Input Parameters:
1882: +  vec - the vector
1883: -  array - the array

1885:    Level: beginner

1887: .seealso: VecGetArray(), VecRestoreArray()
1888: @*/
1889: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1890: {

1895:   if (x->petscnative) {
1896: #if defined(PETSC_HAVE_VIENNACL)
1897:     x->valid_GPU_array = PETSC_VIENNACL_CPU;
1898: #endif
1899:   } else if (x->ops->restorearrayread) {
1900:     (*x->ops->restorearrayread)(x,a);
1901:   } else {
1902:     (*x->ops->restorearray)(x,(PetscScalar**)a);
1903:   }
1904:   if (a) *a = NULL;
1905:   return(0);
1906: }

1910: /*@
1911:    VecPlaceArray - Allows one to replace the array in a vector with an
1912:    array provided by the user. This is useful to avoid copying an array
1913:    into a vector.

1915:    Not Collective

1917:    Input Parameters:
1918: +  vec - the vector
1919: -  array - the array

1921:    Notes:
1922:    You can return to the original array with a call to VecResetArray()

1924:    Level: developer

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

1928: @*/
1929: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
1930: {

1937:   if (vec->ops->placearray) {
1938:     (*vec->ops->placearray)(vec,array);
1939:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
1940:   PetscObjectStateIncrease((PetscObject)vec);
1941:   return(0);
1942: }

1946: /*@C
1947:    VecReplaceArray - Allows one to replace the array in a vector with an
1948:    array provided by the user. This is useful to avoid copying an array
1949:    into a vector.

1951:    Not Collective

1953:    Input Parameters:
1954: +  vec - the vector
1955: -  array - the array

1957:    Notes:
1958:    This permanently replaces the array and frees the memory associated
1959:    with the old array.

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

1964:    Not supported from Fortran

1966:    Level: developer

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

1970: @*/
1971: PetscErrorCode  VecReplaceArray(Vec vec,const PetscScalar array[])
1972: {

1978:   if (vec->ops->replacearray) {
1979:     (*vec->ops->replacearray)(vec,array);
1980:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
1981:   PetscObjectStateIncrease((PetscObject)vec);
1982:   return(0);
1983: }

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

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

1992:     Collective on Vec

1994:     Input Parameters:
1995: +   x - a vector to mimic
1996: -   n - the number of vectors to obtain

1998:     Output Parameters:
1999: +   y - Fortran90 pointer to the array of vectors
2000: -   ierr - error code

2002:     Example of Usage:
2003: .vb
2004: #include <petsc/finclude/petscvec.h90>

2006:     Vec x
2007:     Vec, pointer :: y(:)
2008:     ....
2009:     call VecDuplicateVecsF90(x,2,y,ierr)
2010:     call VecSet(y(2),alpha,ierr)
2011:     call VecSet(y(2),alpha,ierr)
2012:     ....
2013:     call VecDestroyVecsF90(2,y,ierr)
2014: .ve

2016:     Notes:
2017:     Not yet supported for all F90 compilers

2019:     Use VecDestroyVecsF90() to free the space.

2021:     Level: beginner

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

2025: M*/

2027: /*MC
2028:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2029:     VecGetArrayF90().

2031:     Synopsis:
2032:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2034:     Logically Collective on Vec

2036:     Input Parameters:
2037: +   x - vector
2038: -   xx_v - the Fortran90 pointer to the array

2040:     Output Parameter:
2041: .   ierr - error code

2043:     Example of Usage:
2044: .vb
2045: #include <petsc/finclude/petscvec.h90>

2047:     PetscScalar, pointer :: xx_v(:)
2048:     ....
2049:     call VecGetArrayF90(x,xx_v,ierr)
2050:     xx_v(3) = a
2051:     call VecRestoreArrayF90(x,xx_v,ierr)
2052: .ve

2054:     Level: beginner

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

2058: M*/

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

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

2066:     Collective on Vec

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

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

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

2078:     Level: beginner

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

2082: M*/

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

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

2093:     Logically Collective on Vec

2095:     Input Parameter:
2096: .   x - vector

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

2102:     Example of Usage:
2103: .vb
2104: #include <petsc/finclude/petscvec.h90>

2106:     PetscScalar, pointer :: xx_v(:)
2107:     ....
2108:     call VecGetArrayF90(x,xx_v,ierr)
2109:     xx_v(3) = a
2110:     call VecRestoreArrayF90(x,xx_v,ierr)
2111: .ve

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

2115:     Level: beginner

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

2119: M*/

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

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

2130:     Logically Collective on Vec

2132:     Input Parameter:
2133: .   x - vector

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

2139:     Example of Usage:
2140: .vb
2141: #include <petsc/finclude/petscvec.h90>

2143:     PetscScalar, pointer :: xx_v(:)
2144:     ....
2145:     call VecGetArrayReadF90(x,xx_v,ierr)
2146:     a = xx_v(3)
2147:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2148: .ve

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

2152:     Level: beginner

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

2156: M*/

2158: /*MC
2159:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2160:     VecGetArrayReadF90().

2162:     Synopsis:
2163:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2165:     Logically Collective on Vec

2167:     Input Parameters:
2168: +   x - vector
2169: -   xx_v - the Fortran90 pointer to the array

2171:     Output Parameter:
2172: .   ierr - error code

2174:     Example of Usage:
2175: .vb
2176: #include <petsc/finclude/petscvec.h90>

2178:     PetscScalar, pointer :: xx_v(:)
2179:     ....
2180:     call VecGetArrayReadF90(x,xx_v,ierr)
2181:     a = xx_v(3)
2182:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2183: .ve

2185:     Level: beginner

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

2189: M*/

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

2198:    Logically Collective

2200:    Input Parameter:
2201: +  x - the vector
2202: .  m - first dimension of two dimensional array
2203: .  n - second dimension of two dimensional array
2204: .  mstart - first index you will use in first coordinate direction (often 0)
2205: -  nstart - first index in the second coordinate direction (often 0)

2207:    Output Parameter:
2208: .  a - location to put pointer to the array

2210:    Level: developer

2212:   Notes:
2213:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2214:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2215:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2216:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

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

2236:   VecGetLocalSize(x,&N);
2237:   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);
2238:   VecGetArray(x,&aa);

2240:   PetscMalloc1(m,a);
2241:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2242:   *a -= mstart;
2243:   return(0);
2244: }

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

2251:    Logically Collective

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

2261:    Level: developer

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

2269:    This routine actually zeros out the a pointer.

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

2284:   dummy = (void*)(*a + mstart);
2285:   PetscFree(dummy);
2286:   VecRestoreArray(x,NULL);
2287:   return(0);
2288: }

2292: /*@C
2293:    VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2294:    processor's portion of the vector data.  You MUST call VecRestoreArray1d()
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 two dimensional array
2302: -  mstart - first index you will use in first coordinate direction (often 0)

2304:    Output Parameter:
2305: .  a - location to put pointer to the array

2307:    Level: developer

2309:   Notes:
2310:    For a vector obtained from DMCreateLocalVector() mstart are likely
2311:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2312:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2316: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2317:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2318:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2319: @*/
2320: PetscErrorCode  VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2321: {
2323:   PetscInt       N;

2329:   VecGetLocalSize(x,&N);
2330:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2331:   VecGetArray(x,a);
2332:   *a  -= mstart;
2333:   return(0);
2334: }

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

2341:    Logically Collective

2343:    Input Parameters:
2344: +  x - the vector
2345: .  m - first dimension of two dimensional array
2346: .  mstart - first index you will use in first coordinate direction (often 0)
2347: -  a - location of pointer to array obtained from VecGetArray21()

2349:    Level: developer

2351:    Notes:
2352:    For regular PETSc vectors this routine does not involve any copies. For
2353:    any special vectors that do not store local vector data in a contiguous
2354:    array, this routine will copy the data back into the underlying
2355:    vector data structure from the array obtained with VecGetArray1d().

2357:    This routine actually zeros out the a pointer.

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

2361: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2362:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2363:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2364: @*/
2365: PetscErrorCode  VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2366: {

2372:   VecRestoreArray(x,NULL);
2373:   return(0);
2374: }


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

2384:    Logically Collective

2386:    Input Parameter:
2387: +  x - the vector
2388: .  m - first dimension of three dimensional array
2389: .  n - second dimension of three dimensional array
2390: .  p - third dimension of three dimensional array
2391: .  mstart - first index you will use in first coordinate direction (often 0)
2392: .  nstart - first index in the second coordinate direction (often 0)
2393: -  pstart - first index in the third coordinate direction (often 0)

2395:    Output Parameter:
2396: .  a - location to put pointer to the array

2398:    Level: developer

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

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

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

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

2424:   VecGetLocalSize(x,&N);
2425:   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);
2426:   VecGetArray(x,&aa);

2428:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
2429:   b    = (PetscScalar**)((*a) + m);
2430:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2431:   for (i=0; i<m; i++)
2432:     for (j=0; j<n; j++)
2433:       b[i*n+j] = aa + i*n*p + j*p - pstart;

2435:   *a -= mstart;
2436:   return(0);
2437: }

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

2444:    Logically Collective

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

2456:    Level: developer

2458:    Notes:
2459:    For regular PETSc vectors this routine does not involve any copies. For
2460:    any special vectors that do not store local vector data in a contiguous
2461:    array, this routine will copy the data back into the underlying
2462:    vector data structure from the array obtained with VecGetArray().

2464:    This routine actually zeros out the a pointer.

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

2479:   dummy = (void*)(*a + mstart);
2480:   PetscFree(dummy);
2481:   VecRestoreArray(x,NULL);
2482:   return(0);
2483: }

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

2492:    Logically Collective

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

2505:    Output Parameter:
2506: .  a - location to put pointer to the array

2508:    Level: beginner

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

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

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

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

2534:   VecGetLocalSize(x,&N);
2535:   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);
2536:   VecGetArray(x,&aa);

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

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

2558:    Logically Collective

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

2572:    Level: beginner

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

2580:    This routine actually zeros out the a pointer.

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

2595:   dummy = (void*)(*a + mstart);
2596:   PetscFree(dummy);
2597:   VecRestoreArray(x,NULL);
2598:   return(0);
2599: }

2603: /*@C
2604:    VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
2605:    processor's portion of the vector data.  You MUST call VecRestoreArray2dRead()
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: .  n - second dimension of two dimensional array
2614: .  mstart - first index you will use in first coordinate direction (often 0)
2615: -  nstart - first index in the second coordinate direction (often 0)

2617:    Output Parameter:
2618: .  a - location to put pointer to the array

2620:    Level: developer

2622:   Notes:
2623:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2624:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2625:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2626:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

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

2646:   VecGetLocalSize(x,&N);
2647:   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);
2648:   VecGetArrayRead(x,&aa);

2650:   PetscMalloc1(m,a);
2651:   for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
2652:   *a -= mstart;
2653:   return(0);
2654: }

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

2661:    Logically Collective

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

2671:    Level: developer

2673:    Notes:
2674:    For regular PETSc vectors this routine does not involve any copies. For
2675:    any special vectors that do not store local vector data in a contiguous
2676:    array, this routine will copy the data back into the underlying
2677:    vector data structure from the array obtained with VecGetArray().

2679:    This routine actually zeros out the a pointer.

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

2694:   dummy = (void*)(*a + mstart);
2695:   PetscFree(dummy);
2696:   VecRestoreArrayRead(x,NULL);
2697:   return(0);
2698: }

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

2707:    Logically Collective

2709:    Input Parameter:
2710: +  x - the vector
2711: .  m - first dimension of two dimensional array
2712: -  mstart - first index you will use in first coordinate direction (often 0)

2714:    Output Parameter:
2715: .  a - location to put pointer to the array

2717:    Level: developer

2719:   Notes:
2720:    For a vector obtained from DMCreateLocalVector() mstart are likely
2721:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2722:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2726: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2727:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2728:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2729: @*/
2730: PetscErrorCode  VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2731: {
2733:   PetscInt       N;

2739:   VecGetLocalSize(x,&N);
2740:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2741:   VecGetArrayRead(x,(const PetscScalar**)a);
2742:   *a  -= mstart;
2743:   return(0);
2744: }

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

2751:    Logically Collective

2753:    Input Parameters:
2754: +  x - the vector
2755: .  m - first dimension of two dimensional array
2756: .  mstart - first index you will use in first coordinate direction (often 0)
2757: -  a - location of pointer to array obtained from VecGetArray21()

2759:    Level: developer

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

2767:    This routine actually zeros out the a pointer.

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

2771: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2772:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2773:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2774: @*/
2775: PetscErrorCode  VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2776: {

2782:   VecRestoreArrayRead(x,NULL);
2783:   return(0);
2784: }


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

2794:    Logically Collective

2796:    Input Parameter:
2797: +  x - the vector
2798: .  m - first dimension of three dimensional array
2799: .  n - second dimension of three dimensional array
2800: .  p - third dimension of three dimensional array
2801: .  mstart - first index you will use in first coordinate direction (often 0)
2802: .  nstart - first index in the second coordinate direction (often 0)
2803: -  pstart - first index in the third coordinate direction (often 0)

2805:    Output Parameter:
2806: .  a - location to put pointer to the array

2808:    Level: developer

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

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

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

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

2835:   VecGetLocalSize(x,&N);
2836:   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);
2837:   VecGetArrayRead(x,&aa);

2839:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
2840:   b    = (PetscScalar**)((*a) + m);
2841:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2842:   for (i=0; i<m; i++)
2843:     for (j=0; j<n; j++)
2844:       b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;

2846:   *a -= mstart;
2847:   return(0);
2848: }

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

2855:    Logically Collective

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

2867:    Level: developer

2869:    Notes:
2870:    For regular PETSc vectors this routine does not involve any copies. For
2871:    any special vectors that do not store local vector data in a contiguous
2872:    array, this routine will copy the data back into the underlying
2873:    vector data structure from the array obtained with VecGetArray().

2875:    This routine actually zeros out the a pointer.

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

2890:   dummy = (void*)(*a + mstart);
2891:   PetscFree(dummy);
2892:   VecRestoreArrayRead(x,NULL);
2893:   return(0);
2894: }

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

2903:    Logically Collective

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

2916:    Output Parameter:
2917: .  a - location to put pointer to the array

2919:    Level: beginner

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

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

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

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

2946:   VecGetLocalSize(x,&N);
2947:   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);
2948:   VecGetArrayRead(x,&aa);

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

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

2970:    Logically Collective

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

2984:    Level: beginner

2986:    Notes:
2987:    For regular PETSc vectors this routine does not involve any copies. For
2988:    any special vectors that do not store local vector data in a contiguous
2989:    array, this routine will copy the data back into the underlying
2990:    vector data structure from the array obtained with VecGetArray().

2992:    This routine actually zeros out the a pointer.

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

3007:   dummy = (void*)(*a + mstart);
3008:   PetscFree(dummy);
3009:   VecRestoreArrayRead(x,NULL);
3010:   return(0);
3011: }

3013: #if defined(PETSC_USE_DEBUG)

3017: /*@
3018:    VecLockGet  - Gets the current lock status of a vector

3020:    Logically Collective on Vec

3022:    Input Parameter:
3023: .  x - the vector

3025:    Output Parameter:
3026: .  state - greater than zero indicates the vector is still locked

3028:    Level: beginner

3030:    Concepts: vector^accessing local values

3032: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockPush(), VecLockGet()
3033: @*/
3034: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
3035: {
3038:   *state = x->lock;
3039:   return(0);
3040: }

3044: /*@
3045:    VecLockPush  - Lock a vector from writing

3047:    Logically Collective on Vec

3049:    Input Parameter:
3050: .  x - the vector

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

3054:     Call VecLockPop() to remove the latest lock

3056:    Level: beginner

3058:    Concepts: vector^accessing local values

3060: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockPop(), VecLockGet()
3061: @*/
3062: PetscErrorCode VecLockPush(Vec x)
3063: {
3066:   x->lock++;
3067:   return(0);
3068: }

3072: /*@
3073:    VecLockPop  - Unlock a vector from writing

3075:    Logically Collective on Vec

3077:    Input Parameter:
3078: .  x - the vector

3080:    Level: beginner

3082:    Concepts: vector^accessing local values

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

3095: #endif