Actual source code: rvector.c

petsc-3.11.2 2019-05-18
Report Typos and Errors

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

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

 17: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL)
 18:   if ((vec->petscnative || vec->ops->getarray) && (vec->valid_GPU_array == PETSC_OFFLOAD_CPU || vec->valid_GPU_array == PETSC_OFFLOAD_BOTH)) {
 19: #else
 20:   if (vec->petscnative || vec->ops->getarray) {
 21: #endif
 22:     VecGetLocalSize(vec,&n);
 23:     VecGetArrayRead(vec,&x);
 24:     for (i=0; i<n; i++) {
 25:       if (begin) {
 26:         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);
 27:       } else {
 28:         if (PetscIsInfOrNanScalar(x[i])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FP,"Vec entry at local location %D is not-a-number or infinite at end of function: Parameter number %D",i,argnum);
 29:       }
 30:     }
 31:     VecRestoreArrayRead(vec,&x);
 32:   }
 33:   return(0);
 34: #else
 35:   return 0;
 36: #endif
 37: }

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

 42:    Logically Collective on Vec

 44:    Input Parameters:
 45: .  x, y  - the vectors

 47:    Output Parameter:
 48: .  max - the result

 50:    Level: advanced

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

 56: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs()
 57: @*/
 58: PetscErrorCode  VecMaxPointwiseDivide(Vec x,Vec y,PetscReal *max)
 59: {

 69:   VecCheckSameSize(x,1,y,2);
 70:   (*x->ops->maxpointwisedivide)(x,y,max);
 71:   return(0);
 72: }

 74: /*@
 75:    VecDot - Computes the vector dot product.

 77:    Collective on Vec

 79:    Input Parameters:
 80: .  x, y - the vectors

 82:    Output Parameter:
 83: .  val - the dot product

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

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

 97:    Use VecTDot() for the indefinite form
 98: $     val = (x,y) = y^T x,
 99:    where y^T denotes the transpose of y.

101:    Level: intermediate

103:    Concepts: inner product
104:    Concepts: vector^inner product

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

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

121:   PetscLogEventBegin(VEC_Dot,x,y,0,0);
122:   (*x->ops->dot)(x,y,val);
123:   PetscLogEventEnd(VEC_Dot,x,y,0,0);
124:   return(0);
125: }

127: /*@
128:    VecDotRealPart - Computes the real part of the vector dot product.

130:    Collective on Vec

132:    Input Parameters:
133: .  x, y - the vectors

135:    Output Parameter:
136: .  val - the real part of the dot product;

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

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

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

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

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

153:    Level: intermediate

155:    Concepts: inner product
156:    Concepts: vector^inner product

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

166:   VecDot(x,y,&fdot);
167:   *val = PetscRealPart(fdot);
168:   return(0);
169: }

171: /*@
172:    VecNorm  - Computes the vector norm.

174:    Collective on Vec

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

182:    Output Parameter:
183: .  val - the norm

185:    Notes:
186: $     NORM_1 denotes sum_i |x_i|
187: $     NORM_2 denotes sqrt(sum_i |x_i|^2)
188: $     NORM_INFINITY denotes max_i |x_i|

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

195:    Level: intermediate

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

202:    Concepts: norm
203:    Concepts: vector^norm

205: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNormAvailable(),
206:           VecNormBegin(), VecNormEnd()

208: @*/
209: PetscErrorCode  VecNorm(Vec x,NormType type,PetscReal *val)
210: {
211:   PetscBool      flg;


219:   /*
220:    * Cached data?
221:    */
222:   if (type!=NORM_1_AND_2) {
223:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,flg);
224:     if (flg) return(0);
225:   }
226:   PetscLogEventBegin(VEC_Norm,x,0,0,0);
227:   (*x->ops->norm)(x,type,val);
228:   PetscLogEventEnd(VEC_Norm,x,0,0,0);

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

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

239:    Not Collective

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

247:    Output Parameter:
248: +  available - PETSC_TRUE if the val returned is valid
249: -  val - the norm

251:    Notes:
252: $     NORM_1 denotes sum_i |x_i|
253: $     NORM_2 denotes sqrt(sum_i (x_i)^2)
254: $     NORM_INFINITY denotes max_i |x_i|

256:    Level: intermediate

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

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

268:    Concepts: norm
269:    Concepts: vector^norm

271: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNorm()
272:           VecNormBegin(), VecNormEnd()

