Actual source code: rvector.c
petsc-master 2020-11-28
1: /*
2: Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
3: These are the vector functions the user calls.
4: */
5: #include "petsc/private/sfimpl.h"
6: #include "petscsystypes.h"
7: #include <petsc/private/vecimpl.h>
8: #if defined(PETSC_HAVE_CUDA)
9: #include <../src/vec/vec/impls/dvecimpl.h>
10: #include <petsc/private/cudavecimpl.h>
11: #endif
12: static PetscInt VecGetSubVectorSavedStateId = -1;
14: PETSC_EXTERN PetscErrorCode VecValidValues(Vec vec,PetscInt argnum,PetscBool begin)
15: {
16: #if defined(PETSC_USE_DEBUG)
17: PetscErrorCode ierr;
18: PetscInt n,i;
19: const PetscScalar *x;
22: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL)
23: if ((vec->petscnative || vec->ops->getarray) && (vec->offloadmask & PETSC_OFFLOAD_CPU)) {
24: #else
25: if (vec->petscnative || vec->ops->getarray) {
26: #endif
27: VecGetLocalSize(vec,&n);
28: VecGetArrayRead(vec,&x);
29: for (i=0; i<n; i++) {
30: if (begin) {
31: 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);
32: } else {
33: 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);
34: }
35: }
36: VecRestoreArrayRead(vec,&x);
37: }
38: #else
40: #endif
41: return(0);
42: }
44: /*@
45: VecMaxPointwiseDivide - Computes the maximum of the componentwise division max = max_i abs(x_i/y_i).
47: Logically Collective on Vec
49: Input Parameters:
50: . x, y - the vectors
52: Output Parameter:
53: . max - the result
55: Level: advanced
57: Notes:
58: x and y may be the same vector
59: if a particular y_i is zero, it is treated as 1 in the above formula
61: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs()
62: @*/
63: PetscErrorCode VecMaxPointwiseDivide(Vec x,Vec y,PetscReal *max)
64: {
74: VecCheckSameSize(x,1,y,2);
75: (*x->ops->maxpointwisedivide)(x,y,max);
76: return(0);
77: }
79: /*@
80: VecDot - Computes the vector dot product.
82: Collective on Vec
84: Input Parameters:
85: . x, y - the vectors
87: Output Parameter:
88: . val - the dot product
90: Performance Issues:
91: $ per-processor memory bandwidth
92: $ interprocessor latency
93: $ work load inbalance that causes certain processes to arrive much earlier than others
95: Notes for Users of Complex Numbers:
96: For complex vectors, VecDot() computes
97: $ val = (x,y) = y^H x,
98: where y^H denotes the conjugate transpose of y. Note that this corresponds to the usual "mathematicians" complex
99: inner product where the SECOND argument gets the complex conjugate. Since the BLASdot() complex conjugates the first
100: first argument we call the BLASdot() with the arguments reversed.
102: Use VecTDot() for the indefinite form
103: $ val = (x,y) = y^T x,
104: where y^T denotes the transpose of y.
106: Level: intermediate
109: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDotRealPart()
110: @*/
111: PetscErrorCode VecDot(Vec x,Vec y,PetscScalar *val)
112: {
122: VecCheckSameSize(x,1,y,2);
124: PetscLogEventBegin(VEC_Dot,x,y,0,0);
125: (*x->ops->dot)(x,y,val);
126: PetscLogEventEnd(VEC_Dot,x,y,0,0);
127: return(0);
128: }
130: /*@
131: VecDotRealPart - Computes the real part of the vector dot product.
133: Collective on Vec
135: Input Parameters:
136: . x, y - the vectors
138: Output Parameter:
139: . val - the real part of the dot product;
141: Performance Issues:
142: $ per-processor memory bandwidth
143: $ interprocessor latency
144: $ work load inbalance that causes certain processes to arrive much earlier than others
146: Notes for Users of Complex Numbers:
147: See VecDot() for more details on the definition of the dot product for complex numbers
149: For real numbers this returns the same value as VecDot()
151: For complex numbers in C^n (that is a vector of n components with a complex number for each component) this is equal to the usual real dot product on the
152: the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)
154: Developer Note: This is not currently optimized to compute only the real part of the dot product.
156: Level: intermediate
159: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDot(), VecDotNorm2()
160: @*/
161: PetscErrorCode VecDotRealPart(Vec x,Vec y,PetscReal *val)
162: {
164: PetscScalar fdot;
167: VecDot(x,y,&fdot);
168: *val = PetscRealPart(fdot);
169: return(0);
170: }
172: /*@
173: VecNorm - Computes the vector norm.
175: Collective on Vec
177: Input Parameters:
178: + x - the vector
179: - type - one of NORM_1, NORM_2, NORM_INFINITY. Also available
180: NORM_1_AND_2, which computes both norms and stores them
181: in a two element array.
183: Output Parameter:
184: . val - the norm
186: Notes:
187: $ NORM_1 denotes sum_i |x_i|
188: $ NORM_2 denotes sqrt(sum_i |x_i|^2)
189: $ NORM_INFINITY denotes max_i |x_i|
191: For complex numbers NORM_1 will return the traditional 1 norm of the 2 norm of the complex numbers; that is the 1
192: norm of the absolute values of the complex entries. In PETSc 3.6 and earlier releases it returned the 1 norm of
193: the 1 norm of the complex entries (what is returned by the BLAS routine asum()). Both are valid norms but most
194: people expect the former.
196: Level: intermediate
198: Performance Issues:
199: $ per-processor memory bandwidth
200: $ interprocessor latency
201: $ work load inbalance that causes certain processes to arrive much earlier than others
204: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNormAvailable(),
205: VecNormBegin(), VecNormEnd()
207: @*/
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);
229: if (type!=NORM_1_AND_2) {
230: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[type],*val);
231: }
232: return(0);
233: }
235: /*@
236: VecNormAvailable - Returns the vector norm if it is already known.
238: Not Collective
240: Input Parameters:
241: + x - the vector
242: - type - one of NORM_1, NORM_2, NORM_INFINITY. Also available
243: NORM_1_AND_2, which computes both norms and stores them
244: in a two element array.
246: Output Parameter:
247: + available - PETSC_TRUE if the val returned is valid
248: - val - the norm
250: Notes:
251: $ NORM_1 denotes sum_i |x_i|
252: $ NORM_2 denotes sqrt(sum_i (x_i)^2)
253: $ NORM_INFINITY denotes max_i |x_i|
255: Level: intermediate
257: Performance Issues:
258: $ per-processor memory bandwidth
259: $ interprocessor latency
260: $ work load inbalance that causes certain processes to arrive much earlier than others
262: Compile Option:
263: PETSC_HAVE_SLOW_BLAS_NORM2 will cause a C (loop unrolled) version of the norm to be used, rather
264: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
265: (as opposed to vendor provided) because the FORTRAN BLAS NRM2() routine is very slow.
268: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNorm()
269: VecNormBegin(), VecNormEnd()
271: @*/
272: PetscErrorCode VecNormAvailable(Vec x,NormType type,PetscBool *available,PetscReal *val)
273: {
281: *available = PETSC_FALSE;
282: if (type!=NORM_1_AND_2) {
283: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,*available);
284: }
285: return(0);
286: }
288: /*@
289: VecNormalize - Normalizes a vector by 2-norm.
291: Collective on Vec
293: Input Parameters:
294: + x - the vector
296: Output Parameter:
297: . x - the normalized vector
298: - val - the vector norm before normalization
300: Level: intermediate
303: @*/
304: PetscErrorCode VecNormalize(Vec x,PetscReal *val)
305: {
307: PetscReal norm;
312: PetscLogEventBegin(VEC_Normalize,x,0,0,0);
313: VecNorm(x,NORM_2,&norm);
314: if (norm == 0.0) {
315: PetscInfo(x,"Vector of zero norm can not be normalized; Returning only the zero norm\n");
316: } else if (norm != 1.0) {
317: PetscScalar tmp = 1.0/norm;
318: VecScale(x,tmp);
319: }
320: if (val) *val = norm;
321: PetscLogEventEnd(VEC_Normalize,x,0,0,0);
322: return(0);
323: }
325: /*@C
326: VecMax - Determines the vector component with maximum real part and its location.
328: Collective on Vec
330: Input Parameter:
331: . x - the vector
333: Output Parameters:
334: + p - the location of val (pass NULL if you don't want this)
335: - val - the maximum component
337: Notes:
338: Returns the value PETSC_MIN_REAL and p = -1 if the vector is of length 0.
340: Returns the smallest index with the maximum value
341: Level: intermediate
344: .seealso: VecNorm(), VecMin()
345: @*/
346: PetscErrorCode VecMax(Vec x,PetscInt *p,PetscReal *val)
347: {
354: PetscLogEventBegin(VEC_Max,x,0,0,0);
355: (*x->ops->max)(x,p,val);
356: PetscLogEventEnd(VEC_Max,x,0,0,0);
357: return(0);
358: }
360: /*@C
361: VecMin - Determines the vector component with minimum real part and its location.
363: Collective on Vec
365: Input Parameters:
366: . x - the vector
368: Output Parameter:
369: + p - the location of val (pass NULL if you don't want this location)
370: - val - the minimum component
372: Level: intermediate
374: Notes:
375: Returns the value PETSC_MAX_REAL and p = -1 if the vector is of length 0.
377: This returns the smallest index with the minumum value
380: .seealso: VecMax()
381: @*/
382: PetscErrorCode VecMin(Vec x,PetscInt *p,PetscReal *val)
383: {
390: PetscLogEventBegin(VEC_Min,x,0,0,0);
391: (*x->ops->min)(x,p,val);
392: PetscLogEventEnd(VEC_Min,x,0,0,0);
393: return(0);
394: }
396: /*@
397: VecTDot - Computes an indefinite vector dot product. That is, this
398: routine does NOT use the complex conjugate.
400: Collective on Vec
402: Input Parameters:
403: . x, y - the vectors
405: Output Parameter:
406: . val - the dot product
408: Notes for Users of Complex Numbers:
409: For complex vectors, VecTDot() computes the indefinite form
410: $ val = (x,y) = y^T x,
411: where y^T denotes the transpose of y.
413: Use VecDot() for the inner product
414: $ val = (x,y) = y^H x,
415: where y^H denotes the conjugate transpose of y.
417: Level: intermediate
419: .seealso: VecDot(), VecMTDot()
420: @*/
421: PetscErrorCode VecTDot(Vec x,Vec y,PetscScalar *val)
422: {
432: VecCheckSameSize(x,1,y,2);
434: PetscLogEventBegin(VEC_TDot,x,y,0,0);
435: (*x->ops->tdot)(x,y,val);
436: PetscLogEventEnd(VEC_TDot,x,y,0,0);
437: return(0);
438: }
440: /*@
441: VecScale - Scales a vector.
443: Not collective on Vec
445: Input Parameters:
446: + x - the vector
447: - alpha - the scalar
449: Output Parameter:
450: . x - the scaled vector
452: Note:
453: For a vector with n components, VecScale() computes
454: $ x[i] = alpha * x[i], for i=1,...,n.
456: Level: intermediate
459: @*/
460: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
461: {
462: PetscReal norms[4] = {0.0,0.0,0.0, 0.0};
463: PetscBool flgs[4];
465: PetscInt i;
470: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
471: PetscLogEventBegin(VEC_Scale,x,0,0,0);
472: if (alpha != (PetscScalar)1.0) {
473: VecSetErrorIfLocked(x,1);
474: /* get current stashed norms */
475: for (i=0; i<4; i++) {
476: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
477: }
478: (*x->ops->scale)(x,alpha);
479: PetscObjectStateIncrease((PetscObject)x);
480: /* put the scaled stashed norms back into the Vec */
481: for (i=0; i<4; i++) {
482: if (flgs[i]) {
483: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],PetscAbsScalar(alpha)*norms[i]);
484: }
485: }
486: }
487: PetscLogEventEnd(VEC_Scale,x,0,0,0);
488: return(0);
489: }
491: /*@
492: VecSet - Sets all components of a vector to a single scalar value.
494: Logically Collective on Vec
496: Input Parameters:
497: + x - the vector
498: - alpha - the scalar
500: Output Parameter:
501: . x - the vector
503: Note:
504: For a vector of dimension n, VecSet() computes
505: $ x[i] = alpha, for i=1,...,n,
506: so that all vector entries then equal the identical
507: scalar value, alpha. Use the more general routine
508: VecSetValues() to set different vector entries.
510: You CANNOT call this after you have called VecSetValues() but before you call
511: VecAssemblyBegin/End().
513: Level: beginner
515: .seealso VecSetValues(), VecSetValuesBlocked(), VecSetRandom()
517: @*/
518: PetscErrorCode VecSet(Vec x,PetscScalar alpha)
519: {
520: PetscReal val;
526: 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()");
528: VecSetErrorIfLocked(x,1);
530: PetscLogEventBegin(VEC_Set,x,0,0,0);
531: (*x->ops->set)(x,alpha);
532: PetscLogEventEnd(VEC_Set,x,0,0,0);
533: PetscObjectStateIncrease((PetscObject)x);
535: /* norms can be simply set (if |alpha|*N not too large) */
536: val = PetscAbsScalar(alpha);
537: if (x->map->N == 0) {
538: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],0.0l);
539: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],0.0);
540: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],0.0);
541: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],0.0);
542: } else if (val > PETSC_MAX_REAL/x->map->N) {
543: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
544: } else {
545: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],x->map->N * val);
546: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
547: val = PetscSqrtReal((PetscReal)x->map->N) * val;
548: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],val);
549: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],val);
550: }
551: return(0);
552: }
555: /*@
556: VecAXPY - Computes y = alpha x + y.
558: Logically Collective on Vec
560: Input Parameters:
561: + alpha - the scalar
562: - x, y - the vectors
564: Output Parameter:
565: . y - output vector
567: Level: intermediate
569: Notes:
570: x and y MUST be different vectors
571: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
573: $ VecAXPY(y,alpha,x) y = alpha x + y
574: $ VecAYPX(y,beta,x) y = x + beta y
575: $ VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
576: $ VecWAXPY(w,alpha,x,y) w = alpha x + y
577: $ VecAXPBYPCZ(w,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
578: $ VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
581: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPBYPCZ(), VecAXPBY()
582: @*/
583: PetscErrorCode VecAXPY(Vec y,PetscScalar alpha,Vec x)
584: {
593: VecCheckSameSize(x,1,y,3);
594: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
596: if (alpha == (PetscScalar)0.0) return(0);
597: VecSetErrorIfLocked(y,1);
599: VecLockReadPush(x);
600: PetscLogEventBegin(VEC_AXPY,x,y,0,0);
601: (*y->ops->axpy)(y,alpha,x);
602: PetscLogEventEnd(VEC_AXPY,x,y,0,0);
603: VecLockReadPop(x);
604: PetscObjectStateIncrease((PetscObject)y);
605: return(0);
606: }
608: /*@
609: VecAXPBY - Computes y = alpha x + beta y.
611: Logically Collective on Vec
613: Input Parameters:
614: + alpha,beta - the scalars
615: - x, y - the vectors
617: Output Parameter:
618: . y - output vector
620: Level: intermediate
622: Notes:
623: x and y MUST be different vectors
624: The implementation is optimized for alpha and/or beta values of 0.0 and 1.0
627: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ()
628: @*/
629: PetscErrorCode VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
630: {
639: VecCheckSameSize(y,1,x,4);
640: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
643: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) return(0);
644: VecSetErrorIfLocked(y,1);
645: PetscLogEventBegin(VEC_AXPY,x,y,0,0);
646: (*y->ops->axpby)(y,alpha,beta,x);
647: PetscLogEventEnd(VEC_AXPY,x,y,0,0);
648: PetscObjectStateIncrease((PetscObject)y);
649: return(0);
650: }
652: /*@
653: VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z
655: Logically Collective on Vec
657: Input Parameters:
658: + alpha,beta, gamma - the scalars
659: - x, y, z - the vectors
661: Output Parameter:
662: . z - output vector
664: Level: intermediate
666: Notes:
667: x, y and z must be different vectors
668: The implementation is optimized for alpha of 1.0 and gamma of 1.0 or 0.0
671: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBY()
672: @*/
673: PetscErrorCode VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
674: {
686: VecCheckSameSize(x,1,y,5);
687: VecCheckSameSize(x,1,z,6);
688: if (x == y || x == z) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
689: if (y == z) SETERRQ(PetscObjectComm((PetscObject)y),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
693: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) return(0);
694: VecSetErrorIfLocked(z,1);
696: PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
697: (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
698: PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
699: PetscObjectStateIncrease((PetscObject)z);
700: return(0);
701: }
703: /*@
704: VecAYPX - Computes y = x + beta y.
706: Logically Collective on Vec
708: Input Parameters:
709: + beta - the scalar
710: - x, y - the vectors
712: Output Parameter:
713: . y - output vector
715: Level: intermediate
717: Notes:
718: x and y MUST be different vectors
719: The implementation is optimized for beta of -1.0, 0.0, and 1.0
722: .seealso: VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
723: @*/
724: PetscErrorCode VecAYPX(Vec y,PetscScalar beta,Vec x)
725: {
734: VecCheckSameSize(x,1,y,3);
735: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
737: VecSetErrorIfLocked(y,1);
739: PetscLogEventBegin(VEC_AYPX,x,y,0,0);
740: (*y->ops->aypx)(y,beta,x);
741: PetscLogEventEnd(VEC_AYPX,x,y,0,0);
742: PetscObjectStateIncrease((PetscObject)y);
743: return(0);
744: }
747: /*@
748: VecWAXPY - Computes w = alpha x + y.
750: Logically Collective on Vec
752: Input Parameters:
753: + alpha - the scalar
754: - x, y - the vectors
756: Output Parameter:
757: . w - the result
759: Level: intermediate
761: Notes:
762: w cannot be either x or y, but x and y can be the same
763: The implementation is optimzed for alpha of -1.0, 0.0, and 1.0
766: .seealso: VecAXPY(), VecAYPX(), VecAXPBY(), VecMAXPY(), VecAXPBYPCZ()
767: @*/
768: PetscErrorCode VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
769: {
781: VecCheckSameSize(x,3,y,4);
782: VecCheckSameSize(x,3,w,1);
783: if (w == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector y, suggest VecAXPY()");
784: if (w == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector x, suggest VecAYPX()");
786: VecSetErrorIfLocked(w,1);
788: PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
789: (*w->ops->waxpy)(w,alpha,x,y);
790: PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
791: PetscObjectStateIncrease((PetscObject)w);
792: return(0);
793: }
796: /*@C
797: VecSetValues - Inserts or adds values into certain locations of a vector.
799: Not Collective
801: Input Parameters:
802: + x - vector to insert in
803: . ni - number of elements to add
804: . ix - indices where to add
805: . y - array of values
806: - iora - either INSERT_VALUES or ADD_VALUES, where
807: ADD_VALUES adds values to any existing entries, and
808: INSERT_VALUES replaces existing entries with new values
810: Notes:
811: VecSetValues() sets x[ix[i]] = y[i], for i=0,...,ni-1.
813: Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
814: options cannot be mixed without intervening calls to the assembly
815: routines.
817: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
818: MUST be called after all calls to VecSetValues() have been completed.
820: VecSetValues() uses 0-based indices in Fortran as well as in C.
822: If you call VecSetOption(x, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE),
823: negative indices may be passed in ix. These rows are
824: simply ignored. This allows easily inserting element load matrices
825: with homogeneous Dirchlet boundary conditions that you don't want represented
826: in the vector.
828: Level: beginner
830: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
831: VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
832: @*/
833: PetscErrorCode VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
834: {
839: if (!ni) return(0);
844: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
845: (*x->ops->setvalues)(x,ni,ix,y,iora);
846: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
847: PetscObjectStateIncrease((PetscObject)x);
848: return(0);
849: }
851: /*@C
852: VecGetValues - Gets values from certain locations of a vector. Currently
853: can only get values on the same processor
855: Not Collective
857: Input Parameters:
858: + x - vector to get values from
859: . ni - number of elements to get
860: - ix - indices where to get them from (in global 1d numbering)
862: Output Parameter:
863: . y - array of values
865: Notes:
866: The user provides the allocated array y; it is NOT allocated in this routine
868: VecGetValues() gets y[i] = x[ix[i]], for i=0,...,ni-1.
870: VecAssemblyBegin() and VecAssemblyEnd() MUST be called before calling this
872: VecGetValues() uses 0-based indices in Fortran as well as in C.
874: If you call VecSetOption(x, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE),
875: negative indices may be passed in ix. These rows are
876: simply ignored.
878: Level: beginner
880: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues()
881: @*/
882: PetscErrorCode VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
883: {
888: if (!ni) return(0);
892: (*x->ops->getvalues)(x,ni,ix,y);
893: return(0);
894: }
896: /*@C
897: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
899: Not Collective
901: Input Parameters:
902: + x - vector to insert in
903: . ni - number of blocks to add
904: . ix - indices where to add in block count, rather than element count
905: . y - array of values
906: - iora - either INSERT_VALUES or ADD_VALUES, where
907: ADD_VALUES adds values to any existing entries, and
908: INSERT_VALUES replaces existing entries with new values
910: Notes:
911: VecSetValuesBlocked() sets x[bs*ix[i]+j] = y[bs*i+j],
912: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
914: Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
915: options cannot be mixed without intervening calls to the assembly
916: routines.
918: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
919: MUST be called after all calls to VecSetValuesBlocked() have been completed.
921: VecSetValuesBlocked() uses 0-based indices in Fortran as well as in C.
923: Negative indices may be passed in ix, these rows are
924: simply ignored. This allows easily inserting element load matrices
925: with homogeneous Dirchlet boundary conditions that you don't want represented
926: in the vector.
928: Level: intermediate
930: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
931: VecSetValues()
932: @*/
933: PetscErrorCode VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
934: {
939: if (!ni) return(0);
944: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
945: (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
946: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
947: PetscObjectStateIncrease((PetscObject)x);
948: return(0);
949: }
952: /*@C
953: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
954: using a local ordering of the nodes.
956: Not Collective
958: Input Parameters:
959: + x - vector to insert in
960: . ni - number of elements to add
961: . ix - indices where to add
962: . y - array of values
963: - iora - either INSERT_VALUES or ADD_VALUES, where
964: ADD_VALUES adds values to any existing entries, and
965: INSERT_VALUES replaces existing entries with new values
967: Level: intermediate
969: Notes:
970: VecSetValuesLocal() sets x[ix[i]] = y[i], for i=0,...,ni-1.
972: Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
973: options cannot be mixed without intervening calls to the assembly
974: routines.
976: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
977: MUST be called after all calls to VecSetValuesLocal() have been completed.
979: VecSetValuesLocal() uses 0-based indices in Fortran as well as in C.
981: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetLocalToGlobalMapping(),
982: VecSetValuesBlockedLocal()
983: @*/
984: PetscErrorCode VecSetValuesLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
985: {
987: PetscInt lixp[128],*lix = lixp;
991: if (!ni) return(0);
996: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
997: if (!x->ops->setvalueslocal) {
998: if (!x->map->mapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMapping()");
999: if (ni > 128) {
1000: PetscMalloc1(ni,&lix);
1001: }
1002: ISLocalToGlobalMappingApply(x->map->mapping,ni,(PetscInt*)ix,lix);
1003: (*x->ops->setvalues)(x,ni,lix,y,iora);
1004: if (ni > 128) {
1005: PetscFree(lix);
1006: }
1007: } else {
1008: (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
1009: }
1010: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1011: PetscObjectStateIncrease((PetscObject)x);
1012: return(0);
1013: }
1015: /*@
1016: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1017: using a local ordering of the nodes.
1019: Not Collective
1021: Input Parameters:
1022: + x - vector to insert in
1023: . ni - number of blocks to add
1024: . ix - indices where to add in block count, not element count
1025: . y - array of values
1026: - iora - either INSERT_VALUES or ADD_VALUES, where
1027: ADD_VALUES adds values to any existing entries, and
1028: INSERT_VALUES replaces existing entries with new values
1030: Level: intermediate
1032: Notes:
1033: VecSetValuesBlockedLocal() sets x[bs*ix[i]+j] = y[bs*i+j],
1034: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with VecSetBlockSize().
1036: Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1037: options cannot be mixed without intervening calls to the assembly
1038: routines.
1040: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
1041: MUST be called after all calls to VecSetValuesBlockedLocal() have been completed.
1043: VecSetValuesBlockedLocal() uses 0-based indices in Fortran as well as in C.
1046: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
1047: VecSetLocalToGlobalMapping()
1048: @*/
1049: PetscErrorCode VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1050: {
1052: PetscInt lixp[128],*lix = lixp;
1056: if (!ni) return(0);
1060: if (ni > 128) {
1061: PetscMalloc1(ni,&lix);
1062: }
1064: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1065: ISLocalToGlobalMappingApplyBlock(x->map->mapping,ni,(PetscInt*)ix,lix);
1066: (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1067: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1068: if (ni > 128) {
1069: PetscFree(lix);
1070: }
1071: PetscObjectStateIncrease((PetscObject)x);
1072: return(0);
1073: }
1075: /*@
1076: VecMTDot - Computes indefinite vector multiple dot products.
1077: That is, it does NOT use the complex conjugate.
1079: Collective on Vec
1081: Input Parameters:
1082: + x - one vector
1083: . nv - number of vectors
1084: - y - array of vectors. Note that vectors are pointers
1086: Output Parameter:
1087: . val - array of the dot products
1089: Notes for Users of Complex Numbers:
1090: For complex vectors, VecMTDot() computes the indefinite form
1091: $ val = (x,y) = y^T x,
1092: where y^T denotes the transpose of y.
1094: Use VecMDot() for the inner product
1095: $ val = (x,y) = y^H x,
1096: where y^H denotes the conjugate transpose of y.
1098: Level: intermediate
1101: .seealso: VecMDot(), VecTDot()
1102: @*/
1103: PetscErrorCode VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1104: {
1110: if (!nv) return(0);
1117: VecCheckSameSize(x,1,*y,3);
1119: PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1120: (*x->ops->mtdot)(x,nv,y,val);
1121: PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1122: return(0);
1123: }
1125: /*@
1126: VecMDot - Computes vector multiple dot products.
1128: Collective on Vec
1130: Input Parameters:
1131: + x - one vector
1132: . nv - number of vectors
1133: - y - array of vectors.
1135: Output Parameter:
1136: . val - array of the dot products (does not allocate the array)
1138: Notes for Users of Complex Numbers:
1139: For complex vectors, VecMDot() computes
1140: $ val = (x,y) = y^H x,
1141: where y^H denotes the conjugate transpose of y.
1143: Use VecMTDot() for the indefinite form
1144: $ val = (x,y) = y^T x,
1145: where y^T denotes the transpose of y.
1147: Level: intermediate
1150: .seealso: VecMTDot(), VecDot()
1151: @*/
1152: PetscErrorCode VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1153: {
1159: if (!nv) return(0);
1160: if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1167: VecCheckSameSize(x,1,*y,3);
1169: PetscLogEventBegin(VEC_MDot,x,*y,0,0);
1170: (*x->ops->mdot)(x,nv,y,val);
1171: PetscLogEventEnd(VEC_MDot,x,*y,0,0);
1172: return(0);
1173: }
1175: /*@
1176: VecMAXPY - Computes y = y + sum alpha[i] x[i]
1178: Logically Collective on Vec
1180: Input Parameters:
1181: + nv - number of scalars and x-vectors
1182: . alpha - array of scalars
1183: . y - one vector
1184: - x - array of vectors
1186: Level: intermediate
1188: Notes:
1189: y cannot be any of the x vectors
1191: .seealso: VecAYPX(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
1192: @*/
1193: PetscErrorCode VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1194: {
1196: PetscInt i;
1197: PetscBool nonzero;
1202: if (!nv) return(0);
1203: if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1210: VecCheckSameSize(y,1,*x,4);
1212: for (i=0, nonzero = PETSC_FALSE; i<nv && !nonzero; i++) nonzero = (PetscBool)(nonzero || alpha[i] != (PetscScalar)0.0);
1213: if (!nonzero) return(0);
1214: VecSetErrorIfLocked(y,1);
1215: PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1216: (*y->ops->maxpy)(y,nv,alpha,x);
1217: PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1218: PetscObjectStateIncrease((PetscObject)y);
1219: return(0);
1220: }
1222: /*@
1223: VecGetSubVector - Gets a vector representing part of another vector
1225: Collective on X and IS
1227: Input Arguments:
1228: + X - vector from which to extract a subvector
1229: - is - index set representing portion of X to extract
1231: Output Arguments:
1232: . Y - subvector corresponding to is
1234: Level: advanced
1236: Notes:
1237: The subvector Y should be returned with VecRestoreSubVector().
1238: X and is must be defined on the same communicator
1240: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1241: modifying the subvector. Other non-overlapping subvectors can still be obtained from X using this function.
1242: The resulting subvector inherits the block size from the IS if greater than one. Otherwise, the block size is guessed from the block size of the original vec.
1244: .seealso: MatCreateSubMatrix()
1245: @*/
1246: PetscErrorCode VecGetSubVector(Vec X,IS is,Vec *Y)
1247: {
1248: PetscErrorCode ierr;
1249: Vec Z;
1256: if (X->ops->getsubvector) {
1257: (*X->ops->getsubvector)(X,is,&Z);
1258: } else { /* Default implementation currently does no caching */
1259: PetscInt gstart,gend,start;
1260: PetscBool red[2] = { PETSC_TRUE, PETSC_TRUE };
1261: PetscInt n,N,ibs,vbs,bs = -1;
1263: ISGetLocalSize(is,&n);
1264: ISGetSize(is,&N);
1265: ISGetBlockSize(is,&ibs);
1266: VecGetBlockSize(X,&vbs);
1267: VecGetOwnershipRange(X,&gstart,&gend);
1268: ISContiguousLocal(is,gstart,gend,&start,&red[0]);
1269: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1270: if (ibs > 1) {
1271: MPIU_Allreduce(MPI_IN_PLACE,red,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1272: bs = ibs;
1273: } else {
1274: if (n%vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1275: MPIU_Allreduce(MPI_IN_PLACE,red,2,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1276: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1277: }
1278: if (red[0]) { /* We can do a no-copy implementation */
1279: const PetscScalar *x;
1280: PetscInt state = 0;
1281: PetscBool isstd;
1282: #if defined(PETSC_HAVE_CUDA)
1283: PetscBool iscuda;
1284: #endif
1286: PetscObjectTypeCompareAny((PetscObject)X,&isstd,VECSEQ,VECMPI,VECSTANDARD,"");
1287: #if defined(PETSC_HAVE_CUDA)
1288: PetscObjectTypeCompareAny((PetscObject)X,&iscuda,VECSEQCUDA,VECMPICUDA,"");
1289: if (iscuda) {
1290: const PetscScalar *x_d;
1291: PetscMPIInt size;
1292: PetscOffloadMask flg;
1294: VecCUDAGetArrays_Private(X,&x,&x_d,&flg);
1295: if (flg == PETSC_OFFLOAD_UNALLOCATED) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for PETSC_OFFLOAD_UNALLOCATED");
1296: if (n && !x && !x_d) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Missing vector data");
1297: if (x) x += start;
1298: if (x_d) x_d += start;
1299: MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1300: if (size == 1) {
1301: VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X),bs,n,x,x_d,&Z);
1302: } else {
1303: VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X),bs,n,N,x,x_d,&Z);
1304: }
1305: Z->offloadmask = flg;
1306: } else if (isstd) {
1307: #else
1308: if (isstd) { /* standard CPU: use CreateWithArray pattern */
1309: #endif
1310: PetscMPIInt size;
1312: MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1313: VecGetArrayRead(X,&x);
1314: if (x) x += start;
1315: if (size == 1) {
1316: VecCreateSeqWithArray(PetscObjectComm((PetscObject)X),bs,n,x,&Z);
1317: } else {
1318: VecCreateMPIWithArray(PetscObjectComm((PetscObject)X),bs,n,N,x,&Z);
1319: }
1320: VecRestoreArrayRead(X,&x);
1321: } else { /* default implementation: use place array */
1322: VecGetArrayRead(X,&x);
1323: VecCreate(PetscObjectComm((PetscObject)X),&Z);
1324: VecSetType(Z,((PetscObject)X)->type_name);
1325: VecSetSizes(Z,n,N);
1326: VecSetBlockSize(Z,bs);
1327: VecPlaceArray(Z,x ? x+start : NULL);
1328: VecRestoreArrayRead(X,&x);
1329: }
1331: /* this is relevant only in debug mode */
1332: VecLockGet(X,&state);
1333: if (state) {
1334: VecLockReadPush(Z);
1335: }
1336: Z->ops->placearray = NULL;
1337: Z->ops->replacearray = NULL;
1338: } else { /* Have to create a scatter and do a copy */
1339: VecScatter scatter;
1341: VecCreate(PetscObjectComm((PetscObject)is),&Z);
1342: VecSetSizes(Z,n,N);
1343: VecSetBlockSize(Z,bs);
1344: VecSetType(Z,((PetscObject)X)->type_name);
1345: VecScatterCreate(X,is,Z,NULL,&scatter);
1346: VecScatterBegin(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1347: VecScatterEnd(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1348: PetscObjectCompose((PetscObject)Z,"VecGetSubVector_Scatter",(PetscObject)scatter);
1349: VecScatterDestroy(&scatter);
1350: }
1351: }
1352: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1353: if (VecGetSubVectorSavedStateId < 0) {PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);}
1354: PetscObjectComposedDataSetInt((PetscObject)Z,VecGetSubVectorSavedStateId,1);
1355: *Y = Z;
1356: return(0);
1357: }
1359: /*@
1360: VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()
1362: Collective on IS
1364: Input Arguments:
1365: + X - vector from which subvector was obtained
1366: . is - index set representing the subset of X
1367: - Y - subvector being restored
1369: Level: advanced
1371: .seealso: VecGetSubVector()
1372: @*/
1373: PetscErrorCode VecRestoreSubVector(Vec X,IS is,Vec *Y)
1374: {
1383: if (X->ops->restoresubvector) {
1384: (*X->ops->restoresubvector)(X,is,Y);
1385: } else {
1386: PETSC_UNUSED PetscObjectState dummystate = 0;
1387: PetscBool valid;
1389: PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,dummystate,valid);
1390: if (!valid) {
1391: VecScatter scatter;
1392: PetscInt state;
1394: VecLockGet(X,&state);
1395: if (state != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vec X is locked for read-only or read/write access");
1397: PetscObjectQuery((PetscObject)*Y,"VecGetSubVector_Scatter",(PetscObject*)&scatter);
1398: if (scatter) {
1399: VecScatterBegin(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1400: VecScatterEnd(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1401: } else {
1402: #if defined(PETSC_HAVE_CUDA)
1403: PetscBool iscuda;
1405: PetscObjectTypeCompareAny((PetscObject)*Y,&iscuda,VECSEQCUDA,VECMPICUDA,"");
1406: if (iscuda) {
1407: PetscOffloadMask ymask = (*Y)->offloadmask;
1409: /* The offloadmask of X dictates where to move memory
1410: If X GPU data is valid, then move Y data on GPU if needed
1411: Otherwise, move back to the CPU */
1412: switch (X->offloadmask) {
1413: case PETSC_OFFLOAD_BOTH:
1414: if (ymask == PETSC_OFFLOAD_CPU) {
1415: VecCUDAResetArray(*Y);
1416: } else if (ymask == PETSC_OFFLOAD_GPU) {
1417: X->offloadmask = PETSC_OFFLOAD_GPU;
1418: }
1419: break;
1420: case PETSC_OFFLOAD_GPU:
1421: if (ymask == PETSC_OFFLOAD_CPU) {
1422: VecCUDAResetArray(*Y);
1423: }
1424: break;
1425: case PETSC_OFFLOAD_CPU:
1426: if (ymask == PETSC_OFFLOAD_GPU) {
1427: VecResetArray(*Y);
1428: }
1429: break;
1430: case PETSC_OFFLOAD_UNALLOCATED:
1431: case PETSC_OFFLOAD_VECKOKKOS:
1432: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"This should not happen");
1433: break;
1434: }
1435: } else {
1436: #else
1437: {
1438: #endif
1439: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1440: VecResetArray(*Y);
1441: }
1442: PetscObjectStateIncrease((PetscObject)X);
1443: }
1444: }
1445: VecDestroy(Y);
1446: }
1447: return(0);
1448: }
1450: /*@
1451: VecGetLocalVectorRead - Maps the local portion of a vector into a
1452: vector. You must call VecRestoreLocalVectorRead() when the local
1453: vector is no longer needed.
1455: Not collective.
1457: Input parameter:
1458: . v - The vector for which the local vector is desired.
1460: Output parameter:
1461: . w - Upon exit this contains the local vector.
1463: Level: beginner
1465: Notes:
1466: This function is similar to VecGetArrayRead() which maps the local
1467: portion into a raw pointer. VecGetLocalVectorRead() is usually
1468: almost as efficient as VecGetArrayRead() but in certain circumstances
1469: VecGetLocalVectorRead() can be much more efficient than
1470: VecGetArrayRead(). This is because the construction of a contiguous
1471: array representing the vector data required by VecGetArrayRead() can
1472: be an expensive operation for certain vector types. For example, for
1473: GPU vectors VecGetArrayRead() requires that the data between device
1474: and host is synchronized.
1476: Unlike VecGetLocalVector(), this routine is not collective and
1477: preserves cached information.
1479: .seealso: VecRestoreLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1480: @*/
1481: PetscErrorCode VecGetLocalVectorRead(Vec v,Vec w)
1482: {
1484: PetscScalar *a;
1489: VecCheckSameLocalSize(v,1,w,2);
1490: if (v->ops->getlocalvectorread) {
1491: (*v->ops->getlocalvectorread)(v,w);
1492: } else {
1493: VecGetArrayRead(v,(const PetscScalar**)&a);
1494: VecPlaceArray(w,a);
1495: }
1496: PetscObjectStateIncrease((PetscObject)w);
1497: VecLockReadPush(v);
1498: VecLockReadPush(w);
1499: return(0);
1500: }
1502: /*@
1503: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1504: previously mapped into a vector using VecGetLocalVectorRead().
1506: Not collective.
1508: Input parameter:
1509: + v - The local portion of this vector was previously mapped into w using VecGetLocalVectorRead().
1510: - w - The vector into which the local portion of v was mapped.
1512: Level: beginner
1514: .seealso: VecGetLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1515: @*/
1516: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec w)
1517: {
1519: PetscScalar *a;
1524: if (v->ops->restorelocalvectorread) {
1525: (*v->ops->restorelocalvectorread)(v,w);
1526: } else {
1527: VecGetArrayRead(w,(const PetscScalar**)&a);
1528: VecRestoreArrayRead(v,(const PetscScalar**)&a);
1529: VecResetArray(w);
1530: }
1531: VecLockReadPop(v);
1532: VecLockReadPop(w);
1533: PetscObjectStateIncrease((PetscObject)w);
1534: return(0);
1535: }
1537: /*@
1538: VecGetLocalVector - Maps the local portion of a vector into a
1539: vector.
1541: Collective on v, not collective on w.
1543: Input parameter:
1544: . v - The vector for which the local vector is desired.
1546: Output parameter:
1547: . w - Upon exit this contains the local vector.
1549: Level: beginner
1551: Notes:
1552: This function is similar to VecGetArray() which maps the local
1553: portion into a raw pointer. VecGetLocalVector() is usually about as
1554: efficient as VecGetArray() but in certain circumstances
1555: VecGetLocalVector() can be much more efficient than VecGetArray().
1556: This is because the construction of a contiguous array representing
1557: the vector data required by VecGetArray() can be an expensive
1558: operation for certain vector types. For example, for GPU vectors
1559: VecGetArray() requires that the data between device and host is
1560: synchronized.
1562: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1563: @*/
1564: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1565: {
1567: PetscScalar *a;
1572: VecCheckSameLocalSize(v,1,w,2);
1573: if (v->ops->getlocalvector) {
1574: (*v->ops->getlocalvector)(v,w);
1575: } else {
1576: VecGetArray(v,&a);
1577: VecPlaceArray(w,a);
1578: }
1579: PetscObjectStateIncrease((PetscObject)w);
1580: return(0);
1581: }
1583: /*@
1584: VecRestoreLocalVector - Unmaps the local portion of a vector
1585: previously mapped into a vector using VecGetLocalVector().
1587: Logically collective.
1589: Input parameter:
1590: + v - The local portion of this vector was previously mapped into w using VecGetLocalVector().
1591: - w - The vector into which the local portion of v was mapped.
1593: Level: beginner
1595: .seealso: VecGetLocalVector(), VecGetLocalVectorRead(), VecRestoreLocalVectorRead(), LocalVectorRead(), VecGetArrayRead(), VecGetArray()
1596: @*/
1597: PetscErrorCode VecRestoreLocalVector(Vec v,Vec w)
1598: {
1600: PetscScalar *a;
1605: if (v->ops->restorelocalvector) {
1606: (*v->ops->restorelocalvector)(v,w);
1607: } else {
1608: VecGetArray(w,&a);
1609: VecRestoreArray(v,&a);
1610: VecResetArray(w);
1611: }
1612: PetscObjectStateIncrease((PetscObject)w);
1613: PetscObjectStateIncrease((PetscObject)v);
1614: return(0);
1615: }
1617: /*@C
1618: VecGetArray - Returns a pointer to a contiguous array that contains this
1619: processor's portion of the vector data. For the standard PETSc
1620: vectors, VecGetArray() returns a pointer to the local data array and
1621: does not use any copies. If the underlying vector data is not stored
1622: in a contiguous array this routine will copy the data to a contiguous
1623: array and return a pointer to that. You MUST call VecRestoreArray()
1624: when you no longer need access to the array.
1626: Logically Collective on Vec
1628: Input Parameter:
1629: . x - the vector
1631: Output Parameter:
1632: . a - location to put pointer to the array
1634: Fortran Note:
1635: This routine is used differently from Fortran 77
1636: $ Vec x
1637: $ PetscScalar x_array(1)
1638: $ PetscOffset i_x
1639: $ PetscErrorCode ierr
1640: $ call VecGetArray(x,x_array,i_x,ierr)
1641: $
1642: $ Access first local entry in vector with
1643: $ value = x_array(i_x + 1)
1644: $
1645: $ ...... other code
1646: $ call VecRestoreArray(x,x_array,i_x,ierr)
1647: For Fortran 90 see VecGetArrayF90()
1649: See the Fortran chapter of the users manual and
1650: petsc/src/snes/tutorials/ex5f.F for details.
1652: Level: beginner
1654: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1655: VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1656: @*/
1657: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1658: {
1663: VecSetErrorIfLocked(x,1);
1664: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
1665: (*x->ops->getarray)(x,a);
1666: } else if (x->petscnative) { /* VECSTANDARD */
1667: *a = *((PetscScalar**)x->data);
1668: } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array for vector type \"%s\"",((PetscObject)x)->type_name);
1669: return(0);
1670: }
1672: /*@C
1673: VecRestoreArray - Restores a vector after VecGetArray() has been called.
1675: Logically Collective on Vec
1677: Input Parameters:
1678: + x - the vector
1679: - a - location of pointer to array obtained from VecGetArray()
1681: Level: beginner
1683: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1684: VecGetArrayPair(), VecRestoreArrayPair()
1685: @*/
1686: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1687: {
1692: if (x->ops->restorearray) { /* VECNEST, VECCUDA etc */
1693: (*x->ops->restorearray)(x,a);
1694: } else if (x->petscnative) { /* VECSTANDARD */
1695: /* nothing */
1696: } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot restore array for vector type \"%s\"",((PetscObject)x)->type_name);
1697: if (a) *a = NULL;
1698: PetscObjectStateIncrease((PetscObject)x);
1699: return(0);
1700: }
1701: /*@C
1702: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
1704: Not Collective
1706: Input Parameter:
1707: . x - the vector
1709: Output Parameter:
1710: . a - the array
1712: Level: beginner
1714: Notes:
1715: The array must be returned using a matching call to VecRestoreArrayRead().
1717: Unlike VecGetArray(), this routine is not collective and preserves cached information like vector norms.
1719: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
1720: implementations may require a copy, but must such implementations should cache the contiguous representation so that
1721: only one copy is performed when this routine is called multiple times in sequence.
1723: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1724: @*/
1725: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1726: {
1731: if (x->ops->getarray) { /* VECNEST, VECCUDA, VECKOKKOS etc */
1732: (*x->ops->getarray)(x,(PetscScalar**)a);
1733: } else if (x->petscnative) { /* VECSTANDARD */
1734: *a = *((PetscScalar**)x->data);
1735: } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array read for vector type \"%s\"",((PetscObject)x)->type_name);
1736: return(0);
1737: }
1739: /*@C
1740: VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()
1742: Not Collective
1744: Input Parameters:
1745: + vec - the vector
1746: - array - the array
1748: Level: beginner
1750: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1751: @*/
1752: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1753: {
1758: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
1759: /* nothing */
1760: } else if (x->ops->restorearrayread) { /* VECNEST */
1761: (*x->ops->restorearrayread)(x,a);
1762: } else { /* No one? */
1763: (*x->ops->restorearray)(x,(PetscScalar**)a);
1764: }
1765: if (a) *a = NULL;
1766: return(0);
1767: }
1769: /*@C
1770: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contains this
1771: processor's portion of the vector data. The values in this array are NOT valid, the routine calling this
1772: routine is responsible for putting values into the array; any values it does not set will be invalid
1774: Logically Collective on Vec
1776: Input Parameter:
1777: . x - the vector
1779: Output Parameter:
1780: . a - location to put pointer to the array
1782: Level: intermediate
1784: This is for vectors associate with GPUs, the vector is not copied up before giving access. If you need correct
1785: values in the array use VecGetArray()
1787: Concepts: vector^accessing local values
1789: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1790: VecGetArrayPair(), VecRestoreArrayPair(), VecGetArray(), VecRestoreArrayWrite()
1791: @*/
1792: PetscErrorCode VecGetArrayWrite(Vec x,PetscScalar **a)
1793: {
1798: VecSetErrorIfLocked(x,1);
1799: if (x->ops->getarraywrite) {
1800: (*x->ops->getarraywrite)(x,a);
1801: } else {
1802: VecGetArray(x,a);
1803: }
1804: return(0);
1805: }
1807: /*@C
1808: VecRestoreArrayWrite - Restores a vector after VecGetArrayWrite() has been called.
1810: Logically Collective on Vec
1812: Input Parameters:
1813: + x - the vector
1814: - a - location of pointer to array obtained from VecGetArray()
1816: Level: beginner
1818: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1819: VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite()
1820: @*/
1821: PetscErrorCode VecRestoreArrayWrite(Vec x,PetscScalar **a)
1822: {
1827: if (x->ops->restorearraywrite) {
1828: (*x->ops->restorearraywrite)(x,a);
1829: } else if (x->ops->restorearray) {
1830: (*x->ops->restorearray)(x,a);
1831: }
1832: if (a) *a = NULL;
1833: PetscObjectStateIncrease((PetscObject)x);
1834: return(0);
1835: }
1837: /*@C
1838: VecGetArrays - Returns a pointer to the arrays in a set of vectors
1839: that were created by a call to VecDuplicateVecs(). You MUST call
1840: VecRestoreArrays() when you no longer need access to the array.
1842: Logically Collective on Vec
1844: Input Parameters:
1845: + x - the vectors
1846: - n - the number of vectors
1848: Output Parameter:
1849: . a - location to put pointer to the array
1851: Fortran Note:
1852: This routine is not supported in Fortran.
1854: Level: intermediate
1856: .seealso: VecGetArray(), VecRestoreArrays()
1857: @*/
1858: PetscErrorCode VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1859: {
1861: PetscInt i;
1862: PetscScalar **q;
1868: if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1869: PetscMalloc1(n,&q);
1870: for (i=0; i<n; ++i) {
1871: VecGetArray(x[i],&q[i]);
1872: }
1873: *a = q;
1874: return(0);
1875: }
1877: /*@C
1878: VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1879: has been called.
1881: Logically Collective on Vec
1883: Input Parameters:
1884: + x - the vector
1885: . n - the number of vectors
1886: - a - location of pointer to arrays obtained from VecGetArrays()
1888: Notes:
1889: For regular PETSc vectors this routine does not involve any copies. For
1890: any special vectors that do not store local vector data in a contiguous
1891: array, this routine will copy the data back into the underlying
1892: vector data structure from the arrays obtained with VecGetArrays().
1894: Fortran Note:
1895: This routine is not supported in Fortran.
1897: Level: intermediate
1899: .seealso: VecGetArrays(), VecRestoreArray()
1900: @*/
1901: PetscErrorCode VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1902: {
1904: PetscInt i;
1905: PetscScalar **q = *a;
1912: for (i=0; i<n; ++i) {
1913: VecRestoreArray(x[i],&q[i]);
1914: }
1915: PetscFree(q);
1916: return(0);
1917: }
1919: /*@C
1920: VecGetArrayAndMemType - Like VecGetArray(), but if this is a GPU vector and it is currently offloaded to GPU,
1921: the returned pointer will be a GPU pointer to the GPU memory that contains this processor's portion of the
1922: vector data. Otherwise, it functions as VecGetArray().
1924: Logically Collective on Vec
1926: Input Parameter:
1927: . x - the vector
1929: Output Parameters:
1930: + a - location to put pointer to the array
1931: - mtype - memory type of the array
1933: Level: beginner
1935: .seealso: VecRestoreArrayAndMemType(), VecRestoreArrayAndMemType(), VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(),
1936: VecPlaceArray(), VecGetArray2d(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1937: @*/
1938: PetscErrorCode VecGetArrayAndMemType(Vec x,PetscScalar **a,PetscMemType *mtype)
1939: {
1945: VecSetErrorIfLocked(x,1);
1946: if (x->ops->getarrayandmemtype) { /* VECCUDA, VECKOKKOS etc */
1947: (*x->ops->getarrayandmemtype)(x,a,mtype);
1948: } else { /* VECSTANDARD, VECNEST, VECVIENNACL */
1949: VecGetArray(x,a);
1950: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
1951: }
1952: return(0);
1953: }
1955: /*@C
1956: VecRestoreArrayAndMemType - Restores a vector after VecGetArrayAndMemType() has been called.
1958: Logically Collective on Vec
1960: Input Parameters:
1961: + x - the vector
1962: - a - location of pointer to array obtained from VecGetArrayAndMemType()
1964: Level: beginner
1966: .seealso: VecGetArrayAndMemType(), VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(),
1967: VecPlaceArray(), VecRestoreArray2d(), VecGetArrayPair(), VecRestoreArrayPair()
1968: @*/
1969: PetscErrorCode VecRestoreArrayAndMemType(Vec x,PetscScalar **a)
1970: {
1976: if (x->ops->restorearrayandmemtype) { /* VECCUDA, VECKOKKOS etc */
1977: (*x->ops->restorearrayandmemtype)(x,a);
1978: } else if (x->ops->restorearray) { /* VECNEST, VECVIENNACL */
1979: (*x->ops->restorearray)(x,a);
1980: } /* VECSTANDARD does nothing */
1981: if (a) *a = NULL;
1982: PetscObjectStateIncrease((PetscObject)x);
1983: return(0);
1984: }
1986: /*@C
1987: VecGetArrayReadAndMemType - Like VecGetArrayRead(), but if this is a CUDA vector and it is currently offloaded to GPU,
1988: the returned pointer will be a GPU pointer to the GPU memory that contains this processor's portion of the
1989: vector data. Otherwise, it functions as VecGetArrayRead().
1991: Not Collective
1993: Input Parameter:
1994: . x - the vector
1996: Output Parameters:
1997: + a - the array
1998: - mtype - memory type of the array
2000: Level: beginner
2002: Notes:
2003: The array must be returned using a matching call to VecRestoreArrayReadAndMemType().
2006: .seealso: VecRestoreArrayReadAndMemType(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayAndMemType()
2007: @*/
2008: PetscErrorCode VecGetArrayReadAndMemType(Vec x,const PetscScalar **a,PetscMemType *mtype)
2009: {
2015: #if defined(PETSC_USE_DEBUG)
2016: if (x->ops->getarrayreadandmemtype) SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Not expected vector type \"%s\" has ops->getarrayreadandmemtype",((PetscObject)x)->type_name);
2017: #endif
2019: if (x->ops->getarrayandmemtype) { /* VECCUDA, VECKOKKOS etc, though they are also petscnative */
2020: (*x->ops->getarrayandmemtype)(x,(PetscScalar**)a,mtype);
2021: } else if (x->ops->getarray) { /* VECNEST, VECVIENNACL */
2022: (*x->ops->getarray)(x,(PetscScalar**)a);
2023: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2024: } else if (x->petscnative) { /* VECSTANDARD */
2025: *a = *((PetscScalar**)x->data);
2026: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2027: } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array read in place for vector type \"%s\"",((PetscObject)x)->type_name);
2028: return(0);
2029: }
2031: /*@C
2032: VecRestoreArrayReadAndMemType - Restore array obtained with VecGetArrayReadAndMemType()
2034: Not Collective
2036: Input Parameters:
2037: + vec - the vector
2038: - array - the array
2040: Level: beginner
2042: .seealso: VecGetArrayReadAndMemType(), VecRestoreArrayAndMemType(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
2043: @*/
2044: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x,const PetscScalar **a)
2045: {
2051: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS, VECVIENNACL etc */
2052: /* nothing */
2053: } else if (x->ops->restorearrayread) { /* VECNEST */
2054: (*x->ops->restorearrayread)(x,a);
2055: } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot restore array read in place for vector type \"%s\"",((PetscObject)x)->type_name);
2056: if (a) *a = NULL;
2057: return(0);
2058: }
2060: /*@
2061: VecPlaceArray - Allows one to replace the array in a vector with an
2062: array provided by the user. This is useful to avoid copying an array
2063: into a vector.
2065: Not Collective
2067: Input Parameters:
2068: + vec - the vector
2069: - array - the array
2071: Notes:
2072: You can return to the original array with a call to VecResetArray()
2074: Level: developer
2076: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecResetArray()
2078: @*/
2079: PetscErrorCode VecPlaceArray(Vec vec,const PetscScalar array[])
2080: {
2087: if (vec->ops->placearray) {
2088: (*vec->ops->placearray)(vec,array);
2089: } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
2090: PetscObjectStateIncrease((PetscObject)vec);
2091: return(0);
2092: }
2094: /*@C
2095: VecReplaceArray - Allows one to replace the array in a vector with an
2096: array provided by the user. This is useful to avoid copying an array
2097: into a vector.
2099: Not Collective
2101: Input Parameters:
2102: + vec - the vector
2103: - array - the array
2105: Notes:
2106: This permanently replaces the array and frees the memory associated
2107: with the old array.
2109: The memory passed in MUST be obtained with PetscMalloc() and CANNOT be
2110: freed by the user. It will be freed when the vector is destroyed.
2112: Not supported from Fortran
2114: Level: developer
2116: .seealso: VecGetArray(), VecRestoreArray(), VecPlaceArray(), VecResetArray()
2118: @*/
2119: PetscErrorCode VecReplaceArray(Vec vec,const PetscScalar array[])
2120: {
2126: if (vec->ops->replacearray) {
2127: (*vec->ops->replacearray)(vec,array);
2128: } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
2129: PetscObjectStateIncrease((PetscObject)vec);
2130: return(0);
2131: }
2134: /*@C
2135: VecCUDAGetArray - Provides access to the CUDA buffer inside a vector.
2137: This function has semantics similar to VecGetArray(): the pointer
2138: returned by this function points to a consistent view of the vector
2139: data. This may involve a copy operation of data from the host to the
2140: device if the data on the device is out of date. If the device
2141: memory hasn't been allocated previously it will be allocated as part
2142: of this function call. VecCUDAGetArray() assumes that
2143: the user will modify the vector data. This is similar to
2144: intent(inout) in fortran.
2146: The CUDA device pointer has to be released by calling
2147: VecCUDARestoreArray(). Upon restoring the vector data
2148: the data on the host will be marked as out of date. A subsequent
2149: access of the host data will thus incur a data transfer from the
2150: device to the host.
2153: Input Parameter:
2154: . v - the vector
2156: Output Parameter:
2157: . a - the CUDA device pointer
2159: Fortran note:
2160: This function is not currently available from Fortran.
2162: Level: intermediate
2164: .seealso: VecCUDARestoreArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2165: @*/
2166: PETSC_EXTERN PetscErrorCode VecCUDAGetArray(Vec v, PetscScalar **a)
2167: {
2170: #if defined(PETSC_HAVE_CUDA)
2171: {
2173: VecCUDACopyToGPU(v);
2174: *a = ((Vec_CUDA*)v->spptr)->GPUarray;
2175: }
2176: #endif
2177: return(0);
2178: }
2180: /*@C
2181: VecCUDARestoreArray - Restore a CUDA device pointer previously acquired with VecCUDAGetArray().
2183: This marks the host data as out of date. Subsequent access to the
2184: vector data on the host side with for instance VecGetArray() incurs a
2185: data transfer.
2187: Input Parameter:
2188: + v - the vector
2189: - a - the CUDA device pointer. This pointer is invalid after
2190: VecCUDARestoreArray() returns.
2192: Fortran note:
2193: This function is not currently available from Fortran.
2195: Level: intermediate
2197: .seealso: VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2198: @*/
2199: PETSC_EXTERN PetscErrorCode VecCUDARestoreArray(Vec v, PetscScalar **a)
2200: {
2205: #if defined(PETSC_HAVE_CUDA)
2206: v->offloadmask = PETSC_OFFLOAD_GPU;
2207: #endif
2208: PetscObjectStateIncrease((PetscObject)v);
2209: return(0);
2210: }
2212: /*@C
2213: VecCUDAGetArrayRead - Provides read access to the CUDA buffer inside a vector.
2215: This function is analogous to VecGetArrayRead(): The pointer
2216: returned by this function points to a consistent view of the vector
2217: data. This may involve a copy operation of data from the host to the
2218: device if the data on the device is out of date. If the device
2219: memory hasn't been allocated previously it will be allocated as part
2220: of this function call. VecCUDAGetArrayRead() assumes that the
2221: user will not modify the vector data. This is analgogous to
2222: intent(in) in Fortran.
2224: The CUDA device pointer has to be released by calling
2225: VecCUDARestoreArrayRead(). If the data on the host side was
2226: previously up to date it will remain so, i.e. data on both the device
2227: and the host is up to date. Accessing data on the host side does not
2228: incur a device to host data transfer.
2230: Input Parameter:
2231: . v - the vector
2233: Output Parameter:
2234: . a - the CUDA pointer.
2236: Fortran note:
2237: This function is not currently available from Fortran.
2239: Level: intermediate
2241: .seealso: VecCUDARestoreArrayRead(), VecCUDAGetArray(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2242: @*/
2243: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayRead(Vec v,const PetscScalar** a)
2244: {
2247: VecCUDAGetArray(v,(PetscScalar**)a);
2248: return(0);
2249: }
2251: /*@C
2252: VecCUDARestoreArrayRead - Restore a CUDA device pointer previously acquired with VecCUDAGetArrayRead().
2254: If the data on the host side was previously up to date it will remain
2255: so, i.e. data on both the device and the host is up to date.
2256: Accessing data on the host side e.g. with VecGetArray() does not
2257: incur a device to host data transfer.
2259: Input Parameter:
2260: + v - the vector
2261: - a - the CUDA device pointer. This pointer is invalid after
2262: VecCUDARestoreArrayRead() returns.
2264: Fortran note:
2265: This function is not currently available from Fortran.
2267: Level: intermediate
2269: .seealso: VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2270: @*/
2271: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayRead(Vec v, const PetscScalar **a)
2272: {
2275: *a = NULL;
2276: return(0);
2277: }
2279: /*@C
2280: VecCUDAGetArrayWrite - Provides write access to the CUDA buffer inside a vector.
2282: The data pointed to by the device pointer is uninitialized. The user
2283: may not read from this data. Furthermore, the entire array needs to
2284: be filled by the user to obtain well-defined behaviour. The device
2285: memory will be allocated by this function if it hasn't been allocated
2286: previously. This is analogous to intent(out) in Fortran.
2288: The device pointer needs to be released with
2289: VecCUDARestoreArrayWrite(). When the pointer is released the
2290: host data of the vector is marked as out of data. Subsequent access
2291: of the host data with e.g. VecGetArray() incurs a device to host data
2292: transfer.
2295: Input Parameter:
2296: . v - the vector
2298: Output Parameter:
2299: . a - the CUDA pointer
2301: Fortran note:
2302: This function is not currently available from Fortran.
2304: Level: advanced
2306: .seealso: VecCUDARestoreArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2307: @*/
2308: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayWrite(Vec v, PetscScalar **a)
2309: {
2312: #if defined(PETSC_HAVE_CUDA)
2313: {
2315: VecCUDAAllocateCheck(v);
2316: *a = ((Vec_CUDA*)v->spptr)->GPUarray;
2317: }
2318: #endif
2319: return(0);
2320: }
2322: /*@C
2323: VecCUDARestoreArrayWrite - Restore a CUDA device pointer previously acquired with VecCUDAGetArrayWrite().
2325: Data on the host will be marked as out of date. Subsequent access of
2326: the data on the host side e.g. with VecGetArray() will incur a device
2327: to host data transfer.
2329: Input Parameter:
2330: + v - the vector
2331: - a - the CUDA device pointer. This pointer is invalid after
2332: VecCUDARestoreArrayWrite() returns.
2334: Fortran note:
2335: This function is not currently available from Fortran.
2337: Level: intermediate
2339: .seealso: VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2340: @*/
2341: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayWrite(Vec v, PetscScalar **a)
2342: {
2347: #if defined(PETSC_HAVE_CUDA)
2348: v->offloadmask = PETSC_OFFLOAD_GPU;
2349: a = NULL;
2350: #endif
2351: PetscObjectStateIncrease((PetscObject)v);
2352: return(0);
2353: }
2355: /*@C
2356: VecCUDAPlaceArray - Allows one to replace the GPU array in a vector with a
2357: GPU array provided by the user. This is useful to avoid copying an
2358: array into a vector.
2360: Not Collective
2362: Input Parameters:
2363: + vec - the vector
2364: - array - the GPU array
2366: Notes:
2367: You can return to the original GPU array with a call to VecCUDAResetArray()
2368: It is not possible to use VecCUDAPlaceArray() and VecPlaceArray() at the
2369: same time on the same vector.
2371: Level: developer
2373: .seealso: VecPlaceArray(), VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecResetArray(), VecCUDAResetArray(), VecCUDAReplaceArray()
2375: @*/
2376: PetscErrorCode VecCUDAPlaceArray(Vec vin,const PetscScalar a[])
2377: {
2382: #if defined(PETSC_HAVE_CUDA)
2383: VecCUDACopyToGPU(vin);
2384: if (((Vec_Seq*)vin->data)->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"VecCUDAPlaceArray()/VecPlaceArray() was already called on this vector, without a call to VecCUDAResetArray()/VecResetArray()");
2385: ((Vec_Seq*)vin->data)->unplacedarray = (PetscScalar *) ((Vec_CUDA*)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2386: ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2387: vin->offloadmask = PETSC_OFFLOAD_GPU;
2388: #endif
2389: PetscObjectStateIncrease((PetscObject)vin);
2390: return(0);
2391: }
2393: /*@C
2394: VecCUDAReplaceArray - Allows one to replace the GPU array in a vector
2395: with a GPU array provided by the user. This is useful to avoid copying
2396: a GPU array into a vector.
2398: Not Collective
2400: Input Parameters:
2401: + vec - the vector
2402: - array - the GPU array
2404: Notes:
2405: This permanently replaces the GPU array and frees the memory associated
2406: with the old GPU array.
2408: The memory passed in CANNOT be freed by the user. It will be freed
2409: when the vector is destroyed.
2411: Not supported from Fortran
2413: Level: developer
2415: .seealso: VecGetArray(), VecRestoreArray(), VecPlaceArray(), VecResetArray(), VecCUDAResetArray(), VecCUDAPlaceArray(), VecReplaceArray()
2417: @*/
2418: PetscErrorCode VecCUDAReplaceArray(Vec vin,const PetscScalar a[])
2419: {
2420: #if defined(PETSC_HAVE_CUDA)
2421: cudaError_t err;
2422: #endif
2427: #if defined(PETSC_HAVE_CUDA)
2428: err = cudaFree(((Vec_CUDA*)vin->spptr)->GPUarray);CHKERRCUDA(err);
2429: ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2430: vin->offloadmask = PETSC_OFFLOAD_GPU;
2431: #endif
2432: PetscObjectStateIncrease((PetscObject)vin);
2433: return(0);
2434: }
2436: /*@C
2437: VecCUDAResetArray - Resets a vector to use its default memory. Call this
2438: after the use of VecCUDAPlaceArray().
2440: Not Collective
2442: Input Parameters:
2443: . vec - the vector
2445: Level: developer
2447: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecPlaceArray(), VecResetArray(), VecCUDAPlaceArray(), VecCUDAReplaceArray()
2449: @*/
2450: PetscErrorCode VecCUDAResetArray(Vec vin)
2451: {
2456: #if defined(PETSC_HAVE_CUDA)
2457: VecCUDACopyToGPU(vin);
2458: ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar *) ((Vec_Seq*)vin->data)->unplacedarray;
2459: ((Vec_Seq*)vin->data)->unplacedarray = 0;
2460: vin->offloadmask = PETSC_OFFLOAD_GPU;
2461: #endif
2462: PetscObjectStateIncrease((PetscObject)vin);
2463: return(0);
2464: }
2466: /*MC
2467: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
2468: and makes them accessible via a Fortran90 pointer.
2470: Synopsis:
2471: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
2473: Collective on Vec
2475: Input Parameters:
2476: + x - a vector to mimic
2477: - n - the number of vectors to obtain
2479: Output Parameters:
2480: + y - Fortran90 pointer to the array of vectors
2481: - ierr - error code
2483: Example of Usage:
2484: .vb
2485: #include <petsc/finclude/petscvec.h>
2486: use petscvec
2488: Vec x
2489: Vec, pointer :: y(:)
2490: ....
2491: call VecDuplicateVecsF90(x,2,y,ierr)
2492: call VecSet(y(2),alpha,ierr)
2493: call VecSet(y(2),alpha,ierr)
2494: ....
2495: call VecDestroyVecsF90(2,y,ierr)
2496: .ve
2498: Notes:
2499: Not yet supported for all F90 compilers
2501: Use VecDestroyVecsF90() to free the space.
2503: Level: beginner
2505: .seealso: VecDestroyVecsF90(), VecDuplicateVecs()
2507: M*/
2509: /*MC
2510: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2511: VecGetArrayF90().
2513: Synopsis:
2514: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2516: Logically Collective on Vec
2518: Input Parameters:
2519: + x - vector
2520: - xx_v - the Fortran90 pointer to the array
2522: Output Parameter:
2523: . ierr - error code
2525: Example of Usage:
2526: .vb
2527: #include <petsc/finclude/petscvec.h>
2528: use petscvec
2530: PetscScalar, pointer :: xx_v(:)
2531: ....
2532: call VecGetArrayF90(x,xx_v,ierr)
2533: xx_v(3) = a
2534: call VecRestoreArrayF90(x,xx_v,ierr)
2535: .ve
2537: Level: beginner
2539: .seealso: VecGetArrayF90(), VecGetArray(), VecRestoreArray(), UsingFortran, VecRestoreArrayReadF90()
2541: M*/
2543: /*MC
2544: VecDestroyVecsF90 - Frees a block of vectors obtained with VecDuplicateVecsF90().
2546: Synopsis:
2547: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
2549: Collective on Vec
2551: Input Parameters:
2552: + n - the number of vectors previously obtained
2553: - x - pointer to array of vector pointers
2555: Output Parameter:
2556: . ierr - error code
2558: Notes:
2559: Not yet supported for all F90 compilers
2561: Level: beginner
2563: .seealso: VecDestroyVecs(), VecDuplicateVecsF90()
2565: M*/
2567: /*MC
2568: VecGetArrayF90 - Accesses a vector array from Fortran90. For default PETSc
2569: vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
2570: this routine is implementation dependent. You MUST call VecRestoreArrayF90()
2571: when you no longer need access to the array.
2573: Synopsis:
2574: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2576: Logically Collective on Vec
2578: Input Parameter:
2579: . x - vector
2581: Output Parameters:
2582: + xx_v - the Fortran90 pointer to the array
2583: - ierr - error code
2585: Example of Usage:
2586: .vb
2587: #include <petsc/finclude/petscvec.h>
2588: use petscvec
2590: PetscScalar, pointer :: xx_v(:)
2591: ....
2592: call VecGetArrayF90(x,xx_v,ierr)
2593: xx_v(3) = a
2594: call VecRestoreArrayF90(x,xx_v,ierr)
2595: .ve
2597: If you ONLY intend to read entries from the array and not change any entries you should use VecGetArrayReadF90().
2599: Level: beginner
2601: .seealso: VecRestoreArrayF90(), VecGetArray(), VecRestoreArray(), VecGetArrayReadF90(), UsingFortran
2603: M*/
2605: /*MC
2606: VecGetArrayReadF90 - Accesses a read only array from Fortran90. For default PETSc
2607: vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
2608: this routine is implementation dependent. You MUST call VecRestoreArrayReadF90()
2609: when you no longer need access to the array.
2611: Synopsis:
2612: VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2614: Logically Collective on Vec
2616: Input Parameter:
2617: . x - vector
2619: Output Parameters:
2620: + xx_v - the Fortran90 pointer to the array
2621: - ierr - error code
2623: Example of Usage:
2624: .vb
2625: #include <petsc/finclude/petscvec.h>
2626: use petscvec
2628: PetscScalar, pointer :: xx_v(:)
2629: ....
2630: call VecGetArrayReadF90(x,xx_v,ierr)
2631: a = xx_v(3)
2632: call VecRestoreArrayReadF90(x,xx_v,ierr)
2633: .ve
2635: If you intend to write entries into the array you must use VecGetArrayF90().
2637: Level: beginner
2639: .seealso: VecRestoreArrayReadF90(), VecGetArray(), VecRestoreArray(), VecGetArrayRead(), VecRestoreArrayRead(), VecGetArrayF90(), UsingFortran
2641: M*/
2643: /*MC
2644: VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2645: VecGetArrayReadF90().
2647: Synopsis:
2648: VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2650: Logically Collective on Vec
2652: Input Parameters:
2653: + x - vector
2654: - xx_v - the Fortran90 pointer to the array
2656: Output Parameter:
2657: . ierr - error code
2659: Example of Usage:
2660: .vb
2661: #include <petsc/finclude/petscvec.h>
2662: use petscvec
2664: PetscScalar, pointer :: xx_v(:)
2665: ....
2666: call VecGetArrayReadF90(x,xx_v,ierr)
2667: a = xx_v(3)
2668: call VecRestoreArrayReadF90(x,xx_v,ierr)
2669: .ve
2671: Level: beginner
2673: .seealso: VecGetArrayReadF90(), VecGetArray(), VecRestoreArray(), VecGetArrayRead(), VecRestoreArrayRead(),UsingFortran, VecRestoreArrayF90()
2675: M*/
2677: /*@C
2678: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2679: processor's portion of the vector data. You MUST call VecRestoreArray2d()
2680: when you no longer need access to the array.
2682: Logically Collective
2684: Input Parameter:
2685: + x - the vector
2686: . m - first dimension of two dimensional array
2687: . n - second dimension of two dimensional array
2688: . mstart - first index you will use in first coordinate direction (often 0)
2689: - nstart - first index in the second coordinate direction (often 0)
2691: Output Parameter:
2692: . a - location to put pointer to the array
2694: Level: developer
2696: Notes:
2697: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2698: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2699: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2700: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
2702: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2704: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2705: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2706: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2707: @*/
2708: PetscErrorCode VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2709: {
2711: PetscInt i,N;
2712: PetscScalar *aa;
2718: VecGetLocalSize(x,&N);
2719: 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);
2720: VecGetArray(x,&aa);
2722: PetscMalloc1(m,a);
2723: for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2724: *a -= mstart;
2725: return(0);
2726: }
2728: /*@C
2729: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2730: processor's portion of the vector data. You MUST call VecRestoreArray2dWrite()
2731: when you no longer need access to the array.
2733: Logically Collective
2735: Input Parameter:
2736: + x - the vector
2737: . m - first dimension of two dimensional array
2738: . n - second dimension of two dimensional array
2739: . mstart - first index you will use in first coordinate direction (often 0)
2740: - nstart - first index in the second coordinate direction (often 0)
2742: Output Parameter:
2743: . a - location to put pointer to the array
2745: Level: developer
2747: Notes:
2748: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2749: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2750: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2751: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
2753: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2755: Concepts: vector^accessing local values as 2d array
2757: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2758: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2759: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2760: @*/
2761: PetscErrorCode VecGetArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2762: {
2764: PetscInt i,N;
2765: PetscScalar *aa;
2771: VecGetLocalSize(x,&N);
2772: 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);
2773: VecGetArrayWrite(x,&aa);
2775: PetscMalloc1(m,a);
2776: for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2777: *a -= mstart;
2778: return(0);
2779: }
2781: /*@C
2782: VecRestoreArray2d - Restores a vector after VecGetArray2d() has been called.
2784: Logically Collective
2786: Input Parameters:
2787: + x - the vector
2788: . m - first dimension of two dimensional array
2789: . n - second dimension of the two dimensional array
2790: . mstart - first index you will use in first coordinate direction (often 0)
2791: . nstart - first index in the second coordinate direction (often 0)
2792: - a - location of pointer to array obtained from VecGetArray2d()
2794: Level: developer
2796: Notes:
2797: For regular PETSc vectors this routine does not involve any copies. For
2798: any special vectors that do not store local vector data in a contiguous
2799: array, this routine will copy the data back into the underlying
2800: vector data structure from the array obtained with VecGetArray().
2802: This routine actually zeros out the a pointer.
2804: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2805: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2806: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2807: @*/
2808: PetscErrorCode VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2809: {
2811: void *dummy;
2817: dummy = (void*)(*a + mstart);
2818: PetscFree(dummy);
2819: VecRestoreArray(x,NULL);
2820: return(0);
2821: }
2823: /*@C
2824: VecRestoreArray2dWrite - Restores a vector after VecGetArray2dWrite() has been called.
2826: Logically Collective
2828: Input Parameters:
2829: + x - the vector
2830: . m - first dimension of two dimensional array
2831: . n - second dimension of the two dimensional array
2832: . mstart - first index you will use in first coordinate direction (often 0)
2833: . nstart - first index in the second coordinate direction (often 0)
2834: - a - location of pointer to array obtained from VecGetArray2d()
2836: Level: developer
2838: Notes:
2839: For regular PETSc vectors this routine does not involve any copies. For
2840: any special vectors that do not store local vector data in a contiguous
2841: array, this routine will copy the data back into the underlying
2842: vector data structure from the array obtained with VecGetArray().
2844: This routine actually zeros out the a pointer.
2846: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2847: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2848: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2849: @*/
2850: PetscErrorCode VecRestoreArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2851: {
2853: void *dummy;
2859: dummy = (void*)(*a + mstart);
2860: PetscFree(dummy);
2861: VecRestoreArrayWrite(x,NULL);
2862: return(0);
2863: }
2865: /*@C
2866: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2867: processor's portion of the vector data. You MUST call VecRestoreArray1d()
2868: when you no longer need access to the array.
2870: Logically Collective
2872: Input Parameter:
2873: + x - the vector
2874: . m - first dimension of two dimensional array
2875: - mstart - first index you will use in first coordinate direction (often 0)
2877: Output Parameter:
2878: . a - location to put pointer to the array
2880: Level: developer
2882: Notes:
2883: For a vector obtained from DMCreateLocalVector() mstart are likely
2884: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2885: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
2887: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2889: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2890: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2891: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2892: @*/
2893: PetscErrorCode VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2894: {
2896: PetscInt N;
2902: VecGetLocalSize(x,&N);
2903: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2904: VecGetArray(x,a);
2905: *a -= mstart;
2906: return(0);
2907: }
2909: /*@C
2910: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
2911: processor's portion of the vector data. You MUST call VecRestoreArray1dWrite()
2912: when you no longer need access to the array.
2914: Logically Collective
2916: Input Parameter:
2917: + x - the vector
2918: . m - first dimension of two dimensional array
2919: - mstart - first index you will use in first coordinate direction (often 0)
2921: Output Parameter:
2922: . a - location to put pointer to the array
2924: Level: developer
2926: Notes:
2927: For a vector obtained from DMCreateLocalVector() mstart are likely
2928: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2929: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
2931: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2933: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2934: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2935: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2936: @*/
2937: PetscErrorCode VecGetArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2938: {
2940: PetscInt N;
2946: VecGetLocalSize(x,&N);
2947: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2948: VecGetArrayWrite(x,a);
2949: *a -= mstart;
2950: return(0);
2951: }
2953: /*@C
2954: VecRestoreArray1d - Restores a vector after VecGetArray1d() has been called.
2956: Logically Collective
2958: Input Parameters:
2959: + x - the vector
2960: . m - first dimension of two dimensional array
2961: . mstart - first index you will use in first coordinate direction (often 0)
2962: - a - location of pointer to array obtained from VecGetArray21()
2964: Level: developer
2966: Notes:
2967: For regular PETSc vectors this routine does not involve any copies. For
2968: any special vectors that do not store local vector data in a contiguous
2969: array, this routine will copy the data back into the underlying
2970: vector data structure from the array obtained with VecGetArray1d().
2972: This routine actually zeros out the a pointer.
2974: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2975: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2976: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2977: @*/
2978: PetscErrorCode VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2979: {
2985: VecRestoreArray(x,NULL);
2986: return(0);
2987: }
2989: /*@C
2990: VecRestoreArray1dWrite - Restores a vector after VecGetArray1dWrite() has been called.
2992: Logically Collective
2994: Input Parameters:
2995: + x - the vector
2996: . m - first dimension of two dimensional array
2997: . mstart - first index you will use in first coordinate direction (often 0)
2998: - a - location of pointer to array obtained from VecGetArray21()
3000: Level: developer
3002: Notes:
3003: For regular PETSc vectors this routine does not involve any copies. For
3004: any special vectors that do not store local vector data in a contiguous
3005: array, this routine will copy the data back into the underlying
3006: vector data structure from the array obtained with VecGetArray1d().
3008: This routine actually zeros out the a pointer.
3010: Concepts: vector^accessing local values as 1d array
3012: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3013: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3014: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3015: @*/
3016: PetscErrorCode VecRestoreArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3017: {
3023: VecRestoreArrayWrite(x,NULL);
3024: return(0);
3025: }
3027: /*@C
3028: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
3029: processor's portion of the vector data. You MUST call VecRestoreArray3d()
3030: when you no longer need access to the array.
3032: Logically Collective
3034: Input Parameter:
3035: + x - the vector
3036: . m - first dimension of three dimensional array
3037: . n - second dimension of three dimensional array
3038: . p - third dimension of three dimensional array
3039: . mstart - first index you will use in first coordinate direction (often 0)
3040: . nstart - first index in the second coordinate direction (often 0)
3041: - pstart - first index in the third coordinate direction (often 0)
3043: Output Parameter:
3044: . a - location to put pointer to the array
3046: Level: developer
3048: Notes:
3049: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3050: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3051: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3052: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
3054: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3056: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3057: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3058: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3059: @*/
3060: PetscErrorCode VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3061: {
3063: PetscInt i,N,j;
3064: PetscScalar *aa,**b;
3070: VecGetLocalSize(x,&N);
3071: 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);
3072: VecGetArray(x,&aa);
3074: PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3075: b = (PetscScalar**)((*a) + m);
3076: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3077: for (i=0; i<m; i++)
3078: for (j=0; j<n; j++)
3079: b[i*n+j] = aa + i*n*p + j*p - pstart;
3081: *a -= mstart;
3082: return(0);
3083: }
3085: /*@C
3086: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
3087: processor's portion of the vector data. You MUST call VecRestoreArray3dWrite()
3088: when you no longer need access to the array.
3090: Logically Collective
3092: Input Parameter:
3093: + x - the vector
3094: . m - first dimension of three dimensional array
3095: . n - second dimension of three dimensional array
3096: . p - third dimension of three dimensional array
3097: . mstart - first index you will use in first coordinate direction (often 0)
3098: . nstart - first index in the second coordinate direction (often 0)
3099: - pstart - first index in the third coordinate direction (often 0)
3101: Output Parameter:
3102: . a - location to put pointer to the array
3104: Level: developer
3106: Notes:
3107: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3108: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3109: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3110: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
3112: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3114: Concepts: vector^accessing local values as 3d array
3116: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3117: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3118: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3119: @*/
3120: PetscErrorCode VecGetArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3121: {
3123: PetscInt i,N,j;
3124: PetscScalar *aa,**b;
3130: VecGetLocalSize(x,&N);
3131: 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);
3132: VecGetArrayWrite(x,&aa);
3134: PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3135: b = (PetscScalar**)((*a) + m);
3136: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3137: for (i=0; i<m; i++)
3138: for (j=0; j<n; j++)
3139: b[i*n+j] = aa + i*n*p + j*p - pstart;
3141: *a -= mstart;
3142: return(0);
3143: }
3145: /*@C
3146: VecRestoreArray3d - Restores a vector after VecGetArray3d() has been called.
3148: Logically Collective
3150: Input Parameters:
3151: + x - the vector
3152: . m - first dimension of three dimensional array
3153: . n - second dimension of the three dimensional array
3154: . p - third dimension of the three dimensional array
3155: . mstart - first index you will use in first coordinate direction (often 0)
3156: . nstart - first index in the second coordinate direction (often 0)
3157: . pstart - first index in the third coordinate direction (often 0)
3158: - a - location of pointer to array obtained from VecGetArray3d()
3160: Level: developer
3162: Notes:
3163: For regular PETSc vectors this routine does not involve any copies. For
3164: any special vectors that do not store local vector data in a contiguous
3165: array, this routine will copy the data back into the underlying
3166: vector data structure from the array obtained with VecGetArray().
3168: This routine actually zeros out the a pointer.
3170: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3171: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3172: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3173: @*/
3174: PetscErrorCode VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3175: {
3177: void *dummy;
3183: dummy = (void*)(*a + mstart);
3184: PetscFree(dummy);
3185: VecRestoreArray(x,NULL);
3186: return(0);
3187: }
3189: /*@C
3190: VecRestoreArray3dWrite - Restores a vector after VecGetArray3dWrite() has been called.
3192: Logically Collective
3194: Input Parameters:
3195: + x - the vector
3196: . m - first dimension of three dimensional array
3197: . n - second dimension of the three dimensional array
3198: . p - third dimension of the three dimensional array
3199: . mstart - first index you will use in first coordinate direction (often 0)
3200: . nstart - first index in the second coordinate direction (often 0)
3201: . pstart - first index in the third coordinate direction (often 0)
3202: - a - location of pointer to array obtained from VecGetArray3d()
3204: Level: developer
3206: Notes:
3207: For regular PETSc vectors this routine does not involve any copies. For
3208: any special vectors that do not store local vector data in a contiguous
3209: array, this routine will copy the data back into the underlying
3210: vector data structure from the array obtained with VecGetArray().
3212: This routine actually zeros out the a pointer.
3214: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3215: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3216: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3217: @*/
3218: PetscErrorCode VecRestoreArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3219: {
3221: void *dummy;
3227: dummy = (void*)(*a + mstart);
3228: PetscFree(dummy);
3229: VecRestoreArrayWrite(x,NULL);
3230: return(0);
3231: }
3233: /*@C
3234: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3235: processor's portion of the vector data. You MUST call VecRestoreArray4d()
3236: when you no longer need access to the array.
3238: Logically Collective
3240: Input Parameter:
3241: + x - the vector
3242: . m - first dimension of four dimensional array
3243: . n - second dimension of four dimensional array
3244: . p - third dimension of four dimensional array
3245: . q - fourth dimension of four dimensional array
3246: . mstart - first index you will use in first coordinate direction (often 0)
3247: . nstart - first index in the second coordinate direction (often 0)
3248: . pstart - first index in the third coordinate direction (often 0)
3249: - qstart - first index in the fourth coordinate direction (often 0)
3251: Output Parameter:
3252: . a - location to put pointer to the array
3254: Level: beginner
3256: Notes:
3257: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3258: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3259: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3260: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
3262: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3264: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3265: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3266: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3267: @*/
3268: PetscErrorCode VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3269: {
3271: PetscInt i,N,j,k;
3272: PetscScalar *aa,***b,**c;
3278: VecGetLocalSize(x,&N);
3279: 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);
3280: VecGetArray(x,&aa);
3282: PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3283: b = (PetscScalar***)((*a) + m);
3284: c = (PetscScalar**)(b + m*n);
3285: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3286: for (i=0; i<m; i++)
3287: for (j=0; j<n; j++)
3288: b[i*n+j] = c + i*n*p + j*p - pstart;
3289: for (i=0; i<m; i++)
3290: for (j=0; j<n; j++)
3291: for (k=0; k<p; k++)
3292: c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3293: *a -= mstart;
3294: return(0);
3295: }
3297: /*@C
3298: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3299: processor's portion of the vector data. You MUST call VecRestoreArray4dWrite()
3300: when you no longer need access to the array.
3302: Logically Collective
3304: Input Parameter:
3305: + x - the vector
3306: . m - first dimension of four dimensional array
3307: . n - second dimension of four dimensional array
3308: . p - third dimension of four dimensional array
3309: . q - fourth dimension of four dimensional array
3310: . mstart - first index you will use in first coordinate direction (often 0)
3311: . nstart - first index in the second coordinate direction (often 0)
3312: . pstart - first index in the third coordinate direction (often 0)
3313: - qstart - first index in the fourth coordinate direction (often 0)
3315: Output Parameter:
3316: . a - location to put pointer to the array
3318: Level: beginner
3320: Notes:
3321: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3322: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3323: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3324: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
3326: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3328: Concepts: vector^accessing local values as 3d array
3330: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3331: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3332: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3333: @*/
3334: PetscErrorCode VecGetArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3335: {
3337: PetscInt i,N,j,k;
3338: PetscScalar *aa,***b,**c;
3344: VecGetLocalSize(x,&N);
3345: 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);
3346: VecGetArrayWrite(x,&aa);
3348: PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3349: b = (PetscScalar***)((*a) + m);
3350: c = (PetscScalar**)(b + m*n);
3351: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3352: for (i=0; i<m; i++)
3353: for (j=0; j<n; j++)
3354: b[i*n+j] = c + i*n*p + j*p - pstart;
3355: for (i=0; i<m; i++)
3356: for (j=0; j<n; j++)
3357: for (k=0; k<p; k++)
3358: c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3359: *a -= mstart;
3360: return(0);
3361: }
3363: /*@C
3364: VecRestoreArray4d - Restores a vector after VecGetArray3d() has been called.
3366: Logically Collective
3368: Input Parameters:
3369: + x - the vector
3370: . m - first dimension of four dimensional array
3371: . n - second dimension of the four dimensional array
3372: . p - third dimension of the four dimensional array
3373: . q - fourth dimension of the four dimensional array
3374: . mstart - first index you will use in first coordinate direction (often 0)
3375: . nstart - first index in the second coordinate direction (often 0)
3376: . pstart - first index in the third coordinate direction (often 0)
3377: . qstart - first index in the fourth coordinate direction (often 0)
3378: - a - location of pointer to array obtained from VecGetArray4d()
3380: Level: beginner
3382: Notes:
3383: For regular PETSc vectors this routine does not involve any copies. For
3384: any special vectors that do not store local vector data in a contiguous
3385: array, this routine will copy the data back into the underlying
3386: vector data structure from the array obtained with VecGetArray().
3388: This routine actually zeros out the a pointer.
3390: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3391: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3392: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3393: @*/
3394: PetscErrorCode VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3395: {
3397: void *dummy;
3403: dummy = (void*)(*a + mstart);
3404: PetscFree(dummy);
3405: VecRestoreArray(x,NULL);
3406: return(0);
3407: }
3409: /*@C
3410: VecRestoreArray4dWrite - Restores a vector after VecGetArray3dWrite() has been called.
3412: Logically Collective
3414: Input Parameters:
3415: + x - the vector
3416: . m - first dimension of four dimensional array
3417: . n - second dimension of the four dimensional array
3418: . p - third dimension of the four dimensional array
3419: . q - fourth dimension of the four dimensional array
3420: . mstart - first index you will use in first coordinate direction (often 0)
3421: . nstart - first index in the second coordinate direction (often 0)
3422: . pstart - first index in the third coordinate direction (often 0)
3423: . qstart - first index in the fourth coordinate direction (often 0)
3424: - a - location of pointer to array obtained from VecGetArray4d()
3426: Level: beginner
3428: Notes:
3429: For regular PETSc vectors this routine does not involve any copies. For
3430: any special vectors that do not store local vector data in a contiguous
3431: array, this routine will copy the data back into the underlying
3432: vector data structure from the array obtained with VecGetArray().
3434: This routine actually zeros out the a pointer.
3436: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3437: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3438: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3439: @*/
3440: PetscErrorCode VecRestoreArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3441: {
3443: void *dummy;
3449: dummy = (void*)(*a + mstart);
3450: PetscFree(dummy);
3451: VecRestoreArrayWrite(x,NULL);
3452: return(0);
3453: }
3455: /*@C
3456: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3457: processor's portion of the vector data. You MUST call VecRestoreArray2dRead()
3458: when you no longer need access to the array.
3460: Logically Collective
3462: Input Parameter:
3463: + x - the vector
3464: . m - first dimension of two dimensional array
3465: . n - second dimension of two dimensional array
3466: . mstart - first index you will use in first coordinate direction (often 0)
3467: - nstart - first index in the second coordinate direction (often 0)
3469: Output Parameter:
3470: . a - location to put pointer to the array
3472: Level: developer
3474: Notes:
3475: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
3476: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3477: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3478: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
3480: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3482: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3483: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3484: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3485: @*/
3486: PetscErrorCode VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3487: {
3488: PetscErrorCode ierr;
3489: PetscInt i,N;
3490: const PetscScalar *aa;
3496: VecGetLocalSize(x,&N);
3497: 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);
3498: VecGetArrayRead(x,&aa);
3500: PetscMalloc1(m,a);
3501: for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
3502: *a -= mstart;
3503: return(0);
3504: }
3506: /*@C
3507: VecRestoreArray2dRead - Restores a vector after VecGetArray2dRead() has been called.
3509: Logically Collective
3511: Input Parameters:
3512: + x - the vector
3513: . m - first dimension of two dimensional array
3514: . n - second dimension of the two dimensional array
3515: . mstart - first index you will use in first coordinate direction (often 0)
3516: . nstart - first index in the second coordinate direction (often 0)
3517: - a - location of pointer to array obtained from VecGetArray2d()
3519: Level: developer
3521: Notes:
3522: For regular PETSc vectors this routine does not involve any copies. For
3523: any special vectors that do not store local vector data in a contiguous
3524: array, this routine will copy the data back into the underlying
3525: vector data structure from the array obtained with VecGetArray().
3527: This routine actually zeros out the a pointer.
3529: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3530: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3531: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3532: @*/
3533: PetscErrorCode VecRestoreArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3534: {
3536: void *dummy;
3542: dummy = (void*)(*a + mstart);
3543: PetscFree(dummy);
3544: VecRestoreArrayRead(x,NULL);
3545: return(0);
3546: }
3548: /*@C
3549: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3550: processor's portion of the vector data. You MUST call VecRestoreArray1dRead()
3551: when you no longer need access to the array.
3553: Logically Collective
3555: Input Parameter:
3556: + x - the vector
3557: . m - first dimension of two dimensional array
3558: - mstart - first index you will use in first coordinate direction (often 0)
3560: Output Parameter:
3561: . a - location to put pointer to the array
3563: Level: developer
3565: Notes:
3566: For a vector obtained from DMCreateLocalVector() mstart are likely
3567: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3568: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
3570: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3572: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3573: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3574: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3575: @*/
3576: PetscErrorCode VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3577: {
3579: PetscInt N;
3585: VecGetLocalSize(x,&N);
3586: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3587: VecGetArrayRead(x,(const PetscScalar**)a);
3588: *a -= mstart;
3589: return(0);
3590: }
3592: /*@C
3593: VecRestoreArray1dRead - Restores a vector after VecGetArray1dRead() has been called.
3595: Logically Collective
3597: Input Parameters:
3598: + x - the vector
3599: . m - first dimension of two dimensional array
3600: . mstart - first index you will use in first coordinate direction (often 0)
3601: - a - location of pointer to array obtained from VecGetArray21()
3603: Level: developer
3605: Notes:
3606: For regular PETSc vectors this routine does not involve any copies. For
3607: any special vectors that do not store local vector data in a contiguous
3608: array, this routine will copy the data back into the underlying
3609: vector data structure from the array obtained with VecGetArray1dRead().
3611: This routine actually zeros out the a pointer.
3613: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3614: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3615: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3616: @*/
3617: PetscErrorCode VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3618: {
3624: VecRestoreArrayRead(x,NULL);
3625: return(0);
3626: }
3629: /*@C
3630: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3631: processor's portion of the vector data. You MUST call VecRestoreArray3dRead()
3632: when you no longer need access to the array.
3634: Logically Collective
3636: Input Parameter:
3637: + x - the vector
3638: . m - first dimension of three dimensional array
3639: . n - second dimension of three dimensional array
3640: . p - third dimension of three dimensional array
3641: . mstart - first index you will use in first coordinate direction (often 0)
3642: . nstart - first index in the second coordinate direction (often 0)
3643: - pstart - first index in the third coordinate direction (often 0)
3645: Output Parameter:
3646: . a - location to put pointer to the array
3648: Level: developer
3650: Notes:
3651: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3652: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3653: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3654: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3dRead().
3656: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3658: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3659: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3660: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3661: @*/
3662: PetscErrorCode VecGetArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3663: {
3664: PetscErrorCode ierr;
3665: PetscInt i,N,j;
3666: const PetscScalar *aa;
3667: PetscScalar **b;
3673: VecGetLocalSize(x,&N);
3674: 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);
3675: VecGetArrayRead(x,&aa);
3677: PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3678: b = (PetscScalar**)((*a) + m);
3679: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3680: for (i=0; i<m; i++)
3681: for (j=0; j<n; j++)
3682: b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;
3684: *a -= mstart;
3685: return(0);
3686: }
3688: /*@C
3689: VecRestoreArray3dRead - Restores a vector after VecGetArray3dRead() has been called.
3691: Logically Collective
3693: Input Parameters:
3694: + x - the vector
3695: . m - first dimension of three dimensional array
3696: . n - second dimension of the three dimensional array
3697: . p - third dimension of the three dimensional array
3698: . mstart - first index you will use in first coordinate direction (often 0)
3699: . nstart - first index in the second coordinate direction (often 0)
3700: . pstart - first index in the third coordinate direction (often 0)
3701: - a - location of pointer to array obtained from VecGetArray3dRead()
3703: Level: developer
3705: Notes:
3706: For regular PETSc vectors this routine does not involve any copies. For
3707: any special vectors that do not store local vector data in a contiguous
3708: array, this routine will copy the data back into the underlying
3709: vector data structure from the array obtained with VecGetArray().
3711: This routine actually zeros out the a pointer.
3713: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3714: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3715: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3716: @*/
3717: PetscErrorCode VecRestoreArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3718: {
3720: void *dummy;
3726: dummy = (void*)(*a + mstart);
3727: PetscFree(dummy);
3728: VecRestoreArrayRead(x,NULL);
3729: return(0);
3730: }
3732: /*@C
3733: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3734: processor's portion of the vector data. You MUST call VecRestoreArray4dRead()
3735: when you no longer need access to the array.
3737: Logically Collective
3739: Input Parameter:
3740: + x - the vector
3741: . m - first dimension of four dimensional array
3742: . n - second dimension of four dimensional array
3743: . p - third dimension of four dimensional array
3744: . q - fourth dimension of four dimensional array
3745: . mstart - first index you will use in first coordinate direction (often 0)
3746: . nstart - first index in the second coordinate direction (often 0)
3747: . pstart - first index in the third coordinate direction (often 0)
3748: - qstart - first index in the fourth coordinate direction (often 0)
3750: Output Parameter:
3751: . a - location to put pointer to the array
3753: Level: beginner
3755: Notes:
3756: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3757: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3758: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3759: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
3761: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3763: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3764: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3765: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3766: @*/
3767: PetscErrorCode VecGetArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3768: {
3769: PetscErrorCode ierr;
3770: PetscInt i,N,j,k;
3771: const PetscScalar *aa;
3772: PetscScalar ***b,**c;
3778: VecGetLocalSize(x,&N);
3779: 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);
3780: VecGetArrayRead(x,&aa);
3782: PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3783: b = (PetscScalar***)((*a) + m);
3784: c = (PetscScalar**)(b + m*n);
3785: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3786: for (i=0; i<m; i++)
3787: for (j=0; j<n; j++)
3788: b[i*n+j] = c + i*n*p + j*p - pstart;
3789: for (i=0; i<m; i++)
3790: for (j=0; j<n; j++)
3791: for (k=0; k<p; k++)
3792: c[i*n*p+j*p+k] = (PetscScalar*) aa + i*n*p*q + j*p*q + k*q - qstart;
3793: *a -= mstart;
3794: return(0);
3795: }
3797: /*@C
3798: VecRestoreArray4dRead - Restores a vector after VecGetArray3d() has been called.
3800: Logically Collective
3802: Input Parameters:
3803: + x - the vector
3804: . m - first dimension of four dimensional array
3805: . n - second dimension of the four dimensional array
3806: . p - third dimension of the four dimensional array
3807: . q - fourth dimension of the four dimensional array
3808: . mstart - first index you will use in first coordinate direction (often 0)
3809: . nstart - first index in the second coordinate direction (often 0)
3810: . pstart - first index in the third coordinate direction (often 0)
3811: . qstart - first index in the fourth coordinate direction (often 0)
3812: - a - location of pointer to array obtained from VecGetArray4dRead()
3814: Level: beginner
3816: Notes:
3817: For regular PETSc vectors this routine does not involve any copies. For
3818: any special vectors that do not store local vector data in a contiguous
3819: array, this routine will copy the data back into the underlying
3820: vector data structure from the array obtained with VecGetArray().
3822: This routine actually zeros out the a pointer.
3824: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3825: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3826: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3827: @*/
3828: PetscErrorCode VecRestoreArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3829: {
3831: void *dummy;
3837: dummy = (void*)(*a + mstart);
3838: PetscFree(dummy);
3839: VecRestoreArrayRead(x,NULL);
3840: return(0);
3841: }
3843: #if defined(PETSC_USE_DEBUG)
3845: /*@
3846: VecLockGet - Gets the current lock status of a vector
3848: Logically Collective on Vec
3850: Input Parameter:
3851: . x - the vector
3853: Output Parameter:
3854: . state - greater than zero indicates the vector is locked for read; less then zero indicates the vector is
3855: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3857: Level: beginner
3859: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop()
3860: @*/
3861: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
3862: {
3865: *state = x->lock;
3866: return(0);
3867: }
3869: /*@
3870: VecLockReadPush - Pushes a read-only lock on a vector to prevent it from writing
3872: Logically Collective on Vec
3874: Input Parameter:
3875: . x - the vector
3877: Notes:
3878: If this is set then calls to VecGetArray() or VecSetValues() or any other routines that change the vectors values will fail.
3880: The call can be nested, i.e., called multiple times on the same vector, but each VecLockReadPush(x) has to have one matching
3881: VecLockReadPop(x), which removes the latest read-only lock.
3883: Level: beginner
3885: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPop(), VecLockGet()
3886: @*/
3887: PetscErrorCode VecLockReadPush(Vec x)
3888: {
3891: 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");
3892: x->lock++;
3893: return(0);
3894: }
3896: /*@
3897: VecLockReadPop - Pops a read-only lock from a vector
3899: Logically Collective on Vec
3901: Input Parameter:
3902: . x - the vector
3904: Level: beginner
3906: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockGet()
3907: @*/
3908: PetscErrorCode VecLockReadPop(Vec x)
3909: {
3912: x->lock--;
3913: if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector has been unlocked from read-only access too many times");
3914: return(0);
3915: }
3917: /*@C
3918: VecLockWriteSet_Private - Lock or unlock a vector for exclusive read/write access
3920: Logically Collective on Vec
3922: Input Parameter:
3923: + x - the vector
3924: - flg - PETSC_TRUE to lock the vector for writing; PETSC_FALSE to unlock it.
3926: Notes:
3927: The function is usefull in split-phase computations, which usually have a begin phase and an end phase.
3928: One can call VecLockWriteSet_Private(x,PETSC_TRUE) in the begin phase to lock a vector for exclusive
3929: access, and call VecLockWriteSet_Private(x,PETSC_FALSE) in the end phase to unlock the vector from exclusive
3930: access. In this way, one is ensured no other operations can access the vector in between. The code may like
3933: VecGetArray(x,&xdata); // begin phase
3934: VecLockWriteSet_Private(v,PETSC_TRUE);
3936: Other operations, which can not acceess x anymore (they can access xdata, of course)
3938: VecRestoreArray(x,&vdata); // end phase
3939: VecLockWriteSet_Private(v,PETSC_FALSE);
3941: The call can not be nested on the same vector, in other words, one can not call VecLockWriteSet_Private(x,PETSC_TRUE)
3942: again before calling VecLockWriteSet_Private(v,PETSC_FALSE).
3944: Level: beginner
3946: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop(), VecLockGet()
3947: @*/
3948: PetscErrorCode VecLockWriteSet_Private(Vec x,PetscBool flg)
3949: {
3952: if (flg) {
3953: 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");
3954: 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");
3955: else x->lock = -1;
3956: } else {
3957: 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");
3958: x->lock = 0;
3959: }
3960: return(0);
3961: }
3963: /*@
3964: VecLockPush - Pushes a read-only lock on a vector to prevent it from writing
3966: Level: deprecated
3968: .seealso: VecLockReadPush()
3969: @*/
3970: PetscErrorCode VecLockPush(Vec x)
3971: {
3974: VecLockReadPush(x);
3975: return(0);
3976: }
3978: /*@
3979: VecLockPop - Pops a read-only lock from a vector
3981: Level: deprecated
3983: .seealso: VecLockReadPop()
3984: @*/
3985: PetscErrorCode VecLockPop(Vec x)
3986: {
3989: VecLockReadPop(x);
3990: return(0);
3991: }
3993: #endif