274: @*/
275: PetscErrorCode  VecNormAvailable(Vec x,NormType type,PetscBool  *available,PetscReal *val)
276: {


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

291: /*@
292:    VecNormalize - Normalizes a vector by 2-norm.

294:    Collective on Vec

296:    Input Parameters:
297: +  x - the vector

299:    Output Parameter:
300: .  x - the normalized vector
301: -  val - the vector norm before normalization

303:    Level: intermediate

305:    Concepts: vector^normalizing
306:    Concepts: normalizing^vector

308: @*/
309: PetscErrorCode  VecNormalize(Vec x,PetscReal *val)
310: {
312:   PetscReal      norm;

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

330: /*@C
331:    VecMax - Determines the vector component with maximum real part and its location.

333:    Collective on Vec

335:    Input Parameter:
336: .  x - the vector

338:    Output Parameters:
339: +  p - the location of val (pass NULL if you don't want this)
340: -  val - the maximum component

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

345:    Returns the smallest index with the maximum value
346:    Level: intermediate

348:    Concepts: maximum^of vector
349:    Concepts: vector^maximum value

351: .seealso: VecNorm(), VecMin()
352: @*/
353: PetscErrorCode  VecMax(Vec x,PetscInt *p,PetscReal *val)
354: {

361:   PetscLogEventBegin(VEC_Max,x,0,0,0);
362:   (*x->ops->max)(x,p,val);
363:   PetscLogEventEnd(VEC_Max,x,0,0,0);
364:   return(0);
365: }

367: /*@C
368:    VecMin - Determines the vector component with minimum real part and its location.

370:    Collective on Vec

372:    Input Parameters:
373: .  x - the vector

375:    Output Parameter:
376: +  p - the location of val (pass NULL if you don't want this location)
377: -  val - the minimum component

379:    Level: intermediate

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

384:    This returns the smallest index with the minumum value

386:    Concepts: minimum^of vector
387:    Concepts: vector^minimum entry

389: .seealso: VecMax()
390: @*/
391: PetscErrorCode  VecMin(Vec x,PetscInt *p,PetscReal *val)
392: {

399:   PetscLogEventBegin(VEC_Min,x,0,0,0);
400:   (*x->ops->min)(x,p,val);
401:   PetscLogEventEnd(VEC_Min,x,0,0,0);
402:   return(0);
403: }

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

409:    Collective on Vec

411:    Input Parameters:
412: .  x, y - the vectors

414:    Output Parameter:
415: .  val - the dot product

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

422:    Use VecDot() for the inner product
423: $     val = (x,y) = y^H x,
424:    where y^H denotes the conjugate transpose of y.

426:    Level: intermediate

428:    Concepts: inner product^non-Hermitian
429:    Concepts: vector^inner product
430:    Concepts: non-Hermitian inner product

432: .seealso: VecDot(), VecMTDot()
433: @*/
434: PetscErrorCode  VecTDot(Vec x,Vec y,PetscScalar *val)
435: {

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

447:   PetscLogEventBegin(VEC_TDot,x,y,0,0);
448:   (*x->ops->tdot)(x,y,val);
449:   PetscLogEventEnd(VEC_TDot,x,y,0,0);
450:   return(0);
451: }

453: /*@
454:    VecScale - Scales a vector.

456:    Not collective on Vec

458:    Input Parameters:
459: +  x - the vector
460: -  alpha - the scalar

462:    Output Parameter:
463: .  x - the scaled vector

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

469:    Level: intermediate

471:    Concepts: vector^scaling
472:    Concepts: scaling^vector

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

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

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

508:    Logically Collective on Vec

510:    Input Parameters:
511: +  x  - the vector
512: -  alpha - the scalar

514:    Output Parameter:
515: .  x  - the vector

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

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

527:    Level: beginner

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

531:    Concepts: vector^setting to constant

533: @*/
534: PetscErrorCode  VecSet(Vec x,PetscScalar alpha)
535: {
536:   PetscReal      val;

542:   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()");
544:   VecSetErrorIfLocked(x,1);

546:   PetscLogEventBegin(VEC_Set,x,0,0,0);
547:   (*x->ops->set)(x,alpha);
548:   PetscLogEventEnd(VEC_Set,x,0,0,0);
549:   PetscObjectStateIncrease((PetscObject)x);

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


571: /*@
572:    VecAXPY - Computes y = alpha x + y.

574:    Logically Collective on Vec

576:    Input Parameters:
577: +  alpha - the scalar
578: -  x, y  - the vectors

580:    Output Parameter:
581: .  y - output vector

583:    Level: intermediate

585:    Notes:
586:     x and y MUST be different vectors

588:    Concepts: vector^BLAS
589:    Concepts: BLAS

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

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

608:   VecLockReadPush(x);
609:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
610:   (*y->ops->axpy)(y,alpha,x);
611:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
612:   VecLockReadPop(x);
613:   PetscObjectStateIncrease((PetscObject)y);
614:   return(0);
615: }

617: /*@
618:    VecAXPBY - Computes y = alpha x + beta y.

620:    Logically Collective on Vec

622:    Input Parameters:
623: +  alpha,beta - the scalars
624: -  x, y  - the vectors

626:    Output Parameter:
627: .  y - output vector

629:    Level: intermediate

631:    Notes:
632:     x and y MUST be different vectors

634:    Concepts: BLAS
635:    Concepts: vector^BLAS

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

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

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

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

664:    Logically Collective on Vec

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

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

673:    Level: intermediate

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

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

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

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

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

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

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

716:    Logically Collective on Vec

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

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

725:    Level: intermediate

727:    Notes:
728:     x and y MUST be different vectors

730:    Concepts: vector^BLAS
731:    Concepts: BLAS

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

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

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


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

760:    Logically Collective on Vec

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

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

769:    Level: intermediate

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

774:    Concepts: vector^BLAS
775:    Concepts: BLAS

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

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

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


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

809:    Not Collective

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

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

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

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

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

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

838:    Level: beginner

840:    Concepts: vector^setting values

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

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

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

866:     Not Collective

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

873:    Output Parameter:
874: .   y - array of values

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

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

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

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

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

889:    Level: beginner

891:    Concepts: vector^getting values

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

901:   if (!ni) return(0);
905:   (*x->ops->getvalues)(x,ni,ix,y);
906:   return(0);
907: }

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

912:    Not Collective

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

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

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

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

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

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

941:    Level: intermediate

943:    Concepts: vector^setting values blocked

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

954:   if (!ni) return(0);
958:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
959:   (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
960:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
961:   PetscObjectStateIncrease((PetscObject)x);
962:   return(0);
963: }


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

970:    Not Collective

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

981:    Level: intermediate

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

986:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
987:    options cannot be mixed without intervening calls to the assembly
988:    routines.

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

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

995:    Concepts: vector^setting values with local numbering

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

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

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

1031: /*@
1032:    VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1033:    using a local ordering of the nodes.

1035:    Not Collective

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

1046:    Level: intermediate

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

1052:    Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1053:    options cannot be mixed without intervening calls to the assembly
1054:    routines.

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

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


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

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

1074:   if (!ni) return(0);
1078:   if (ni > 128) {
1079:     PetscMalloc1(ni,&lix);
1080:   }

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

1093: /*@
1094:    VecMTDot - Computes indefinite vector multiple dot products.
1095:    That is, it does NOT use the complex conjugate.

1097:    Collective on Vec

1099:    Input Parameters:
1100: +  x - one vector
1101: .  nv - number of vectors
1102: -  y - array of vectors.  Note that vectors are pointers

1104:    Output Parameter:
1105: .  val - array of the dot products

1107:    Notes for Users of Complex Numbers:
1108:    For complex vectors, VecMTDot() computes the indefinite form
1109: $      val = (x,y) = y^T x,
1110:    where y^T denotes the transpose of y.

1112:    Use VecMDot() for the inner product
1113: $      val = (x,y) = y^H x,
1114:    where y^H denotes the conjugate transpose of y.

1116:    Level: intermediate

1118:    Concepts: inner product^multiple
1119:    Concepts: vector^multiple inner products

1121: .seealso: VecMDot(), VecTDot()
1122: @*/
1123: PetscErrorCode  VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1124: {

1130:   if (!nv) return(0);
1137:   VecCheckSameSize(x,1,*y,3);

1139:   PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1140:   (*x->ops->mtdot)(x,nv,y,val);
1141:   PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1142:   return(0);
1143: }

1145: /*@
1146:    VecMDot - Computes vector multiple dot products.

1148:    Collective on Vec

1150:    Input Parameters:
1151: +  x - one vector
1152: .  nv - number of vectors
1153: -  y - array of vectors.

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

1158:    Notes for Users of Complex Numbers:
1159:    For complex vectors, VecMDot() computes
1160: $     val = (x,y) = y^H x,
1161:    where y^H denotes the conjugate transpose of y.

1163:    Use VecMTDot() for the indefinite form
1164: $     val = (x,y) = y^T x,
1165:    where y^T denotes the transpose of y.

1167:    Level: intermediate

1169:    Concepts: inner product^multiple
1170:    Concepts: vector^multiple inner products

1172: .seealso: VecMTDot(), VecDot()
1173: @*/
1174: PetscErrorCode  VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1175: {

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

1191:   PetscLogEventBegin(VEC_MDot,x,*y,0,0);
1192:   (*x->ops->mdot)(x,nv,y,val);
1193:   PetscLogEventEnd(VEC_MDot,x,*y,0,0);
1194:   return(0);
1195: }

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

1200:    Logically Collective on Vec

1202:    Input Parameters:
1203: +  nv - number of scalars and x-vectors
1204: .  alpha - array of scalars
1205: .  y - one vector
1206: -  x - array of vectors

1208:    Level: intermediate

1210:    Notes:
1211:     y cannot be any of the x vectors

1213:    Concepts: BLAS

1215: .seealso: VecAXPY(), VecWAXPY(), VecAYPX()
1216: @*/
1217: PetscErrorCode  VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1218: {
1220:   PetscInt       i;

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

1236:   PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1237:   (*y->ops->maxpy)(y,nv,alpha,x);
1238:   PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1239:   PetscObjectStateIncrease((PetscObject)y);
1240:   return(0);
1241: }

1243: /*@
1244:    VecGetSubVector - Gets a vector representing part of another vector

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

1248:    Input Arguments:
1249: + X - vector from which to extract a subvector
1250: - is - index set representing portion of X to extract

1252:    Output Arguments:
1253: . Y - subvector corresponding to is

1255:    Level: advanced

1257:    Notes:
1258:    The subvector Y should be returned with VecRestoreSubVector().

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

1263: .seealso: MatCreateSubMatrix()
1264: @*/
1265: PetscErrorCode  VecGetSubVector(Vec X,IS is,Vec *Y)
1266: {
1267:   PetscErrorCode   ierr;
1268:   Vec              Z;

1274:   if (X->ops->getsubvector) {
1275:     (*X->ops->getsubvector)(X,is,&Z);
1276:   } else { /* Default implementation currently does no caching */
1277:     PetscInt  gstart,gend,start;
1278:     PetscBool contiguous,gcontiguous;
1279:     VecGetOwnershipRange(X,&gstart,&gend);
1280:     ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1281:     MPIU_Allreduce(&contiguous,&gcontiguous,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1282:     if (gcontiguous) { /* We can do a no-copy implementation */
1283:       PetscInt n,N,bs;
1284:       PetscInt state;

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

1333: /*@
1334:    VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()

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

1338:    Input Arguments:
1339: + X - vector from which subvector was obtained
1340: . is - index set representing the subset of X
1341: - Y - subvector being restored

1343:    Level: advanced

1345: .seealso: VecGetSubVector()
1346: @*/
1347: PetscErrorCode  VecRestoreSubVector(Vec X,IS is,Vec *Y)
1348: {

1356:   if (X->ops->restoresubvector) {
1357:     (*X->ops->restoresubvector)(X,is,Y);
1358:   } else {
1359:     PETSC_UNUSED PetscObjectState dummystate = 0;
1360:     PetscBool valid;
1361:     PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,dummystate,valid);
1362:     if (!valid) {
1363:       VecScatter scatter;

1365:       PetscObjectQuery((PetscObject)*Y,"VecGetSubVector_Scatter",(PetscObject*)&scatter);
1366:       if (scatter) {
1367:         VecScatterBegin(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1368:         VecScatterEnd(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1369:       }
1370:     }
1371:     VecDestroy(Y);
1372:   }
1373:   return(0);
1374: }

1376: /*@
1377:    VecGetLocalVectorRead - Maps the local portion of a vector into a
1378:    vector.  You must call VecRestoreLocalVectorRead() when the local
1379:    vector is no longer needed.

1381:    Not collective.

1383:    Input parameter:
1384: .  v - The vector for which the local vector is desired.

1386:    Output parameter:
1387: .  w - Upon exit this contains the local vector.

1389:    Level: beginner

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

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

1405: .seealso: VecRestoreLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1406: @*/
1407: PetscErrorCode VecGetLocalVectorRead(Vec v,Vec w)
1408: {
1410:   PetscScalar    *a;

1415:   VecCheckSameLocalSize(v,1,w,2);
1416:   if (v->ops->getlocalvectorread) {
1417:     (*v->ops->getlocalvectorread)(v,w);
1418:   } else {
1419:     VecGetArrayRead(v,(const PetscScalar**)&a);
1420:     VecPlaceArray(w,a);
1421:   }
1422:   return(0);
1423: }

1425: /*@
1426:    VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1427:    previously mapped into a vector using VecGetLocalVectorRead().

1429:    Not collective.

1431:    Input parameter:
1432: .  v - The local portion of this vector was previously mapped into w using VecGetLocalVectorRead().
1433: .  w - The vector into which the local portion of v was mapped.

1435:    Level: beginner

1437: .seealso: VecGetLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1438: @*/
1439: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec w)
1440: {
1442:   PetscScalar    *a;

1447:   if (v->ops->restorelocalvectorread) {
1448:     (*v->ops->restorelocalvectorread)(v,w);
1449:   } else {
1450:     VecGetArrayRead(w,(const PetscScalar**)&a);
1451:     VecRestoreArrayRead(v,(const PetscScalar**)&a);
1452:     VecResetArray(w);
1453:   }
1454:   return(0);
1455: }

1457: /*@
1458:    VecGetLocalVector - Maps the local portion of a vector into a
1459:    vector.

1461:    Collective on v, not collective on w.

1463:    Input parameter:
1464: .  v - The vector for which the local vector is desired.

1466:    Output parameter:
1467: .  w - Upon exit this contains the local vector.

1469:    Level: beginner

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

1482: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1483: @*/
1484: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1485: {
1487:   PetscScalar    *a;

1492:   VecCheckSameLocalSize(v,1,w,2);
1493:   if (v->ops->getlocalvector) {
1494:     (*v->ops->getlocalvector)(v,w);
1495:   } else {
1496:     VecGetArray(v,&a);
1497:     VecPlaceArray(w,a);
1498:   }
1499:   return(0);
1500: }

1502: /*@
1503:    VecRestoreLocalVector - Unmaps the local portion of a vector
1504:    previously mapped into a vector using VecGetLocalVector().

1506:    Logically collective.

1508:    Input parameter:
1509: .  v - The local portion of this vector was previously mapped into w using VecGetLocalVector().
1510: .  w - The vector into which the local portion of v was mapped.

1512:    Level: beginner

1514: .seealso: VecGetLocalVector(), VecGetLocalVectorRead(), VecRestoreLocalVectorRead(), LocalVectorRead(), VecGetArrayRead(), VecGetArray()
1515: @*/
1516: PetscErrorCode VecRestoreLocalVector(Vec v,Vec w)
1517: {
1519:   PetscScalar    *a;

1524:   if (v->ops->restorelocalvector) {
1525:     (*v->ops->restorelocalvector)(v,w);
1526:   } else {
1527:     VecGetArray(w,&a);
1528:     VecRestoreArray(v,&a);
1529:     VecResetArray(w);
1530:   }
1531:   return(0);
1532: }

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

1543:    Logically Collective on Vec

1545:    Input Parameter:
1546: .  x - the vector

1548:    Output Parameter:
1549: .  a - location to put pointer to the array

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

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

1569:    Level: beginner

1571:    Concepts: vector^accessing local values

1573: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1574:           VecGetArrayPair(), VecRestoreArrayPair()
1575: @*/
1576: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1577: {
1579: #if defined(PETSC_HAVE_VIENNACL)
1580:   PetscBool      is_viennacltype = PETSC_FALSE;
1581: #endif

1585:   VecSetErrorIfLocked(x,1);
1586:   if (x->petscnative) {
1587: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1588:     if (x->valid_GPU_array == PETSC_OFFLOAD_GPU) {
1589: #if defined(PETSC_HAVE_VIENNACL)
1590:       PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1591:       if (is_viennacltype) {
1592:         VecViennaCLCopyFromGPU(x);
1593:       } else
1594: #endif
1595:       {
1596: #if defined(PETSC_HAVE_CUDA)
1597:         VecCUDACopyFromGPU(x);
1598: #endif
1599:       }
1600:     } else if (x->valid_GPU_array == PETSC_OFFLOAD_UNALLOCATED) {
1601: #if defined(PETSC_HAVE_VIENNACL)
1602:       PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1603:       if (is_viennacltype) {
1604:         VecViennaCLAllocateCheckHost(x);
1605:       } else
1606: #endif
1607:       {
1608: #if defined(PETSC_HAVE_CUDA)
1609:         VecCUDAAllocateCheckHost(x);
1610: #endif
1611:       }
1612:     }
1613: #endif
1614:     *a = *((PetscScalar**)x->data);
1615:   } else {
1616:     if (x->ops->getarray) {
1617:       (*x->ops->getarray)(x,a);
1618:     } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array for vector type \"%s\"",((PetscObject)x)->type_name);
1619:   }
1620:   return(0);
1621: }

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

1626:    Not Collective

1628:    Input Parameters:
1629: .  x - the vector

1631:    Output Parameter:
1632: .  a - the array

1634:    Level: beginner

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

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

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

1645: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1646: @*/
1647: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1648: {
1650: #if defined(PETSC_HAVE_VIENNACL)
1651:   PetscBool      is_viennacltype = PETSC_FALSE;
1652: #endif

1656:   if (x->petscnative) {
1657: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1658:     if (x->valid_GPU_array == PETSC_OFFLOAD_GPU) {
1659: #if defined(PETSC_HAVE_VIENNACL)
1660:       PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1661:       if (is_viennacltype) {
1662:         VecViennaCLCopyFromGPU(x);
1663:       } else
1664: #endif
1665:       {
1666: #if defined(PETSC_HAVE_CUDA)
1667:         VecCUDACopyFromGPU(x);
1668: #endif
1669:       }
1670:     }
1671: #endif
1672:     *a = *((PetscScalar **)x->data);
1673:   } else if (x->ops->getarrayread) {
1674:     (*x->ops->getarrayread)(x,a);
1675:   } else {
1676:     (*x->ops->getarray)(x,(PetscScalar**)a);
1677:   }
1678:   return(0);
1679: }

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

1686:    Logically Collective on Vec

1688:    Input Parameter:
1689: +  x - the vectors
1690: -  n - the number of vectors

1692:    Output Parameter:
1693: .  a - location to put pointer to the array

1695:    Fortran Note:
1696:    This routine is not supported in Fortran.

1698:    Level: intermediate

1700: .seealso: VecGetArray(), VecRestoreArrays()
1701: @*/
1702: PetscErrorCode  VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1703: {
1705:   PetscInt       i;
1706:   PetscScalar    **q;

1712:   if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1713:   PetscMalloc1(n,&q);
1714:   for (i=0; i<n; ++i) {
1715:     VecGetArray(x[i],&q[i]);
1716:   }
1717:   *a = q;
1718:   return(0);
1719: }

1721: /*@C
1722:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1723:    has been called.

1725:    Logically Collective on Vec

1727:    Input Parameters:
1728: +  x - the vector
1729: .  n - the number of vectors
1730: -  a - location of pointer to arrays obtained from VecGetArrays()

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

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

1741:    Level: intermediate

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


1756:   for (i=0; i<n; ++i) {
1757:     VecRestoreArray(x[i],&q[i]);
1758:   }
1759:   PetscFree(q);
1760:   return(0);
1761: }

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

1766:    Logically Collective on Vec

1768:    Input Parameters:
1769: +  x - the vector
1770: -  a - location of pointer to array obtained from VecGetArray()

1772:    Level: beginner

1774:    Notes:
1775:    For regular PETSc vectors this routine does not involve any copies. For
1776:    any special vectors that do not store local vector data in a contiguous
1777:    array, this routine will copy the data back into the underlying
1778:    vector data structure from the array obtained with VecGetArray().

1780:    This routine actually zeros out the a pointer. This is to prevent accidental
1781:    use of the array after it has been restored. If you pass null for a it will
1782:    not zero the array pointer a.

1784:    Fortran Note:
1785:    This routine is used differently from Fortran 77
1786: $    Vec         x
1787: $    PetscScalar x_array(1)
1788: $    PetscOffset i_x
1789: $    PetscErrorCode ierr
1790: $       call VecGetArray(x,x_array,i_x,ierr)
1791: $
1792: $   Access first local entry in vector with
1793: $      value = x_array(i_x + 1)
1794: $
1795: $      ...... other code
1796: $       call VecRestoreArray(x,x_array,i_x,ierr)

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

1802: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1803:           VecGetArrayPair(), VecRestoreArrayPair()
1804: @*/
1805: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1806: {

1811:   if (x->petscnative) {
1812: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1813:     x->valid_GPU_array = PETSC_OFFLOAD_CPU;
1814: #endif
1815:   } else {
1816:     (*x->ops->restorearray)(x,a);
1817:   }
1818:   if (a) *a = NULL;
1819:   PetscObjectStateIncrease((PetscObject)x);
1820:   return(0);
1821: }

1823: /*@C
1824:    VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()

1826:    Not Collective

1828:    Input Parameters:
1829: +  vec - the vector
1830: -  array - the array

1832:    Level: beginner

1834: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1835: @*/
1836: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1837: {

1842:   if (x->petscnative) {
1843:     /* nothing */
1844:   } else if (x->ops->restorearrayread) {
1845:     (*x->ops->restorearrayread)(x,a);
1846:   } else {
1847:     (*x->ops->restorearray)(x,(PetscScalar**)a);
1848:   }
1849:   if (a) *a = NULL;
1850:   return(0);
1851: }

1853: /*@
1854:    VecPlaceArray - Allows one to replace the array in a vector with an
1855:    array provided by the user. This is useful to avoid copying an array
1856:    into a vector.

1858:    Not Collective

1860:    Input Parameters:
1861: +  vec - the vector
1862: -  array - the array

1864:    Notes:
1865:    You can return to the original array with a call to VecResetArray()

1867:    Level: developer

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

1871: @*/
1872: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
1873: {

1880:   if (vec->ops->placearray) {
1881:     (*vec->ops->placearray)(vec,array);
1882:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
1883:   PetscObjectStateIncrease((PetscObject)vec);
1884:   return(0);
1885: }

1887: /*@C
1888:    VecReplaceArray - Allows one to replace the array in a vector with an
1889:    array provided by the user. This is useful to avoid copying an array
1890:    into a vector.

1892:    Not Collective

1894:    Input Parameters:
1895: +  vec - the vector
1896: -  array - the array

1898:    Notes:
1899:    This permanently replaces the array and frees the memory associated
1900:    with the old array.

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

1905:    Not supported from Fortran

1907:    Level: developer

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

1911: @*/
1912: PetscErrorCode  VecReplaceArray(Vec vec,const PetscScalar array[])
1913: {

1919:   if (vec->ops->replacearray) {
1920:     (*vec->ops->replacearray)(vec,array);
1921:   } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
1922:   PetscObjectStateIncrease((PetscObject)vec);
1923:   return(0);
1924: }

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

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

1933:     Collective on Vec

1935:     Input Parameters:
1936: +   x - a vector to mimic
1937: -   n - the number of vectors to obtain

1939:     Output Parameters:
1940: +   y - Fortran90 pointer to the array of vectors
1941: -   ierr - error code

1943:     Example of Usage:
1944: .vb
1945: #include <petsc/finclude/petscvec.h>
1946:     use petscvec

1948:     Vec x
1949:     Vec, pointer :: y(:)
1950:     ....
1951:     call VecDuplicateVecsF90(x,2,y,ierr)
1952:     call VecSet(y(2),alpha,ierr)
1953:     call VecSet(y(2),alpha,ierr)
1954:     ....
1955:     call VecDestroyVecsF90(2,y,ierr)
1956: .ve

1958:     Notes:
1959:     Not yet supported for all F90 compilers

1961:     Use VecDestroyVecsF90() to free the space.

1963:     Level: beginner

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

1967: M*/

1969: /*MC
1970:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
1971:     VecGetArrayF90().

1973:     Synopsis:
1974:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

1976:     Logically Collective on Vec

1978:     Input Parameters:
1979: +   x - vector
1980: -   xx_v - the Fortran90 pointer to the array

1982:     Output Parameter:
1983: .   ierr - error code

1985:     Example of Usage:
1986: .vb
1987: #include <petsc/finclude/petscvec.h>
1988:     use petscvec

1990:     PetscScalar, pointer :: xx_v(:)
1991:     ....
1992:     call VecGetArrayF90(x,xx_v,ierr)
1993:     xx_v(3) = a
1994:     call VecRestoreArrayF90(x,xx_v,ierr)
1995: .ve

1997:     Level: beginner

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

2001: M*/

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

2006:     Synopsis:
2007:     VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)

2009:     Collective on Vec

2011:     Input Parameters:
2012: +   n - the number of vectors previously obtained
2013: -   x - pointer to array of vector pointers

2015:     Output Parameter:
2016: .   ierr - error code

2018:     Notes:
2019:     Not yet supported for all F90 compilers

2021:     Level: beginner

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

2025: M*/

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

2033:     Synopsis:
2034:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2036:     Logically Collective on Vec

2038:     Input Parameter:
2039: .   x - vector

2041:     Output Parameters:
2042: +   xx_v - the Fortran90 pointer to the array
2043: -   ierr - error code

2045:     Example of Usage:
2046: .vb
2047: #include <petsc/finclude/petscvec.h>
2048:     use petscvec

2050:     PetscScalar, pointer :: xx_v(:)
2051:     ....
2052:     call VecGetArrayF90(x,xx_v,ierr)
2053:     xx_v(3) = a
2054:     call VecRestoreArrayF90(x,xx_v,ierr)
2055: .ve

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

2059:     Level: beginner

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

2063: M*/

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

2071:     Synopsis:
2072:     VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2074:     Logically Collective on Vec

2076:     Input Parameter:
2077: .   x - vector

2079:     Output Parameters:
2080: +   xx_v - the Fortran90 pointer to the array
2081: -   ierr - error code

2083:     Example of Usage:
2084: .vb
2085: #include <petsc/finclude/petscvec.h>
2086:     use petscvec

2088:     PetscScalar, pointer :: xx_v(:)
2089:     ....
2090:     call VecGetArrayReadF90(x,xx_v,ierr)
2091:     a = xx_v(3)
2092:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2093: .ve

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

2097:     Level: beginner

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

2101: M*/

2103: /*MC
2104:     VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2105:     VecGetArrayReadF90().

2107:     Synopsis:
2108:     VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2110:     Logically Collective on Vec

2112:     Input Parameters:
2113: +   x - vector
2114: -   xx_v - the Fortran90 pointer to the array

2116:     Output Parameter:
2117: .   ierr - error code

2119:     Example of Usage:
2120: .vb
2121: #include <petsc/finclude/petscvec.h>
2122:     use petscvec

2124:     PetscScalar, pointer :: xx_v(:)
2125:     ....
2126:     call VecGetArrayReadF90(x,xx_v,ierr)
2127:     a = xx_v(3)
2128:     call VecRestoreArrayReadF90(x,xx_v,ierr)
2129: .ve

2131:     Level: beginner

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

2135: M*/

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

2142:    Logically Collective

2144:    Input Parameter:
2145: +  x - the vector
2146: .  m - first dimension of two dimensional array
2147: .  n - second dimension of two dimensional array
2148: .  mstart - first index you will use in first coordinate direction (often 0)
2149: -  nstart - first index in the second coordinate direction (often 0)

2151:    Output Parameter:
2152: .  a - location to put pointer to the array

2154:    Level: developer

2156:   Notes:
2157:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2158:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2159:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2160:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

2166: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2167:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2168:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2169: @*/
2170: PetscErrorCode  VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2171: {
2173:   PetscInt       i,N;
2174:   PetscScalar    *aa;

2180:   VecGetLocalSize(x,&N);
2181:   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);
2182:   VecGetArray(x,&aa);

2184:   PetscMalloc1(m,a);
2185:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2186:   *a -= mstart;
2187:   return(0);
2188: }

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

2193:    Logically Collective

2195:    Input Parameters:
2196: +  x - the vector
2197: .  m - first dimension of two dimensional array
2198: .  n - second dimension of the two dimensional array
2199: .  mstart - first index you will use in first coordinate direction (often 0)
2200: .  nstart - first index in the second coordinate direction (often 0)
2201: -  a - location of pointer to array obtained from VecGetArray2d()

2203:    Level: developer

2205:    Notes:
2206:    For regular PETSc vectors this routine does not involve any copies. For
2207:    any special vectors that do not store local vector data in a contiguous
2208:    array, this routine will copy the data back into the underlying
2209:    vector data structure from the array obtained with VecGetArray().

2211:    This routine actually zeros out the a pointer.

2213: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2214:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2215:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2216: @*/
2217: PetscErrorCode  VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2218: {
2220:   void           *dummy;

2226:   dummy = (void*)(*a + mstart);
2227:   PetscFree(dummy);
2228:   VecRestoreArray(x,NULL);
2229:   return(0);
2230: }

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

2237:    Logically Collective

2239:    Input Parameter:
2240: +  x - the vector
2241: .  m - first dimension of two dimensional array
2242: -  mstart - first index you will use in first coordinate direction (often 0)

2244:    Output Parameter:
2245: .  a - location to put pointer to the array

2247:    Level: developer

2249:   Notes:
2250:    For a vector obtained from DMCreateLocalVector() mstart are likely
2251:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2252:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2256: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2257:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2258:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2259: @*/
2260: PetscErrorCode  VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2261: {
2263:   PetscInt       N;

2269:   VecGetLocalSize(x,&N);
2270:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2271:   VecGetArray(x,a);
2272:   *a  -= mstart;
2273:   return(0);
2274: }

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

2279:    Logically Collective

2281:    Input Parameters:
2282: +  x - the vector
2283: .  m - first dimension of two dimensional array
2284: .  mstart - first index you will use in first coordinate direction (often 0)
2285: -  a - location of pointer to array obtained from VecGetArray21()

2287:    Level: developer

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

2295:    This routine actually zeros out the a pointer.

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

2299: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2300:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2301:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2302: @*/
2303: PetscErrorCode  VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2304: {

2310:   VecRestoreArray(x,NULL);
2311:   return(0);
2312: }


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

2320:    Logically Collective

2322:    Input Parameter:
2323: +  x - the vector
2324: .  m - first dimension of three dimensional array
2325: .  n - second dimension of three dimensional array
2326: .  p - third dimension of three dimensional array
2327: .  mstart - first index you will use in first coordinate direction (often 0)
2328: .  nstart - first index in the second coordinate direction (often 0)
2329: -  pstart - first index in the third coordinate direction (often 0)

2331:    Output Parameter:
2332: .  a - location to put pointer to the array

2334:    Level: developer

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

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

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

2346: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2347:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2348:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2349: @*/
2350: PetscErrorCode  VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2351: {
2353:   PetscInt       i,N,j;
2354:   PetscScalar    *aa,**b;

2360:   VecGetLocalSize(x,&N);
2361:   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);
2362:   VecGetArray(x,&aa);

2364:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
2365:   b    = (PetscScalar**)((*a) + m);
2366:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2367:   for (i=0; i<m; i++)
2368:     for (j=0; j<n; j++)
2369:       b[i*n+j] = aa + i*n*p + j*p - pstart;

2371:   *a -= mstart;
2372:   return(0);
2373: }

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

2378:    Logically Collective

2380:    Input Parameters:
2381: +  x - the vector
2382: .  m - first dimension of three dimensional array
2383: .  n - second dimension of the three dimensional array
2384: .  p - third dimension of the three dimensional array
2385: .  mstart - first index you will use in first coordinate direction (often 0)
2386: .  nstart - first index in the second coordinate direction (often 0)
2387: .  pstart - first index in the third coordinate direction (often 0)
2388: -  a - location of pointer to array obtained from VecGetArray3d()

2390:    Level: developer

2392:    Notes:
2393:    For regular PETSc vectors this routine does not involve any copies. For
2394:    any special vectors that do not store local vector data in a contiguous
2395:    array, this routine will copy the data back into the underlying
2396:    vector data structure from the array obtained with VecGetArray().

2398:    This routine actually zeros out the a pointer.

2400: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2401:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2402:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2403: @*/
2404: PetscErrorCode  VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2405: {
2407:   void           *dummy;

2413:   dummy = (void*)(*a + mstart);
2414:   PetscFree(dummy);
2415:   VecRestoreArray(x,NULL);
2416:   return(0);
2417: }

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

2424:    Logically Collective

2426:    Input Parameter:
2427: +  x - the vector
2428: .  m - first dimension of four dimensional array
2429: .  n - second dimension of four dimensional array
2430: .  p - third dimension of four dimensional array
2431: .  q - fourth dimension of four dimensional array
2432: .  mstart - first index you will use in first coordinate direction (often 0)
2433: .  nstart - first index in the second coordinate direction (often 0)
2434: .  pstart - first index in the third coordinate direction (often 0)
2435: -  qstart - first index in the fourth coordinate direction (often 0)

2437:    Output Parameter:
2438: .  a - location to put pointer to the array

2440:    Level: beginner

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

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

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

2452: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2453:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2454:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2455: @*/
2456: PetscErrorCode  VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2457: {
2459:   PetscInt       i,N,j,k;
2460:   PetscScalar    *aa,***b,**c;

2466:   VecGetLocalSize(x,&N);
2467:   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);
2468:   VecGetArray(x,&aa);

2470:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
2471:   b    = (PetscScalar***)((*a) + m);
2472:   c    = (PetscScalar**)(b + m*n);
2473:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2474:   for (i=0; i<m; i++)
2475:     for (j=0; j<n; j++)
2476:       b[i*n+j] = c + i*n*p + j*p - pstart;
2477:   for (i=0; i<m; i++)
2478:     for (j=0; j<n; j++)
2479:       for (k=0; k<p; k++)
2480:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
2481:   *a -= mstart;
2482:   return(0);
2483: }

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

2488:    Logically Collective

2490:    Input Parameters:
2491: +  x - the vector
2492: .  m - first dimension of four dimensional array
2493: .  n - second dimension of the four dimensional array
2494: .  p - third dimension of the four dimensional array
2495: .  q - fourth dimension of the four dimensional array
2496: .  mstart - first index you will use in first coordinate direction (often 0)
2497: .  nstart - first index in the second coordinate direction (often 0)
2498: .  pstart - first index in the third coordinate direction (often 0)
2499: .  qstart - first index in the fourth coordinate direction (often 0)
2500: -  a - location of pointer to array obtained from VecGetArray4d()

2502:    Level: beginner

2504:    Notes:
2505:    For regular PETSc vectors this routine does not involve any copies. For
2506:    any special vectors that do not store local vector data in a contiguous
2507:    array, this routine will copy the data back into the underlying
2508:    vector data structure from the array obtained with VecGetArray().

2510:    This routine actually zeros out the a pointer.

2512: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2513:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2514:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2515: @*/
2516: PetscErrorCode  VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2517: {
2519:   void           *dummy;

2525:   dummy = (void*)(*a + mstart);
2526:   PetscFree(dummy);
2527:   VecRestoreArray(x,NULL);
2528:   return(0);
2529: }

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

2536:    Logically Collective

2538:    Input Parameter:
2539: +  x - the vector
2540: .  m - first dimension of two dimensional array
2541: .  n - second dimension of two dimensional array
2542: .  mstart - first index you will use in first coordinate direction (often 0)
2543: -  nstart - first index in the second coordinate direction (often 0)

2545:    Output Parameter:
2546: .  a - location to put pointer to the array

2548:    Level: developer

2550:   Notes:
2551:    For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2552:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2553:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2554:    the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().

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

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

2560: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2561:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2562:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2563: @*/
2564: PetscErrorCode  VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2565: {
2566:   PetscErrorCode    ierr;
2567:   PetscInt          i,N;
2568:   const PetscScalar *aa;

2574:   VecGetLocalSize(x,&N);
2575:   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);
2576:   VecGetArrayRead(x,&aa);

2578:   PetscMalloc1(m,a);
2579:   for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
2580:   *a -= mstart;
2581:   return(0);
2582: }

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

2587:    Logically Collective

2589:    Input Parameters:
2590: +  x - the vector
2591: .  m - first dimension of two dimensional array
2592: .  n - second dimension of the two dimensional array
2593: .  mstart - first index you will use in first coordinate direction (often 0)
2594: .  nstart - first index in the second coordinate direction (often 0)
2595: -  a - location of pointer to array obtained from VecGetArray2d()

2597:    Level: developer

2599:    Notes:
2600:    For regular PETSc vectors this routine does not involve any copies. For
2601:    any special vectors that do not store local vector data in a contiguous
2602:    array, this routine will copy the data back into the underlying
2603:    vector data structure from the array obtained with VecGetArray().

2605:    This routine actually zeros out the a pointer.

2607: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2608:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2609:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2610: @*/
2611: PetscErrorCode  VecRestoreArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2612: {
2614:   void           *dummy;

2620:   dummy = (void*)(*a + mstart);
2621:   PetscFree(dummy);
2622:   VecRestoreArrayRead(x,NULL);
2623:   return(0);
2624: }

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

2631:    Logically Collective

2633:    Input Parameter:
2634: +  x - the vector
2635: .  m - first dimension of two dimensional array
2636: -  mstart - first index you will use in first coordinate direction (often 0)

2638:    Output Parameter:
2639: .  a - location to put pointer to the array

2641:    Level: developer

2643:   Notes:
2644:    For a vector obtained from DMCreateLocalVector() mstart are likely
2645:    obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2646:    DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().

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

2650: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2651:           VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2652:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2653: @*/
2654: PetscErrorCode  VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2655: {
2657:   PetscInt       N;

2663:   VecGetLocalSize(x,&N);
2664:   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2665:   VecGetArrayRead(x,(const PetscScalar**)a);
2666:   *a  -= mstart;
2667:   return(0);
2668: }

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

2673:    Logically Collective

2675:    Input Parameters:
2676: +  x - the vector
2677: .  m - first dimension of two dimensional array
2678: .  mstart - first index you will use in first coordinate direction (often 0)
2679: -  a - location of pointer to array obtained from VecGetArray21()

2681:    Level: developer

2683:    Notes:
2684:    For regular PETSc vectors this routine does not involve any copies. For
2685:    any special vectors that do not store local vector data in a contiguous
2686:    array, this routine will copy the data back into the underlying
2687:    vector data structure from the array obtained with VecGetArray1dRead().

2689:    This routine actually zeros out the a pointer.

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

2693: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2694:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2695:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2696: @*/
2697: PetscErrorCode  VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2698: {

2704:   VecRestoreArrayRead(x,NULL);
2705:   return(0);
2706: }


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

2714:    Logically Collective

2716:    Input Parameter:
2717: +  x - the vector
2718: .  m - first dimension of three dimensional array
2719: .  n - second dimension of three dimensional array
2720: .  p - third dimension of three dimensional array
2721: .  mstart - first index you will use in first coordinate direction (often 0)
2722: .  nstart - first index in the second coordinate direction (often 0)
2723: -  pstart - first index in the third coordinate direction (often 0)

2725:    Output Parameter:
2726: .  a - location to put pointer to the array

2728:    Level: developer

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

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

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

2740: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2741:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2742:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2743: @*/
2744: PetscErrorCode  VecGetArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2745: {
2746:   PetscErrorCode    ierr;
2747:   PetscInt          i,N,j;
2748:   const PetscScalar *aa;
2749:   PetscScalar       **b;

2755:   VecGetLocalSize(x,&N);
2756:   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);
2757:   VecGetArrayRead(x,&aa);

2759:   PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
2760:   b    = (PetscScalar**)((*a) + m);
2761:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2762:   for (i=0; i<m; i++)
2763:     for (j=0; j<n; j++)
2764:       b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;

2766:   *a -= mstart;
2767:   return(0);
2768: }

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

2773:    Logically Collective

2775:    Input Parameters:
2776: +  x - the vector
2777: .  m - first dimension of three dimensional array
2778: .  n - second dimension of the three dimensional array
2779: .  p - third dimension of the three dimensional array
2780: .  mstart - first index you will use in first coordinate direction (often 0)
2781: .  nstart - first index in the second coordinate direction (often 0)
2782: .  pstart - first index in the third coordinate direction (often 0)
2783: -  a - location of pointer to array obtained from VecGetArray3dRead()

2785:    Level: developer

2787:    Notes:
2788:    For regular PETSc vectors this routine does not involve any copies. For
2789:    any special vectors that do not store local vector data in a contiguous
2790:    array, this routine will copy the data back into the underlying
2791:    vector data structure from the array obtained with VecGetArray().

2793:    This routine actually zeros out the a pointer.

2795: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2796:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2797:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2798: @*/
2799: PetscErrorCode  VecRestoreArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2800: {
2802:   void           *dummy;

2808:   dummy = (void*)(*a + mstart);
2809:   PetscFree(dummy);
2810:   VecRestoreArrayRead(x,NULL);
2811:   return(0);
2812: }

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

2819:    Logically Collective

2821:    Input Parameter:
2822: +  x - the vector
2823: .  m - first dimension of four dimensional array
2824: .  n - second dimension of four dimensional array
2825: .  p - third dimension of four dimensional array
2826: .  q - fourth dimension of four dimensional array
2827: .  mstart - first index you will use in first coordinate direction (often 0)
2828: .  nstart - first index in the second coordinate direction (often 0)
2829: .  pstart - first index in the third coordinate direction (often 0)
2830: -  qstart - first index in the fourth coordinate direction (often 0)

2832:    Output Parameter:
2833: .  a - location to put pointer to the array

2835:    Level: beginner

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

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

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

2847: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2848:           VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2849:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2850: @*/
2851: PetscErrorCode  VecGetArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2852: {
2853:   PetscErrorCode    ierr;
2854:   PetscInt          i,N,j,k;
2855:   const PetscScalar *aa;
2856:   PetscScalar       ***b,**c;

2862:   VecGetLocalSize(x,&N);
2863:   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);
2864:   VecGetArrayRead(x,&aa);

2866:   PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
2867:   b    = (PetscScalar***)((*a) + m);
2868:   c    = (PetscScalar**)(b + m*n);
2869:   for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2870:   for (i=0; i<m; i++)
2871:     for (j=0; j<n; j++)
2872:       b[i*n+j] = c + i*n*p + j*p - pstart;
2873:   for (i=0; i<m; i++)
2874:     for (j=0; j<n; j++)
2875:       for (k=0; k<p; k++)
2876:         c[i*n*p+j*p+k] = (PetscScalar*) aa + i*n*p*q + j*p*q + k*q - qstart;
2877:   *a -= mstart;
2878:   return(0);
2879: }

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

2884:    Logically Collective

2886:    Input Parameters:
2887: +  x - the vector
2888: .  m - first dimension of four dimensional array
2889: .  n - second dimension of the four dimensional array
2890: .  p - third dimension of the four dimensional array
2891: .  q - fourth dimension of the four dimensional array
2892: .  mstart - first index you will use in first coordinate direction (often 0)
2893: .  nstart - first index in the second coordinate direction (often 0)
2894: .  pstart - first index in the third coordinate direction (often 0)
2895: .  qstart - first index in the fourth coordinate direction (often 0)
2896: -  a - location of pointer to array obtained from VecGetArray4dRead()

2898:    Level: beginner

2900:    Notes:
2901:    For regular PETSc vectors this routine does not involve any copies. For
2902:    any special vectors that do not store local vector data in a contiguous
2903:    array, this routine will copy the data back into the underlying
2904:    vector data structure from the array obtained with VecGetArray().

2906:    This routine actually zeros out the a pointer.

2908: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2909:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2910:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2911: @*/
2912: PetscErrorCode  VecRestoreArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2913: {
2915:   void           *dummy;

2921:   dummy = (void*)(*a + mstart);
2922:   PetscFree(dummy);
2923:   VecRestoreArrayRead(x,NULL);
2924:   return(0);
2925: }

2927: #if defined(PETSC_USE_DEBUG)

2929: /*@
2930:    VecLockGet  - Gets the current lock status of a vector

2932:    Logically Collective on Vec

2934:    Input Parameter:
2935: .  x - the vector

2937:    Output Parameter:
2938: .  state - greater than zero indicates the vector is locked for read; less then zero indicates the vector is
2939:            locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.

2941:    Level: beginner

2943:    Concepts: vector^accessing local values

2945: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop()
2946: @*/
2947: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
2948: {
2951:   *state = x->lock;
2952:   return(0);
2953: }

2955: /*@
2956:    VecLockReadPush  - Pushes a read-only lock on a vector to prevent it from writing

2958:    Logically Collective on Vec

2960:    Input Parameter:
2961: .  x - the vector

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

2966:     The call can be nested, i.e., called multiple times on the same vector, but each VecLockReadPush(x) has to have one matching
2967:     VecLockReadPop(x), which removes the latest read-only lock.

2969:    Level: beginner

2971:    Concepts: vector^accessing local values

2973: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPop(), VecLockGet()
2974: @*/
2975: PetscErrorCode VecLockReadPush(Vec x)
2976: {
2979:   if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector is already locked for exclusive write access but you want to read it");
2980:   x->lock++;
2981:   return(0);
2982: }

2984: /*@
2985:    VecLockReadPop  - Pops a read-only lock from a vector

2987:    Logically Collective on Vec

2989:    Input Parameter:
2990: .  x - the vector

2992:    Level: beginner

2994:    Concepts: vector^accessing local values

2996: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockGet()
2997: @*/
2998: PetscErrorCode VecLockReadPop(Vec x)
2999: {
3002:   x->lock--;
3003:   if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector has been unlocked from read-only access too many times");
3004:   return(0);
3005: }

3007: /*@C
3008:    VecLockWriteSet_Private  - Lock or unlock a vector for exclusive read/write access

3010:    Logically Collective on Vec

3012:    Input Parameter:
3013: +  x   - the vector
3014: -  flg - PETSC_TRUE to lock the vector for writing; PETSC_FALSE to unlock it.

3016:    Notes:
3017:     The function is usefull in split-phase computations, which usually have a begin phase and an end phase.
3018:     One can call VecLockWriteSet_Private(x,PETSC_TRUE) in the begin phase to lock a vector for exclusive
3019:     access, and call VecLockWriteSet_Private(x,PETSC_FALSE) in the end phase to unlock the vector from exclusive
3020:     access. In this way, one is ensured no other operations can access the vector in between. The code may like


3023:        VecGetArray(x,&xdata); // begin phase
3024:        VecLockWriteSet_Private(v,PETSC_TRUE);

3026:        Other operations, which can not acceess x anymore (they can access xdata, of course)

3028:        VecRestoreArray(x,&vdata); // end phase
3029:        VecLockWriteSet_Private(v,PETSC_FALSE);

3031:     The call can not be nested on the same vector, in other words, one can not call VecLockWriteSet_Private(x,PETSC_TRUE)
3032:     again before calling VecLockWriteSet_Private(v,PETSC_FALSE).

3034:    Level: beginner

3036:    Concepts: vector^accessing local values

3038: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop(), VecLockGet()
3039: @*/
3040: PetscErrorCode VecLockWriteSet_Private(Vec x,PetscBool flg)
3041: {
3044:   if (flg) {
3045:     if (x->lock > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector is already locked for read-only access but you want to write it");
3046:     else if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector is already locked for exclusive write access but you want to write it");
3047:     else x->lock = -1;
3048:   } else {
3049:     if (x->lock != -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector is not locked for exclusive write access but you want to unlock it from that");
3050:     x->lock = 0;
3051:   }
3052:   return(0);
3053: }

3055: /*@
3056:    VecLockPush  - Pushes a read-only lock on a vector to prevent it from writing

3058:    Level: deprecated

3060:    Concepts: vector^accessing local values

3062: .seealso: VecLockReadPush()
3063: @*/
3064: PetscErrorCode VecLockPush(Vec x)
3065: {
3068:   VecLockReadPush(x);
3069:   return(0);
3070: }

3072: /*@
3073:    VecLockPop  - Pops a read-only lock from a vector

3075:    Level: deprecated

3077:    Concepts: vector^accessing local values

3079: .seealso: VecLockReadPop()
3080: @*/
3081: PetscErrorCode VecLockPop(Vec x)
3082: {
3085:   VecLockReadPop(x);
3086:   return(0);
3087: }

3089: #endif