Actual source code: rvector.c
petsc-main 2021-04-16
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: #if defined(PETSC_HAVE_HIP)
13: #include <../src/vec/vec/impls/dvecimpl.h>
14: #include <petsc/private/hipvecimpl.h>
15: #endif
16: static PetscInt VecGetSubVectorSavedStateId = -1;
18: PETSC_EXTERN PetscErrorCode VecValidValues(Vec vec,PetscInt argnum,PetscBool begin)
19: {
20: #if defined(PETSC_USE_DEBUG)
21: PetscErrorCode ierr;
22: PetscInt n,i;
23: const PetscScalar *x;
26: #if defined(PETSC_HAVE_DEVICE)
27: if ((vec->petscnative || vec->ops->getarray) && (vec->offloadmask & PETSC_OFFLOAD_CPU)) {
28: #else
29: if (vec->petscnative || vec->ops->getarray) {
30: #endif
31: VecGetLocalSize(vec,&n);
32: VecGetArrayRead(vec,&x);
33: for (i=0; i<n; i++) {
34: if (begin) {
35: 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);
36: } else {
37: 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);
38: }
39: }
40: VecRestoreArrayRead(vec,&x);
41: }
42: #else
44: #endif
45: return(0);
46: }
48: /*@
49: VecMaxPointwiseDivide - Computes the maximum of the componentwise division max = max_i abs(x_i/y_i).
51: Logically Collective on Vec
53: Input Parameters:
54: . x, y - the vectors
56: Output Parameter:
57: . max - the result
59: Level: advanced
61: Notes:
62: x and y may be the same vector
63: if a particular y_i is zero, it is treated as 1 in the above formula
65: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs()
66: @*/
67: PetscErrorCode VecMaxPointwiseDivide(Vec x,Vec y,PetscReal *max)
68: {
78: VecCheckSameSize(x,1,y,2);
79: (*x->ops->maxpointwisedivide)(x,y,max);
80: return(0);
81: }
83: /*@
84: VecDot - Computes the vector dot product.
86: Collective on Vec
88: Input Parameters:
89: . x, y - the vectors
91: Output Parameter:
92: . val - the dot product
94: Performance Issues:
95: $ per-processor memory bandwidth
96: $ interprocessor latency
97: $ work load inbalance that causes certain processes to arrive much earlier than others
99: Notes for Users of Complex Numbers:
100: For complex vectors, VecDot() computes
101: $ val = (x,y) = y^H x,
102: where y^H denotes the conjugate transpose of y. Note that this corresponds to the usual "mathematicians" complex
103: inner product where the SECOND argument gets the complex conjugate. Since the BLASdot() complex conjugates the first
104: first argument we call the BLASdot() with the arguments reversed.
106: Use VecTDot() for the indefinite form
107: $ val = (x,y) = y^T x,
108: where y^T denotes the transpose of y.
110: Level: intermediate
113: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDotRealPart()
114: @*/
115: PetscErrorCode VecDot(Vec x,Vec y,PetscScalar *val)
116: {
126: VecCheckSameSize(x,1,y,2);
128: PetscLogEventBegin(VEC_Dot,x,y,0,0);
129: (*x->ops->dot)(x,y,val);
130: PetscLogEventEnd(VEC_Dot,x,y,0,0);
131: return(0);
132: }
134: /*@
135: VecDotRealPart - Computes the real part of the vector dot product.
137: Collective on Vec
139: Input Parameters:
140: . x, y - the vectors
142: Output Parameter:
143: . val - the real part of the dot product;
145: Performance Issues:
146: $ per-processor memory bandwidth
147: $ interprocessor latency
148: $ work load inbalance that causes certain processes to arrive much earlier than others
150: Notes for Users of Complex Numbers:
151: See VecDot() for more details on the definition of the dot product for complex numbers
153: For real numbers this returns the same value as VecDot()
155: 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
156: the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)
158: Developer Note: This is not currently optimized to compute only the real part of the dot product.
160: Level: intermediate
163: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDot(), VecDotNorm2()
164: @*/
165: PetscErrorCode VecDotRealPart(Vec x,Vec y,PetscReal *val)
166: {
168: PetscScalar fdot;
171: VecDot(x,y,&fdot);
172: *val = PetscRealPart(fdot);
173: return(0);
174: }
176: /*@
177: VecNorm - Computes the vector norm.
179: Collective on Vec
181: Input Parameters:
182: + x - the vector
183: - type - one of NORM_1, NORM_2, NORM_INFINITY. Also available
184: NORM_1_AND_2, which computes both norms and stores them
185: in a two element array.
187: Output Parameter:
188: . val - the norm
190: Notes:
191: $ NORM_1 denotes sum_i |x_i|
192: $ NORM_2 denotes sqrt(sum_i |x_i|^2)
193: $ NORM_INFINITY denotes max_i |x_i|
195: For complex numbers NORM_1 will return the traditional 1 norm of the 2 norm of the complex numbers; that is the 1
196: norm of the absolute values of the complex entries. In PETSc 3.6 and earlier releases it returned the 1 norm of
197: the 1 norm of the complex entries (what is returned by the BLAS routine asum()). Both are valid norms but most
198: people expect the former.
200: Level: intermediate
202: Performance Issues:
203: $ per-processor memory bandwidth
204: $ interprocessor latency
205: $ work load inbalance that causes certain processes to arrive much earlier than others
208: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNormAvailable(),
209: VecNormBegin(), VecNormEnd()
211: @*/
213: PetscErrorCode VecNorm(Vec x,NormType type,PetscReal *val)
214: {
215: PetscBool flg;
223: /*
224: * Cached data?
225: */
226: if (type!=NORM_1_AND_2) {
227: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,flg);
228: if (flg) return(0);
229: }
230: PetscLogEventBegin(VEC_Norm,x,0,0,0);
231: (*x->ops->norm)(x,type,val);
232: PetscLogEventEnd(VEC_Norm,x,0,0,0);
233: if (type!=NORM_1_AND_2) {
234: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[type],*val);
235: }
236: return(0);
237: }
239: /*@
240: VecNormAvailable - Returns the vector norm if it is already known.
242: Not Collective
244: Input Parameters:
245: + x - the vector
246: - type - one of NORM_1, NORM_2, NORM_INFINITY. Also available
247: NORM_1_AND_2, which computes both norms and stores them
248: in a two element array.
250: Output Parameter:
251: + available - PETSC_TRUE if the val returned is valid
252: - val - the norm
254: Notes:
255: $ NORM_1 denotes sum_i |x_i|
256: $ NORM_2 denotes sqrt(sum_i (x_i)^2)
257: $ NORM_INFINITY denotes max_i |x_i|
259: Level: intermediate
261: Performance Issues:
262: $ per-processor memory bandwidth
263: $ interprocessor latency
264: $ work load inbalance that causes certain processes to arrive much earlier than others
266: Compile Option:
267: PETSC_HAVE_SLOW_BLAS_NORM2 will cause a C (loop unrolled) version of the norm to be used, rather
268: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
269: (as opposed to vendor provided) because the FORTRAN BLAS NRM2() routine is very slow.
272: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNorm()
273: VecNormBegin(), VecNormEnd()
275: @*/
276: PetscErrorCode VecNormAvailable(Vec x,NormType type,PetscBool *available,PetscReal *val)
277: {
285: *available = PETSC_FALSE;
286: if (type!=NORM_1_AND_2) {
287: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,*available);
288: }
289: return(0);
290: }
292: /*@
293: VecNormalize - Normalizes a vector by 2-norm.
295: Collective on Vec
297: Input Parameters:
298: + x - the vector
300: Output Parameter:
301: . x - the normalized vector
302: - val - the vector norm before normalization
304: Level: intermediate
307: @*/
308: PetscErrorCode VecNormalize(Vec x,PetscReal *val)
309: {
311: PetscReal norm;
316: PetscLogEventBegin(VEC_Normalize,x,0,0,0);
317: VecNorm(x,NORM_2,&norm);
318: if (norm == 0.0) {
319: PetscInfo(x,"Vector of zero norm can not be normalized; Returning only the zero norm\n");
320: } else if (norm != 1.0) {
321: PetscScalar tmp = 1.0/norm;
322: VecScale(x,tmp);
323: }
324: if (val) *val = norm;
325: PetscLogEventEnd(VEC_Normalize,x,0,0,0);
326: return(0);
327: }
329: /*@C
330: VecMax - Determines the vector component with maximum real part and its location.
332: Collective on Vec
334: Input Parameter:
335: . x - the vector
337: Output Parameters:
338: + p - the location of val (pass NULL if you don't want this)
339: - val - the maximum component
341: Notes:
342: Returns the value PETSC_MIN_REAL and p = -1 if the vector is of length 0.
344: Returns the smallest index with the maximum value
345: Level: intermediate
348: .seealso: VecNorm(), VecMin()
349: @*/
350: PetscErrorCode VecMax(Vec x,PetscInt *p,PetscReal *val)
351: {
358: PetscLogEventBegin(VEC_Max,x,0,0,0);
359: (*x->ops->max)(x,p,val);
360: PetscLogEventEnd(VEC_Max,x,0,0,0);
361: return(0);
362: }
364: /*@C
365: VecMin - Determines the vector component with minimum real part and its location.
367: Collective on Vec
369: Input Parameters:
370: . x - the vector
372: Output Parameter:
373: + p - the location of val (pass NULL if you don't want this location)
374: - val - the minimum component
376: Level: intermediate
378: Notes:
379: Returns the value PETSC_MAX_REAL and p = -1 if the vector is of length 0.
381: This returns the smallest index with the minumum value
384: .seealso: VecMax()
385: @*/
386: PetscErrorCode VecMin(Vec x,PetscInt *p,PetscReal *val)
387: {
394: PetscLogEventBegin(VEC_Min,x,0,0,0);
395: (*x->ops->min)(x,p,val);
396: PetscLogEventEnd(VEC_Min,x,0,0,0);
397: return(0);
398: }
400: /*@
401: VecTDot - Computes an indefinite vector dot product. That is, this
402: routine does NOT use the complex conjugate.
404: Collective on Vec
406: Input Parameters:
407: . x, y - the vectors
409: Output Parameter:
410: . val - the dot product
412: Notes for Users of Complex Numbers:
413: For complex vectors, VecTDot() computes the indefinite form
414: $ val = (x,y) = y^T x,
415: where y^T denotes the transpose of y.
417: Use VecDot() for the inner product
418: $ val = (x,y) = y^H x,
419: where y^H denotes the conjugate transpose of y.
421: Level: intermediate
423: .seealso: VecDot(), VecMTDot()
424: @*/
425: PetscErrorCode VecTDot(Vec x,Vec y,PetscScalar *val)
426: {
436: VecCheckSameSize(x,1,y,2);
438: PetscLogEventBegin(VEC_TDot,x,y,0,0);
439: (*x->ops->tdot)(x,y,val);
440: PetscLogEventEnd(VEC_TDot,x,y,0,0);
441: return(0);
442: }
444: /*@
445: VecScale - Scales a vector.
447: Not collective on Vec
449: Input Parameters:
450: + x - the vector
451: - alpha - the scalar
453: Output Parameter:
454: . x - the scaled vector
456: Note:
457: For a vector with n components, VecScale() computes
458: $ x[i] = alpha * x[i], for i=1,...,n.
460: Level: intermediate
463: @*/
464: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
465: {
466: PetscReal norms[4] = {0.0,0.0,0.0, 0.0};
467: PetscBool flgs[4];
469: PetscInt i;
474: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
475: PetscLogEventBegin(VEC_Scale,x,0,0,0);
476: if (alpha != (PetscScalar)1.0) {
477: VecSetErrorIfLocked(x,1);
478: /* get current stashed norms */
479: for (i=0; i<4; i++) {
480: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
481: }
482: (*x->ops->scale)(x,alpha);
483: PetscObjectStateIncrease((PetscObject)x);
484: /* put the scaled stashed norms back into the Vec */
485: for (i=0; i<4; i++) {
486: if (flgs[i]) {
487: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],PetscAbsScalar(alpha)*norms[i]);
488: }
489: }
490: }
491: PetscLogEventEnd(VEC_Scale,x,0,0,0);
492: return(0);
493: }
495: /*@
496: VecSet - Sets all components of a vector to a single scalar value.
498: Logically Collective on Vec
500: Input Parameters:
501: + x - the vector
502: - alpha - the scalar
504: Output Parameter:
505: . x - the vector
507: Note:
508: For a vector of dimension n, VecSet() computes
509: $ x[i] = alpha, for i=1,...,n,
510: so that all vector entries then equal the identical
511: scalar value, alpha. Use the more general routine
512: VecSetValues() to set different vector entries.
514: You CANNOT call this after you have called VecSetValues() but before you call
515: VecAssemblyBegin/End().
517: Level: beginner
519: .seealso VecSetValues(), VecSetValuesBlocked(), VecSetRandom()
521: @*/
522: PetscErrorCode VecSet(Vec x,PetscScalar alpha)
523: {
524: PetscReal val;
530: 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()");
532: VecSetErrorIfLocked(x,1);
534: PetscLogEventBegin(VEC_Set,x,0,0,0);
535: (*x->ops->set)(x,alpha);
536: PetscLogEventEnd(VEC_Set,x,0,0,0);
537: PetscObjectStateIncrease((PetscObject)x);
539: /* norms can be simply set (if |alpha|*N not too large) */
540: val = PetscAbsScalar(alpha);
541: if (x->map->N == 0) {
542: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],0.0l);
543: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],0.0);
544: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],0.0);
545: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],0.0);
546: } else if (val > PETSC_MAX_REAL/x->map->N) {
547: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
548: } else {
549: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],x->map->N * val);
550: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
551: val = PetscSqrtReal((PetscReal)x->map->N) * val;
552: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],val);
553: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],val);
554: }
555: return(0);
556: }
559: /*@
560: VecAXPY - Computes y = alpha x + y.
562: Logically Collective on Vec
564: Input Parameters:
565: + alpha - the scalar
566: - x, y - the vectors
568: Output Parameter:
569: . y - output vector
571: Level: intermediate
573: Notes:
574: x and y MUST be different vectors
575: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
577: $ VecAXPY(y,alpha,x) y = alpha x + y
578: $ VecAYPX(y,beta,x) y = x + beta y
579: $ VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
580: $ VecWAXPY(w,alpha,x,y) w = alpha x + y
581: $ VecAXPBYPCZ(w,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
582: $ VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
585: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPBYPCZ(), VecAXPBY()
586: @*/
587: PetscErrorCode VecAXPY(Vec y,PetscScalar alpha,Vec x)
588: {
597: VecCheckSameSize(x,1,y,3);
598: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
600: if (alpha == (PetscScalar)0.0) return(0);
601: VecSetErrorIfLocked(y,1);
603: VecLockReadPush(x);
604: PetscLogEventBegin(VEC_AXPY,x,y,0,0);
605: (*y->ops->axpy)(y,alpha,x);
606: PetscLogEventEnd(VEC_AXPY,x,y,0,0);
607: VecLockReadPop(x);
608: PetscObjectStateIncrease((PetscObject)y);
609: return(0);
610: }
612: /*@
613: VecAXPBY - Computes y = alpha x + beta y.
615: Logically Collective on Vec
617: Input Parameters:
618: + alpha,beta - the scalars
619: - x, y - the vectors
621: Output Parameter:
622: . y - output vector
624: Level: intermediate
626: Notes:
627: x and y MUST be different vectors
628: The implementation is optimized for alpha and/or beta values of 0.0 and 1.0
631: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ()
632: @*/
633: PetscErrorCode VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
634: {
643: VecCheckSameSize(y,1,x,4);
644: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
647: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) return(0);
648: VecSetErrorIfLocked(y,1);
649: PetscLogEventBegin(VEC_AXPY,x,y,0,0);
650: (*y->ops->axpby)(y,alpha,beta,x);
651: PetscLogEventEnd(VEC_AXPY,x,y,0,0);
652: PetscObjectStateIncrease((PetscObject)y);
653: return(0);
654: }
656: /*@
657: VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z
659: Logically Collective on Vec
661: Input Parameters:
662: + alpha,beta, gamma - the scalars
663: - x, y, z - the vectors
665: Output Parameter:
666: . z - output vector
668: Level: intermediate
670: Notes:
671: x, y and z must be different vectors
672: The implementation is optimized for alpha of 1.0 and gamma of 1.0 or 0.0
675: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBY()
676: @*/
677: PetscErrorCode VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
678: {
690: VecCheckSameSize(x,1,y,5);
691: VecCheckSameSize(x,1,z,6);
692: if (x == y || x == z) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
693: if (y == z) SETERRQ(PetscObjectComm((PetscObject)y),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
697: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) return(0);
698: VecSetErrorIfLocked(z,1);
700: PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
701: (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
702: PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
703: PetscObjectStateIncrease((PetscObject)z);
704: return(0);
705: }
707: /*@
708: VecAYPX - Computes y = x + beta y.
710: Logically Collective on Vec
712: Input Parameters:
713: + beta - the scalar
714: - x, y - the vectors
716: Output Parameter:
717: . y - output vector
719: Level: intermediate
721: Notes:
722: x and y MUST be different vectors
723: The implementation is optimized for beta of -1.0, 0.0, and 1.0
726: .seealso: VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
727: @*/
728: PetscErrorCode VecAYPX(Vec y,PetscScalar beta,Vec x)
729: {
738: VecCheckSameSize(x,1,y,3);
739: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
741: VecSetErrorIfLocked(y,1);
743: PetscLogEventBegin(VEC_AYPX,x,y,0,0);
744: (*y->ops->aypx)(y,beta,x);
745: PetscLogEventEnd(VEC_AYPX,x,y,0,0);
746: PetscObjectStateIncrease((PetscObject)y);
747: return(0);
748: }
751: /*@
752: VecWAXPY - Computes w = alpha x + y.
754: Logically Collective on Vec
756: Input Parameters:
757: + alpha - the scalar
758: - x, y - the vectors
760: Output Parameter:
761: . w - the result
763: Level: intermediate
765: Notes:
766: w cannot be either x or y, but x and y can be the same
767: The implementation is optimzed for alpha of -1.0, 0.0, and 1.0
770: .seealso: VecAXPY(), VecAYPX(), VecAXPBY(), VecMAXPY(), VecAXPBYPCZ()
771: @*/
772: PetscErrorCode VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
773: {
785: VecCheckSameSize(x,3,y,4);
786: VecCheckSameSize(x,3,w,1);
787: if (w == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector y, suggest VecAXPY()");
788: if (w == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector x, suggest VecAYPX()");
790: VecSetErrorIfLocked(w,1);
792: PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
793: (*w->ops->waxpy)(w,alpha,x,y);
794: PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
795: PetscObjectStateIncrease((PetscObject)w);
796: return(0);
797: }
800: /*@C
801: VecSetValues - Inserts or adds values into certain locations of a vector.
803: Not Collective
805: Input Parameters:
806: + x - vector to insert in
807: . ni - number of elements to add
808: . ix - indices where to add
809: . y - array of values
810: - iora - either INSERT_VALUES or ADD_VALUES, where
811: ADD_VALUES adds values to any existing entries, and
812: INSERT_VALUES replaces existing entries with new values
814: Notes:
815: VecSetValues() sets x[ix[i]] = y[i], for i=0,...,ni-1.
817: Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
818: options cannot be mixed without intervening calls to the assembly
819: routines.
821: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
822: MUST be called after all calls to VecSetValues() have been completed.
824: VecSetValues() uses 0-based indices in Fortran as well as in C.
826: If you call VecSetOption(x, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE),
827: negative indices may be passed in ix. These rows are
828: simply ignored. This allows easily inserting element load matrices
829: with homogeneous Dirchlet boundary conditions that you don't want represented
830: in the vector.
832: Level: beginner
834: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
835: VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
836: @*/
837: PetscErrorCode VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
838: {
843: if (!ni) return(0);
848: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
849: (*x->ops->setvalues)(x,ni,ix,y,iora);
850: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
851: PetscObjectStateIncrease((PetscObject)x);
852: return(0);
853: }
855: /*@C
856: VecGetValues - Gets values from certain locations of a vector. Currently
857: can only get values on the same processor
859: Not Collective
861: Input Parameters:
862: + x - vector to get values from
863: . ni - number of elements to get
864: - ix - indices where to get them from (in global 1d numbering)
866: Output Parameter:
867: . y - array of values
869: Notes:
870: The user provides the allocated array y; it is NOT allocated in this routine
872: VecGetValues() gets y[i] = x[ix[i]], for i=0,...,ni-1.
874: VecAssemblyBegin() and VecAssemblyEnd() MUST be called before calling this
876: VecGetValues() uses 0-based indices in Fortran as well as in C.
878: If you call VecSetOption(x, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE),
879: negative indices may be passed in ix. These rows are
880: simply ignored.
882: Level: beginner
884: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues()
885: @*/
886: PetscErrorCode VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
887: {
892: if (!ni) return(0);
896: (*x->ops->getvalues)(x,ni,ix,y);
897: return(0);
898: }
900: /*@C
901: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
903: Not Collective
905: Input Parameters:
906: + x - vector to insert in
907: . ni - number of blocks to add
908: . ix - indices where to add in block count, rather than element count
909: . y - array of values
910: - iora - either INSERT_VALUES or ADD_VALUES, where
911: ADD_VALUES adds values to any existing entries, and
912: INSERT_VALUES replaces existing entries with new values
914: Notes:
915: VecSetValuesBlocked() sets x[bs*ix[i]+j] = y[bs*i+j],
916: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
918: Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
919: options cannot be mixed without intervening calls to the assembly
920: routines.
922: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
923: MUST be called after all calls to VecSetValuesBlocked() have been completed.
925: VecSetValuesBlocked() uses 0-based indices in Fortran as well as in C.
927: Negative indices may be passed in ix, these rows are
928: simply ignored. This allows easily inserting element load matrices
929: with homogeneous Dirchlet boundary conditions that you don't want represented
930: in the vector.
932: Level: intermediate
934: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
935: VecSetValues()
936: @*/
937: PetscErrorCode VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
938: {
943: if (!ni) return(0);
948: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
949: (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
950: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
951: PetscObjectStateIncrease((PetscObject)x);
952: return(0);
953: }
956: /*@C
957: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
958: using a local ordering of the nodes.
960: Not Collective
962: Input Parameters:
963: + x - vector to insert in
964: . ni - number of elements to add
965: . ix - indices where to add
966: . y - array of values
967: - iora - either INSERT_VALUES or ADD_VALUES, where
968: ADD_VALUES adds values to any existing entries, and
969: INSERT_VALUES replaces existing entries with new values
971: Level: intermediate
973: Notes:
974: VecSetValuesLocal() sets x[ix[i]] = y[i], for i=0,...,ni-1.
976: Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
977: options cannot be mixed without intervening calls to the assembly
978: routines.
980: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
981: MUST be called after all calls to VecSetValuesLocal() have been completed.
983: VecSetValuesLocal() uses 0-based indices in Fortran as well as in C.
985: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetLocalToGlobalMapping(),
986: VecSetValuesBlockedLocal()
987: @*/
988: PetscErrorCode VecSetValuesLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
989: {
991: PetscInt lixp[128],*lix = lixp;
995: if (!ni) return(0);
1000: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1001: if (!x->ops->setvalueslocal) {
1002: if (!x->map->mapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMapping()");
1003: if (ni > 128) {
1004: PetscMalloc1(ni,&lix);
1005: }
1006: ISLocalToGlobalMappingApply(x->map->mapping,ni,(PetscInt*)ix,lix);
1007: (*x->ops->setvalues)(x,ni,lix,y,iora);
1008: if (ni > 128) {
1009: PetscFree(lix);
1010: }
1011: } else {
1012: (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
1013: }
1014: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1015: PetscObjectStateIncrease((PetscObject)x);
1016: return(0);
1017: }
1019: /*@
1020: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1021: using a local ordering of the nodes.
1023: Not Collective
1025: Input Parameters:
1026: + x - vector to insert in
1027: . ni - number of blocks to add
1028: . ix - indices where to add in block count, not element count
1029: . y - array of values
1030: - iora - either INSERT_VALUES or ADD_VALUES, where
1031: ADD_VALUES adds values to any existing entries, and
1032: INSERT_VALUES replaces existing entries with new values
1034: Level: intermediate
1036: Notes:
1037: VecSetValuesBlockedLocal() sets x[bs*ix[i]+j] = y[bs*i+j],
1038: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with VecSetBlockSize().
1040: Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1041: options cannot be mixed without intervening calls to the assembly
1042: routines.
1044: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
1045: MUST be called after all calls to VecSetValuesBlockedLocal() have been completed.
1047: VecSetValuesBlockedLocal() uses 0-based indices in Fortran as well as in C.
1050: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
1051: VecSetLocalToGlobalMapping()
1052: @*/
1053: PetscErrorCode VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1054: {
1056: PetscInt lixp[128],*lix = lixp;
1060: if (!ni) return(0);
1064: if (ni > 128) {
1065: PetscMalloc1(ni,&lix);
1066: }
1068: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1069: ISLocalToGlobalMappingApplyBlock(x->map->mapping,ni,(PetscInt*)ix,lix);
1070: (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1071: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1072: if (ni > 128) {
1073: PetscFree(lix);
1074: }
1075: PetscObjectStateIncrease((PetscObject)x);
1076: return(0);
1077: }
1079: /*@
1080: VecMTDot - Computes indefinite vector multiple dot products.
1081: That is, it does NOT use the complex conjugate.
1083: Collective on Vec
1085: Input Parameters:
1086: + x - one vector
1087: . nv - number of vectors
1088: - y - array of vectors. Note that vectors are pointers
1090: Output Parameter:
1091: . val - array of the dot products
1093: Notes for Users of Complex Numbers:
1094: For complex vectors, VecMTDot() computes the indefinite form
1095: $ val = (x,y) = y^T x,
1096: where y^T denotes the transpose of y.
1098: Use VecMDot() for the inner product
1099: $ val = (x,y) = y^H x,
1100: where y^H denotes the conjugate transpose of y.
1102: Level: intermediate
1105: .seealso: VecMDot(), VecTDot()
1106: @*/
1107: PetscErrorCode VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1108: {
1114: if (!nv) return(0);
1121: VecCheckSameSize(x,1,*y,3);
1123: PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1124: (*x->ops->mtdot)(x,nv,y,val);
1125: PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1126: return(0);
1127: }
1129: /*@
1130: VecMDot - Computes vector multiple dot products.
1132: Collective on Vec
1134: Input Parameters:
1135: + x - one vector
1136: . nv - number of vectors
1137: - y - array of vectors.
1139: Output Parameter:
1140: . val - array of the dot products (does not allocate the array)
1142: Notes for Users of Complex Numbers:
1143: For complex vectors, VecMDot() computes
1144: $ val = (x,y) = y^H x,
1145: where y^H denotes the conjugate transpose of y.
1147: Use VecMTDot() for the indefinite form
1148: $ val = (x,y) = y^T x,
1149: where y^T denotes the transpose of y.
1151: Level: intermediate
1154: .seealso: VecMTDot(), VecDot()
1155: @*/
1156: PetscErrorCode VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1157: {
1163: if (!nv) return(0);
1164: if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1171: VecCheckSameSize(x,1,*y,3);
1173: PetscLogEventBegin(VEC_MDot,x,*y,0,0);
1174: (*x->ops->mdot)(x,nv,y,val);
1175: PetscLogEventEnd(VEC_MDot,x,*y,0,0);
1176: return(0);
1177: }
1179: /*@
1180: VecMAXPY - Computes y = y + sum alpha[i] x[i]
1182: Logically Collective on Vec
1184: Input Parameters:
1185: + nv - number of scalars and x-vectors
1186: . alpha - array of scalars
1187: . y - one vector
1188: - x - array of vectors
1190: Level: intermediate
1192: Notes:
1193: y cannot be any of the x vectors
1195: .seealso: VecAYPX(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
1196: @*/
1197: PetscErrorCode VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1198: {
1200: PetscInt i;
1201: PetscBool nonzero;
1206: if (!nv) return(0);
1207: if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1214: VecCheckSameSize(y,1,*x,4);
1216: for (i=0, nonzero = PETSC_FALSE; i<nv && !nonzero; i++) nonzero = (PetscBool)(nonzero || alpha[i] != (PetscScalar)0.0);
1217: if (!nonzero) return(0);
1218: VecSetErrorIfLocked(y,1);
1219: PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1220: (*y->ops->maxpy)(y,nv,alpha,x);
1221: PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1222: PetscObjectStateIncrease((PetscObject)y);
1223: return(0);
1224: }
1226: /*@
1227: VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
1228: in the order they appear in the array. The concatenated vector resides on the same
1229: communicator and is the same type as the source vectors.
1231: Collective on X
1233: Input Arguments:
1234: + nx - number of vectors to be concatenated
1235: - X - array containing the vectors to be concatenated in the order of concatenation
1237: Output Arguments:
1238: + Y - concatenated vector
1239: - x_is - array of index sets corresponding to the concatenated components of Y (NULL if not needed)
1241: Notes:
1242: Concatenation is similar to the functionality of a VecNest object; they both represent combination of
1243: different vector spaces. However, concatenated vectors do not store any information about their
1244: sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1245: manipulation of data in the concatenated vector that corresponds to the original components at creation.
1247: This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1248: has to operate on combined vector spaces and cannot utilize VecNest objects due to incompatibility with
1249: bound projections.
1251: Level: advanced
1253: .seealso: VECNEST, VECSCATTER, VecScatterCreate()
1254: @*/
1255: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1256: {
1257: MPI_Comm comm;
1258: VecType vec_type;
1259: Vec Ytmp, Xtmp;
1260: IS *is_tmp;
1261: PetscInt i, shift=0, Xnl, Xng, Xbegin;
1270: if ((*X)->ops->concatenate) {
1271: /* use the dedicated concatenation function if available */
1272: (*(*X)->ops->concatenate)(nx,X,Y,x_is);
1273: } else {
1274: /* loop over vectors and start creating IS */
1275: comm = PetscObjectComm((PetscObject)(*X));
1276: VecGetType(*X, &vec_type);
1277: PetscMalloc1(nx, &is_tmp);
1278: for (i=0; i<nx; i++) {
1279: VecGetSize(X[i], &Xng);
1280: VecGetLocalSize(X[i], &Xnl);
1281: VecGetOwnershipRange(X[i], &Xbegin, NULL);
1282: ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]);
1283: shift += Xng;
1284: }
1285: /* create the concatenated vector */
1286: VecCreate(comm, &Ytmp);
1287: VecSetType(Ytmp, vec_type);
1288: VecSetSizes(Ytmp, PETSC_DECIDE, shift);
1289: VecSetUp(Ytmp);
1290: /* copy data from X array to Y and return */
1291: for (i=0; i<nx; i++) {
1292: VecGetSubVector(Ytmp, is_tmp[i], &Xtmp);
1293: VecCopy(X[i], Xtmp);
1294: VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp);
1295: }
1296: *Y = Ytmp;
1297: if (x_is) {
1298: *x_is = is_tmp;
1299: } else {
1300: for (i=0; i<nx; i++) {
1301: ISDestroy(&is_tmp[i]);
1302: }
1303: PetscFree(is_tmp);
1304: }
1305: }
1306: return(0);
1307: }
1309: /*@
1310: VecGetSubVector - Gets a vector representing part of another vector
1312: Collective on X and IS
1314: Input Arguments:
1315: + X - vector from which to extract a subvector
1316: - is - index set representing portion of X to extract
1318: Output Arguments:
1319: . Y - subvector corresponding to is
1321: Level: advanced
1323: Notes:
1324: The subvector Y should be returned with VecRestoreSubVector().
1325: X and is must be defined on the same communicator
1327: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1328: modifying the subvector. Other non-overlapping subvectors can still be obtained from X using this function.
1329: 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.
1331: .seealso: MatCreateSubMatrix()
1332: @*/
1333: PetscErrorCode VecGetSubVector(Vec X,IS is,Vec *Y)
1334: {
1335: PetscErrorCode ierr;
1336: Vec Z;
1343: if (X->ops->getsubvector) {
1344: (*X->ops->getsubvector)(X,is,&Z);
1345: } else { /* Default implementation currently does no caching */
1346: PetscInt gstart,gend,start;
1347: PetscBool red[2] = { PETSC_TRUE, PETSC_TRUE };
1348: PetscInt n,N,ibs,vbs,bs = -1;
1350: ISGetLocalSize(is,&n);
1351: ISGetSize(is,&N);
1352: ISGetBlockSize(is,&ibs);
1353: VecGetBlockSize(X,&vbs);
1354: VecGetOwnershipRange(X,&gstart,&gend);
1355: ISContiguousLocal(is,gstart,gend,&start,&red[0]);
1356: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1357: if (ibs > 1) {
1358: MPIU_Allreduce(MPI_IN_PLACE,red,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1359: bs = ibs;
1360: } else {
1361: if (n%vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1362: MPIU_Allreduce(MPI_IN_PLACE,red,2,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1363: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1364: }
1365: if (red[0]) { /* We can do a no-copy implementation */
1366: const PetscScalar *x;
1367: PetscInt state = 0;
1368: PetscBool isstd,iscuda,iship;
1370: PetscObjectTypeCompareAny((PetscObject)X,&isstd,VECSEQ,VECMPI,VECSTANDARD,"");
1371: PetscObjectTypeCompareAny((PetscObject)X,&iscuda,VECSEQCUDA,VECMPICUDA,"");
1372: PetscObjectTypeCompareAny((PetscObject)X,&iship,VECSEQHIP,VECMPIHIP,"");
1373: if (iscuda) {
1374: #if defined(PETSC_HAVE_CUDA)
1375: const PetscScalar *x_d;
1376: PetscMPIInt size;
1377: PetscOffloadMask flg;
1379: VecCUDAGetArrays_Private(X,&x,&x_d,&flg);
1380: if (flg == PETSC_OFFLOAD_UNALLOCATED) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for PETSC_OFFLOAD_UNALLOCATED");
1381: if (n && !x && !x_d) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Missing vector data");
1382: if (x) x += start;
1383: if (x_d) x_d += start;
1384: MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1385: if (size == 1) {
1386: VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X),bs,n,x,x_d,&Z);
1387: } else {
1388: VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X),bs,n,N,x,x_d,&Z);
1389: }
1390: Z->offloadmask = flg;
1391: #endif
1392: } else if (iship) {
1393: #if defined(PETSC_HAVE_HIP)
1394: const PetscScalar *x_d;
1395: PetscMPIInt size;
1396: PetscOffloadMask flg;
1398: VecHIPGetArrays_Private(X,&x,&x_d,&flg);
1399: if (flg == PETSC_OFFLOAD_UNALLOCATED) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for PETSC_OFFLOAD_UNALLOCATED");
1400: if (n && !x && !x_d) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Missing vector data");
1401: if (x) x += start;
1402: if (x_d) x_d += start;
1403: MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1404: if (size == 1) {
1405: VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X),bs,n,x,x_d,&Z);
1406: } else {
1407: VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X),bs,n,N,x,x_d,&Z);
1408: }
1409: Z->offloadmask = flg;
1410: #endif
1411: } else if (isstd) {
1412: PetscMPIInt size;
1414: MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1415: VecGetArrayRead(X,&x);
1416: if (x) x += start;
1417: if (size == 1) {
1418: VecCreateSeqWithArray(PetscObjectComm((PetscObject)X),bs,n,x,&Z);
1419: } else {
1420: VecCreateMPIWithArray(PetscObjectComm((PetscObject)X),bs,n,N,x,&Z);
1421: }
1422: VecRestoreArrayRead(X,&x);
1423: } else { /* default implementation: use place array */
1424: VecGetArrayRead(X,&x);
1425: VecCreate(PetscObjectComm((PetscObject)X),&Z);
1426: VecSetType(Z,((PetscObject)X)->type_name);
1427: VecSetSizes(Z,n,N);
1428: VecSetBlockSize(Z,bs);
1429: VecPlaceArray(Z,x ? x+start : NULL);
1430: VecRestoreArrayRead(X,&x);
1431: }
1433: /* this is relevant only in debug mode */
1434: VecLockGet(X,&state);
1435: if (state) {
1436: VecLockReadPush(Z);
1437: }
1438: Z->ops->placearray = NULL;
1439: Z->ops->replacearray = NULL;
1440: } else { /* Have to create a scatter and do a copy */
1441: VecScatter scatter;
1443: VecCreate(PetscObjectComm((PetscObject)is),&Z);
1444: VecSetSizes(Z,n,N);
1445: VecSetBlockSize(Z,bs);
1446: VecSetType(Z,((PetscObject)X)->type_name);
1447: VecScatterCreate(X,is,Z,NULL,&scatter);
1448: VecScatterBegin(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1449: VecScatterEnd(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1450: PetscObjectCompose((PetscObject)Z,"VecGetSubVector_Scatter",(PetscObject)scatter);
1451: VecScatterDestroy(&scatter);
1452: }
1453: }
1454: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1455: if (VecGetSubVectorSavedStateId < 0) {PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);}
1456: PetscObjectComposedDataSetInt((PetscObject)Z,VecGetSubVectorSavedStateId,1);
1457: *Y = Z;
1458: return(0);
1459: }
1461: /*@
1462: VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()
1464: Collective on IS
1466: Input Arguments:
1467: + X - vector from which subvector was obtained
1468: . is - index set representing the subset of X
1469: - Y - subvector being restored
1471: Level: advanced
1473: .seealso: VecGetSubVector()
1474: @*/
1475: PetscErrorCode VecRestoreSubVector(Vec X,IS is,Vec *Y)
1476: {
1485: if (X->ops->restoresubvector) {
1486: (*X->ops->restoresubvector)(X,is,Y);
1487: } else {
1488: PETSC_UNUSED PetscObjectState dummystate = 0;
1489: PetscBool valid;
1491: PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,dummystate,valid);
1492: if (!valid) {
1493: VecScatter scatter;
1494: PetscInt state;
1496: VecLockGet(X,&state);
1497: if (state != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vec X is locked for read-only or read/write access");
1499: PetscObjectQuery((PetscObject)*Y,"VecGetSubVector_Scatter",(PetscObject*)&scatter);
1500: if (scatter) {
1501: VecScatterBegin(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1502: VecScatterEnd(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1503: } else {
1504: PetscBool iscuda,iship;
1505: PetscObjectTypeCompareAny((PetscObject)X,&iscuda,VECSEQCUDA,VECMPICUDA,"");
1506: PetscObjectTypeCompareAny((PetscObject)X,&iship,VECSEQHIP,VECMPIHIP,"");
1508: if (iscuda) {
1509: #if defined(PETSC_HAVE_CUDA)
1510: PetscOffloadMask ymask = (*Y)->offloadmask;
1512: /* The offloadmask of X dictates where to move memory
1513: If X GPU data is valid, then move Y data on GPU if needed
1514: Otherwise, move back to the CPU */
1515: switch (X->offloadmask) {
1516: case PETSC_OFFLOAD_BOTH:
1517: if (ymask == PETSC_OFFLOAD_CPU) {
1518: VecCUDAResetArray(*Y);
1519: } else if (ymask == PETSC_OFFLOAD_GPU) {
1520: X->offloadmask = PETSC_OFFLOAD_GPU;
1521: }
1522: break;
1523: case PETSC_OFFLOAD_GPU:
1524: if (ymask == PETSC_OFFLOAD_CPU) {
1525: VecCUDAResetArray(*Y);
1526: }
1527: break;
1528: case PETSC_OFFLOAD_CPU:
1529: if (ymask == PETSC_OFFLOAD_GPU) {
1530: VecResetArray(*Y);
1531: }
1532: break;
1533: case PETSC_OFFLOAD_UNALLOCATED:
1534: case PETSC_OFFLOAD_VECKOKKOS:
1535: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"This should not happen");
1536: break;
1537: }
1538: #endif
1539: } else if (iship) {
1540: #if defined(PETSC_HAVE_HIP)
1541: PetscOffloadMask ymask = (*Y)->offloadmask;
1543: /* The offloadmask of X dictates where to move memory
1544: If X GPU data is valid, then move Y data on GPU if needed
1545: Otherwise, move back to the CPU */
1546: switch (X->offloadmask) {
1547: case PETSC_OFFLOAD_BOTH:
1548: if (ymask == PETSC_OFFLOAD_CPU) {
1549: VecHIPResetArray(*Y);
1550: } else if (ymask == PETSC_OFFLOAD_GPU) {
1551: X->offloadmask = PETSC_OFFLOAD_GPU;
1552: }
1553: break;
1554: case PETSC_OFFLOAD_GPU:
1555: if (ymask == PETSC_OFFLOAD_CPU) {
1556: VecHIPResetArray(*Y);
1557: }
1558: break;
1559: case PETSC_OFFLOAD_CPU:
1560: if (ymask == PETSC_OFFLOAD_GPU) {
1561: VecResetArray(*Y);
1562: }
1563: break;
1564: case PETSC_OFFLOAD_UNALLOCATED:
1565: case PETSC_OFFLOAD_VECKOKKOS:
1566: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"This should not happen");
1567: break;
1568: }
1569: #endif
1570: } else {
1571: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1572: VecResetArray(*Y);
1573: }
1574: PetscObjectStateIncrease((PetscObject)X);
1575: }
1576: }
1577: VecDestroy(Y);
1578: }
1579: return(0);
1580: }
1582: /*@
1583: VecGetLocalVectorRead - Maps the local portion of a vector into a
1584: vector. You must call VecRestoreLocalVectorRead() when the local
1585: vector is no longer needed.
1587: Not collective.
1589: Input parameter:
1590: . v - The vector for which the local vector is desired.
1592: Output parameter:
1593: . w - Upon exit this contains the local vector.
1595: Level: beginner
1597: Notes:
1598: This function is similar to VecGetArrayRead() which maps the local
1599: portion into a raw pointer. VecGetLocalVectorRead() is usually
1600: almost as efficient as VecGetArrayRead() but in certain circumstances
1601: VecGetLocalVectorRead() can be much more efficient than
1602: VecGetArrayRead(). This is because the construction of a contiguous
1603: array representing the vector data required by VecGetArrayRead() can
1604: be an expensive operation for certain vector types. For example, for
1605: GPU vectors VecGetArrayRead() requires that the data between device
1606: and host is synchronized.
1608: Unlike VecGetLocalVector(), this routine is not collective and
1609: preserves cached information.
1611: .seealso: VecRestoreLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1612: @*/
1613: PetscErrorCode VecGetLocalVectorRead(Vec v,Vec w)
1614: {
1616: PetscScalar *a;
1621: VecCheckSameLocalSize(v,1,w,2);
1622: if (v->ops->getlocalvectorread) {
1623: (*v->ops->getlocalvectorread)(v,w);
1624: } else {
1625: VecGetArrayRead(v,(const PetscScalar**)&a);
1626: VecPlaceArray(w,a);
1627: }
1628: PetscObjectStateIncrease((PetscObject)w);
1629: VecLockReadPush(v);
1630: VecLockReadPush(w);
1631: return(0);
1632: }
1634: /*@
1635: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1636: previously mapped into a vector using VecGetLocalVectorRead().
1638: Not collective.
1640: Input parameter:
1641: + v - The local portion of this vector was previously mapped into w using VecGetLocalVectorRead().
1642: - w - The vector into which the local portion of v was mapped.
1644: Level: beginner
1646: .seealso: VecGetLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1647: @*/
1648: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec w)
1649: {
1651: PetscScalar *a;
1656: if (v->ops->restorelocalvectorread) {
1657: (*v->ops->restorelocalvectorread)(v,w);
1658: } else {
1659: VecGetArrayRead(w,(const PetscScalar**)&a);
1660: VecRestoreArrayRead(v,(const PetscScalar**)&a);
1661: VecResetArray(w);
1662: }
1663: VecLockReadPop(v);
1664: VecLockReadPop(w);
1665: PetscObjectStateIncrease((PetscObject)w);
1666: return(0);
1667: }
1669: /*@
1670: VecGetLocalVector - Maps the local portion of a vector into a
1671: vector.
1673: Collective on v, not collective on w.
1675: Input parameter:
1676: . v - The vector for which the local vector is desired.
1678: Output parameter:
1679: . w - Upon exit this contains the local vector.
1681: Level: beginner
1683: Notes:
1684: This function is similar to VecGetArray() which maps the local
1685: portion into a raw pointer. VecGetLocalVector() is usually about as
1686: efficient as VecGetArray() but in certain circumstances
1687: VecGetLocalVector() can be much more efficient than VecGetArray().
1688: This is because the construction of a contiguous array representing
1689: the vector data required by VecGetArray() can be an expensive
1690: operation for certain vector types. For example, for GPU vectors
1691: VecGetArray() requires that the data between device and host is
1692: synchronized.
1694: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1695: @*/
1696: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1697: {
1699: PetscScalar *a;
1704: VecCheckSameLocalSize(v,1,w,2);
1705: if (v->ops->getlocalvector) {
1706: (*v->ops->getlocalvector)(v,w);
1707: } else {
1708: VecGetArray(v,&a);
1709: VecPlaceArray(w,a);
1710: }
1711: PetscObjectStateIncrease((PetscObject)w);
1712: return(0);
1713: }
1715: /*@
1716: VecRestoreLocalVector - Unmaps the local portion of a vector
1717: previously mapped into a vector using VecGetLocalVector().
1719: Logically collective.
1721: Input parameter:
1722: + v - The local portion of this vector was previously mapped into w using VecGetLocalVector().
1723: - w - The vector into which the local portion of v was mapped.
1725: Level: beginner
1727: .seealso: VecGetLocalVector(), VecGetLocalVectorRead(), VecRestoreLocalVectorRead(), LocalVectorRead(), VecGetArrayRead(), VecGetArray()
1728: @*/
1729: PetscErrorCode VecRestoreLocalVector(Vec v,Vec w)
1730: {
1732: PetscScalar *a;
1737: if (v->ops->restorelocalvector) {
1738: (*v->ops->restorelocalvector)(v,w);
1739: } else {
1740: VecGetArray(w,&a);
1741: VecRestoreArray(v,&a);
1742: VecResetArray(w);
1743: }
1744: PetscObjectStateIncrease((PetscObject)w);
1745: PetscObjectStateIncrease((PetscObject)v);
1746: return(0);
1747: }
1749: /*@C
1750: VecGetArray - Returns a pointer to a contiguous array that contains this
1751: processor's portion of the vector data. For the standard PETSc
1752: vectors, VecGetArray() returns a pointer to the local data array and
1753: does not use any copies. If the underlying vector data is not stored
1754: in a contiguous array this routine will copy the data to a contiguous
1755: array and return a pointer to that. You MUST call VecRestoreArray()
1756: when you no longer need access to the array.
1758: Logically Collective on Vec
1760: Input Parameter:
1761: . x - the vector
1763: Output Parameter:
1764: . a - location to put pointer to the array
1766: Fortran Note:
1767: This routine is used differently from Fortran 77
1768: $ Vec x
1769: $ PetscScalar x_array(1)
1770: $ PetscOffset i_x
1771: $ PetscErrorCode ierr
1772: $ call VecGetArray(x,x_array,i_x,ierr)
1773: $
1774: $ Access first local entry in vector with
1775: $ value = x_array(i_x + 1)
1776: $
1777: $ ...... other code
1778: $ call VecRestoreArray(x,x_array,i_x,ierr)
1779: For Fortran 90 see VecGetArrayF90()
1781: See the Fortran chapter of the users manual and
1782: petsc/src/snes/tutorials/ex5f.F for details.
1784: Level: beginner
1786: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1787: VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1788: @*/
1789: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1790: {
1795: VecSetErrorIfLocked(x,1);
1796: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
1797: (*x->ops->getarray)(x,a);
1798: } else if (x->petscnative) { /* VECSTANDARD */
1799: *a = *((PetscScalar**)x->data);
1800: } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array for vector type \"%s\"",((PetscObject)x)->type_name);
1801: return(0);
1802: }
1804: /*@C
1805: VecRestoreArray - Restores a vector after VecGetArray() has been called.
1807: Logically Collective on Vec
1809: Input Parameters:
1810: + x - the vector
1811: - a - location of pointer to array obtained from VecGetArray()
1813: Level: beginner
1815: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1816: VecGetArrayPair(), VecRestoreArrayPair()
1817: @*/
1818: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1819: {
1824: if (x->ops->restorearray) { /* VECNEST, VECCUDA etc */
1825: (*x->ops->restorearray)(x,a);
1826: } else if (x->petscnative) { /* VECSTANDARD */
1827: /* nothing */
1828: } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot restore array for vector type \"%s\"",((PetscObject)x)->type_name);
1829: if (a) *a = NULL;
1830: PetscObjectStateIncrease((PetscObject)x);
1831: return(0);
1832: }
1833: /*@C
1834: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
1836: Not Collective
1838: Input Parameter:
1839: . x - the vector
1841: Output Parameter:
1842: . a - the array
1844: Level: beginner
1846: Notes:
1847: The array must be returned using a matching call to VecRestoreArrayRead().
1849: Unlike VecGetArray(), this routine is not collective and preserves cached information like vector norms.
1851: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
1852: implementations may require a copy, but must such implementations should cache the contiguous representation so that
1853: only one copy is performed when this routine is called multiple times in sequence.
1855: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1856: @*/
1857: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1858: {
1863: if (x->ops->getarray) { /* VECNEST, VECCUDA, VECKOKKOS etc */
1864: (*x->ops->getarray)(x,(PetscScalar**)a);
1865: } else if (x->petscnative) { /* VECSTANDARD */
1866: *a = *((PetscScalar**)x->data);
1867: } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array read for vector type \"%s\"",((PetscObject)x)->type_name);
1868: return(0);
1869: }
1871: /*@C
1872: VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()
1874: Not Collective
1876: Input Parameters:
1877: + vec - the vector
1878: - array - the array
1880: Level: beginner
1882: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1883: @*/
1884: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1885: {
1890: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
1891: /* nothing */
1892: } else if (x->ops->restorearrayread) { /* VECNEST */
1893: (*x->ops->restorearrayread)(x,a);
1894: } else { /* No one? */
1895: (*x->ops->restorearray)(x,(PetscScalar**)a);
1896: }
1897: if (a) *a = NULL;
1898: return(0);
1899: }
1901: /*@C
1902: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contains this
1903: processor's portion of the vector data. The values in this array are NOT valid, the routine calling this
1904: routine is responsible for putting values into the array; any values it does not set will be invalid
1906: Logically Collective on Vec
1908: Input Parameter:
1909: . x - the vector
1911: Output Parameter:
1912: . a - location to put pointer to the array
1914: Level: intermediate
1916: This is for vectors associate with GPUs, the vector is not copied up before giving access. If you need correct
1917: values in the array use VecGetArray()
1919: Concepts: vector^accessing local values
1921: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1922: VecGetArrayPair(), VecRestoreArrayPair(), VecGetArray(), VecRestoreArrayWrite()
1923: @*/
1924: PetscErrorCode VecGetArrayWrite(Vec x,PetscScalar **a)
1925: {
1930: VecSetErrorIfLocked(x,1);
1931: if (x->ops->getarraywrite) {
1932: (*x->ops->getarraywrite)(x,a);
1933: } else {
1934: VecGetArray(x,a);
1935: }
1936: return(0);
1937: }
1939: /*@C
1940: VecRestoreArrayWrite - Restores a vector after VecGetArrayWrite() has been called.
1942: Logically Collective on Vec
1944: Input Parameters:
1945: + x - the vector
1946: - a - location of pointer to array obtained from VecGetArray()
1948: Level: beginner
1950: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1951: VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite()
1952: @*/
1953: PetscErrorCode VecRestoreArrayWrite(Vec x,PetscScalar **a)
1954: {
1959: if (x->ops->restorearraywrite) {
1960: (*x->ops->restorearraywrite)(x,a);
1961: } else if (x->ops->restorearray) {
1962: (*x->ops->restorearray)(x,a);
1963: }
1964: if (a) *a = NULL;
1965: PetscObjectStateIncrease((PetscObject)x);
1966: return(0);
1967: }
1969: /*@C
1970: VecGetArrays - Returns a pointer to the arrays in a set of vectors
1971: that were created by a call to VecDuplicateVecs(). You MUST call
1972: VecRestoreArrays() when you no longer need access to the array.
1974: Logically Collective on Vec
1976: Input Parameters:
1977: + x - the vectors
1978: - n - the number of vectors
1980: Output Parameter:
1981: . a - location to put pointer to the array
1983: Fortran Note:
1984: This routine is not supported in Fortran.
1986: Level: intermediate
1988: .seealso: VecGetArray(), VecRestoreArrays()
1989: @*/
1990: PetscErrorCode VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1991: {
1993: PetscInt i;
1994: PetscScalar **q;
2000: if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
2001: PetscMalloc1(n,&q);
2002: for (i=0; i<n; ++i) {
2003: VecGetArray(x[i],&q[i]);
2004: }
2005: *a = q;
2006: return(0);
2007: }
2009: /*@C
2010: VecRestoreArrays - Restores a group of vectors after VecGetArrays()
2011: has been called.
2013: Logically Collective on Vec
2015: Input Parameters:
2016: + x - the vector
2017: . n - the number of vectors
2018: - a - location of pointer to arrays obtained from VecGetArrays()
2020: Notes:
2021: For regular PETSc vectors this routine does not involve any copies. For
2022: any special vectors that do not store local vector data in a contiguous
2023: array, this routine will copy the data back into the underlying
2024: vector data structure from the arrays obtained with VecGetArrays().
2026: Fortran Note:
2027: This routine is not supported in Fortran.
2029: Level: intermediate
2031: .seealso: VecGetArrays(), VecRestoreArray()
2032: @*/
2033: PetscErrorCode VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
2034: {
2036: PetscInt i;
2037: PetscScalar **q = *a;
2044: for (i=0; i<n; ++i) {
2045: VecRestoreArray(x[i],&q[i]);
2046: }
2047: PetscFree(q);
2048: return(0);
2049: }
2051: /*@C
2052: VecGetArrayAndMemType - Like VecGetArray(), but if this is a GPU vector and it is currently offloaded to GPU,
2053: the returned pointer will be a GPU pointer to the GPU memory that contains this processor's portion of the
2054: vector data. Otherwise, it functions as VecGetArray().
2056: Logically Collective on Vec
2058: Input Parameter:
2059: . x - the vector
2061: Output Parameters:
2062: + a - location to put pointer to the array
2063: - mtype - memory type of the array
2065: Level: beginner
2067: .seealso: VecRestoreArrayAndMemType(), VecRestoreArrayAndMemType(), VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(),
2068: VecPlaceArray(), VecGetArray2d(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
2069: @*/
2070: PetscErrorCode VecGetArrayAndMemType(Vec x,PetscScalar **a,PetscMemType *mtype)
2071: {
2077: VecSetErrorIfLocked(x,1);
2078: if (x->ops->getarrayandmemtype) { /* VECCUDA, VECKOKKOS etc */
2079: (*x->ops->getarrayandmemtype)(x,a,mtype);
2080: } else { /* VECSTANDARD, VECNEST, VECVIENNACL */
2081: VecGetArray(x,a);
2082: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2083: }
2084: return(0);
2085: }
2087: /*@C
2088: VecRestoreArrayAndMemType - Restores a vector after VecGetArrayAndMemType() has been called.
2090: Logically Collective on Vec
2092: Input Parameters:
2093: + x - the vector
2094: - a - location of pointer to array obtained from VecGetArrayAndMemType()
2096: Level: beginner
2098: .seealso: VecGetArrayAndMemType(), VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(),
2099: VecPlaceArray(), VecRestoreArray2d(), VecGetArrayPair(), VecRestoreArrayPair()
2100: @*/
2101: PetscErrorCode VecRestoreArrayAndMemType(Vec x,PetscScalar **a)
2102: {
2108: if (x->ops->restorearrayandmemtype) { /* VECCUDA, VECKOKKOS etc */
2109: (*x->ops->restorearrayandmemtype)(x,a);
2110: } else if (x->ops->restorearray) { /* VECNEST, VECVIENNACL */
2111: (*x->ops->restorearray)(x,a);
2112: } /* VECSTANDARD does nothing */
2113: if (a) *a = NULL;
2114: PetscObjectStateIncrease((PetscObject)x);
2115: return(0);
2116: }
2118: /*@C
2119: VecGetArrayReadAndMemType - Like VecGetArrayRead(), but if this is a CUDA vector and it is currently offloaded to GPU,
2120: the returned pointer will be a GPU pointer to the GPU memory that contains this processor's portion of the
2121: vector data. Otherwise, it functions as VecGetArrayRead().
2123: Not Collective
2125: Input Parameter:
2126: . x - the vector
2128: Output Parameters:
2129: + a - the array
2130: - mtype - memory type of the array
2132: Level: beginner
2134: Notes:
2135: The array must be returned using a matching call to VecRestoreArrayReadAndMemType().
2138: .seealso: VecRestoreArrayReadAndMemType(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayAndMemType()
2139: @*/
2140: PetscErrorCode VecGetArrayReadAndMemType(Vec x,const PetscScalar **a,PetscMemType *mtype)
2141: {
2147: #if defined(PETSC_USE_DEBUG)
2148: if (x->ops->getarrayreadandmemtype) SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Not expected vector type \"%s\" has ops->getarrayreadandmemtype",((PetscObject)x)->type_name);
2149: #endif
2151: if (x->ops->getarrayandmemtype) { /* VECCUDA, VECKOKKOS etc, though they are also petscnative */
2152: (*x->ops->getarrayandmemtype)(x,(PetscScalar**)a,mtype);
2153: } else if (x->ops->getarray) { /* VECNEST, VECVIENNACL */
2154: (*x->ops->getarray)(x,(PetscScalar**)a);
2155: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2156: } else if (x->petscnative) { /* VECSTANDARD */
2157: *a = *((PetscScalar**)x->data);
2158: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2159: } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array read in place for vector type \"%s\"",((PetscObject)x)->type_name);
2160: return(0);
2161: }
2163: /*@C
2164: VecRestoreArrayReadAndMemType - Restore array obtained with VecGetArrayReadAndMemType()
2166: Not Collective
2168: Input Parameters:
2169: + vec - the vector
2170: - array - the array
2172: Level: beginner
2174: .seealso: VecGetArrayReadAndMemType(), VecRestoreArrayAndMemType(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
2175: @*/
2176: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x,const PetscScalar **a)
2177: {
2183: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS, VECVIENNACL etc */
2184: /* nothing */
2185: } else if (x->ops->restorearrayread) { /* VECNEST */
2186: (*x->ops->restorearrayread)(x,a);
2187: } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot restore array read in place for vector type \"%s\"",((PetscObject)x)->type_name);
2188: if (a) *a = NULL;
2189: return(0);
2190: }
2192: /*@
2193: VecPlaceArray - Allows one to replace the array in a vector with an
2194: array provided by the user. This is useful to avoid copying an array
2195: into a vector.
2197: Not Collective
2199: Input Parameters:
2200: + vec - the vector
2201: - array - the array
2203: Notes:
2204: You can return to the original array with a call to VecResetArray()
2206: Level: developer
2208: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecResetArray()
2210: @*/
2211: PetscErrorCode VecPlaceArray(Vec vec,const PetscScalar array[])
2212: {
2219: if (vec->ops->placearray) {
2220: (*vec->ops->placearray)(vec,array);
2221: } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
2222: PetscObjectStateIncrease((PetscObject)vec);
2223: return(0);
2224: }
2226: /*@C
2227: VecReplaceArray - Allows one to replace the array in a vector with an
2228: array provided by the user. This is useful to avoid copying an array
2229: into a vector.
2231: Not Collective
2233: Input Parameters:
2234: + vec - the vector
2235: - array - the array
2237: Notes:
2238: This permanently replaces the array and frees the memory associated
2239: with the old array.
2241: The memory passed in MUST be obtained with PetscMalloc() and CANNOT be
2242: freed by the user. It will be freed when the vector is destroyed.
2244: Not supported from Fortran
2246: Level: developer
2248: .seealso: VecGetArray(), VecRestoreArray(), VecPlaceArray(), VecResetArray()
2250: @*/
2251: PetscErrorCode VecReplaceArray(Vec vec,const PetscScalar array[])
2252: {
2258: if (vec->ops->replacearray) {
2259: (*vec->ops->replacearray)(vec,array);
2260: } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
2261: PetscObjectStateIncrease((PetscObject)vec);
2262: return(0);
2263: }
2266: /*@C
2267: VecCUDAGetArray - Provides access to the CUDA buffer inside a vector.
2269: This function has semantics similar to VecGetArray(): the pointer
2270: returned by this function points to a consistent view of the vector
2271: data. This may involve a copy operation of data from the host to the
2272: device if the data on the device is out of date. If the device
2273: memory hasn't been allocated previously it will be allocated as part
2274: of this function call. VecCUDAGetArray() assumes that
2275: the user will modify the vector data. This is similar to
2276: intent(inout) in fortran.
2278: The CUDA device pointer has to be released by calling
2279: VecCUDARestoreArray(). Upon restoring the vector data
2280: the data on the host will be marked as out of date. A subsequent
2281: access of the host data will thus incur a data transfer from the
2282: device to the host.
2285: Input Parameter:
2286: . v - the vector
2288: Output Parameter:
2289: . a - the CUDA device pointer
2291: Fortran note:
2292: This function is not currently available from Fortran.
2294: Level: intermediate
2296: .seealso: VecCUDARestoreArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2297: @*/
2298: PETSC_EXTERN PetscErrorCode VecCUDAGetArray(Vec v, PetscScalar **a)
2299: {
2302: #if defined(PETSC_HAVE_CUDA)
2303: {
2305: VecCUDACopyToGPU(v);
2306: *a = ((Vec_CUDA*)v->spptr)->GPUarray;
2307: }
2308: #endif
2309: return(0);
2310: }
2312: /*@C
2313: VecCUDARestoreArray - Restore a CUDA device pointer previously acquired with VecCUDAGetArray().
2315: This marks the host data as out of date. Subsequent access to the
2316: vector data on the host side with for instance VecGetArray() incurs a
2317: data transfer.
2319: Input Parameter:
2320: + v - the vector
2321: - a - the CUDA device pointer. This pointer is invalid after
2322: VecCUDARestoreArray() returns.
2324: Fortran note:
2325: This function is not currently available from Fortran.
2327: Level: intermediate
2329: .seealso: VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2330: @*/
2331: PETSC_EXTERN PetscErrorCode VecCUDARestoreArray(Vec v, PetscScalar **a)
2332: {
2337: #if defined(PETSC_HAVE_CUDA)
2338: v->offloadmask = PETSC_OFFLOAD_GPU;
2339: #endif
2340: PetscObjectStateIncrease((PetscObject)v);
2341: return(0);
2342: }
2344: /*@C
2345: VecCUDAGetArrayRead - Provides read access to the CUDA buffer inside a vector.
2347: This function is analogous to VecGetArrayRead(): The pointer
2348: returned by this function points to a consistent view of the vector
2349: data. This may involve a copy operation of data from the host to the
2350: device if the data on the device is out of date. If the device
2351: memory hasn't been allocated previously it will be allocated as part
2352: of this function call. VecCUDAGetArrayRead() assumes that the
2353: user will not modify the vector data. This is analgogous to
2354: intent(in) in Fortran.
2356: The CUDA device pointer has to be released by calling
2357: VecCUDARestoreArrayRead(). If the data on the host side was
2358: previously up to date it will remain so, i.e. data on both the device
2359: and the host is up to date. Accessing data on the host side does not
2360: incur a device to host data transfer.
2362: Input Parameter:
2363: . v - the vector
2365: Output Parameter:
2366: . a - the CUDA pointer.
2368: Fortran note:
2369: This function is not currently available from Fortran.
2371: Level: intermediate
2373: .seealso: VecCUDARestoreArrayRead(), VecCUDAGetArray(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2374: @*/
2375: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayRead(Vec v,const PetscScalar** a)
2376: {
2379: VecCUDAGetArray(v,(PetscScalar**)a);
2380: return(0);
2381: }
2383: /*@C
2384: VecCUDARestoreArrayRead - Restore a CUDA device pointer previously acquired with VecCUDAGetArrayRead().
2386: If the data on the host side was previously up to date it will remain
2387: so, i.e. data on both the device and the host is up to date.
2388: Accessing data on the host side e.g. with VecGetArray() does not
2389: incur a device to host data transfer.
2391: Input Parameter:
2392: + v - the vector
2393: - a - the CUDA device pointer. This pointer is invalid after
2394: VecCUDARestoreArrayRead() returns.
2396: Fortran note:
2397: This function is not currently available from Fortran.
2399: Level: intermediate
2401: .seealso: VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2402: @*/
2403: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayRead(Vec v, const PetscScalar **a)
2404: {
2407: *a = NULL;
2408: return(0);
2409: }
2411: /*@C
2412: VecCUDAGetArrayWrite - Provides write access to the CUDA buffer inside a vector.
2414: The data pointed to by the device pointer is uninitialized. The user
2415: may not read from this data. Furthermore, the entire array needs to
2416: be filled by the user to obtain well-defined behaviour. The device
2417: memory will be allocated by this function if it hasn't been allocated
2418: previously. This is analogous to intent(out) in Fortran.
2420: The device pointer needs to be released with
2421: VecCUDARestoreArrayWrite(). When the pointer is released the
2422: host data of the vector is marked as out of data. Subsequent access
2423: of the host data with e.g. VecGetArray() incurs a device to host data
2424: transfer.
2427: Input Parameter:
2428: . v - the vector
2430: Output Parameter:
2431: . a - the CUDA pointer
2433: Fortran note:
2434: This function is not currently available from Fortran.
2436: Level: advanced
2438: .seealso: VecCUDARestoreArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2439: @*/
2440: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayWrite(Vec v, PetscScalar **a)
2441: {
2444: #if defined(PETSC_HAVE_CUDA)
2445: {
2447: VecCUDAAllocateCheck(v);
2448: *a = ((Vec_CUDA*)v->spptr)->GPUarray;
2449: }
2450: #endif
2451: return(0);
2452: }
2454: /*@C
2455: VecCUDARestoreArrayWrite - Restore a CUDA device pointer previously acquired with VecCUDAGetArrayWrite().
2457: Data on the host will be marked as out of date. Subsequent access of
2458: the data on the host side e.g. with VecGetArray() will incur a device
2459: to host data transfer.
2461: Input Parameter:
2462: + v - the vector
2463: - a - the CUDA device pointer. This pointer is invalid after
2464: VecCUDARestoreArrayWrite() returns.
2466: Fortran note:
2467: This function is not currently available from Fortran.
2469: Level: intermediate
2471: .seealso: VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2472: @*/
2473: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayWrite(Vec v, PetscScalar **a)
2474: {
2479: #if defined(PETSC_HAVE_CUDA)
2480: v->offloadmask = PETSC_OFFLOAD_GPU;
2481: if (a) *a = NULL;
2482: #endif
2483: PetscObjectStateIncrease((PetscObject)v);
2484: return(0);
2485: }
2487: /*@C
2488: VecCUDAPlaceArray - Allows one to replace the GPU array in a vector with a
2489: GPU array provided by the user. This is useful to avoid copying an
2490: array into a vector.
2492: Not Collective
2494: Input Parameters:
2495: + vec - the vector
2496: - array - the GPU array
2498: Notes:
2499: You can return to the original GPU array with a call to VecCUDAResetArray()
2500: It is not possible to use VecCUDAPlaceArray() and VecPlaceArray() at the
2501: same time on the same vector.
2503: Level: developer
2505: .seealso: VecPlaceArray(), VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecResetArray(), VecCUDAResetArray(), VecCUDAReplaceArray()
2507: @*/
2508: PetscErrorCode VecCUDAPlaceArray(Vec vin,const PetscScalar a[])
2509: {
2514: #if defined(PETSC_HAVE_CUDA)
2515: VecCUDACopyToGPU(vin);
2516: 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()");
2517: ((Vec_Seq*)vin->data)->unplacedarray = (PetscScalar *) ((Vec_CUDA*)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2518: ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2519: vin->offloadmask = PETSC_OFFLOAD_GPU;
2520: #endif
2521: PetscObjectStateIncrease((PetscObject)vin);
2522: return(0);
2523: }
2525: /*@C
2526: VecCUDAReplaceArray - Allows one to replace the GPU array in a vector
2527: with a GPU array provided by the user. This is useful to avoid copying
2528: a GPU array into a vector.
2530: Not Collective
2532: Input Parameters:
2533: + vec - the vector
2534: - array - the GPU array
2536: Notes:
2537: This permanently replaces the GPU array and frees the memory associated
2538: with the old GPU array.
2540: The memory passed in CANNOT be freed by the user. It will be freed
2541: when the vector is destroyed.
2543: Not supported from Fortran
2545: Level: developer
2547: .seealso: VecGetArray(), VecRestoreArray(), VecPlaceArray(), VecResetArray(), VecCUDAResetArray(), VecCUDAPlaceArray(), VecReplaceArray()
2549: @*/
2550: PetscErrorCode VecCUDAReplaceArray(Vec vin,const PetscScalar a[])
2551: {
2552: #if defined(PETSC_HAVE_CUDA)
2553: cudaError_t err;
2554: #endif
2559: #if defined(PETSC_HAVE_CUDA)
2560: if (((Vec_CUDA*)vin->spptr)->GPUarray_allocated) {
2561: err = cudaFree(((Vec_CUDA*)vin->spptr)->GPUarray_allocated);CHKERRCUDA(err);
2562: }
2563: ((Vec_CUDA*)vin->spptr)->GPUarray_allocated = ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2564: vin->offloadmask = PETSC_OFFLOAD_GPU;
2565: #endif
2566: PetscObjectStateIncrease((PetscObject)vin);
2567: return(0);
2568: }
2570: /*@C
2571: VecCUDAResetArray - Resets a vector to use its default memory. Call this
2572: after the use of VecCUDAPlaceArray().
2574: Not Collective
2576: Input Parameters:
2577: . vec - the vector
2579: Level: developer
2581: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecPlaceArray(), VecResetArray(), VecCUDAPlaceArray(), VecCUDAReplaceArray()
2583: @*/
2584: PetscErrorCode VecCUDAResetArray(Vec vin)
2585: {
2590: #if defined(PETSC_HAVE_CUDA)
2591: VecCUDACopyToGPU(vin);
2592: ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar *) ((Vec_Seq*)vin->data)->unplacedarray;
2593: ((Vec_Seq*)vin->data)->unplacedarray = 0;
2594: vin->offloadmask = PETSC_OFFLOAD_GPU;
2595: #endif
2596: PetscObjectStateIncrease((PetscObject)vin);
2597: return(0);
2598: }
2600: /*@C
2601: VecHIPGetArray - Provides access to the HIP buffer inside a vector.
2603: This function has semantics similar to VecGetArray(): the pointer
2604: returned by this function points to a consistent view of the vector
2605: data. This may involve a copy operation of data from the host to the
2606: device if the data on the device is out of date. If the device
2607: memory hasn't been allocated previously it will be allocated as part
2608: of this function call. VecHIPGetArray() assumes that
2609: the user will modify the vector data. This is similar to
2610: intent(inout) in fortran.
2612: The HIP device pointer has to be released by calling
2613: VecHIPRestoreArray(). Upon restoring the vector data
2614: the data on the host will be marked as out of date. A subsequent
2615: access of the host data will thus incur a data transfer from the
2616: device to the host.
2619: Input Parameter:
2620: . v - the vector
2622: Output Parameter:
2623: . a - the HIP device pointer
2625: Fortran note:
2626: This function is not currently available from Fortran.
2628: Level: intermediate
2630: .seealso: VecHIPRestoreArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2631: @*/
2632: PETSC_EXTERN PetscErrorCode VecHIPGetArray(Vec v, PetscScalar **a)
2633: {
2634: #if defined(PETSC_HAVE_HIP)
2636: #endif
2640: #if defined(PETSC_HAVE_HIP)
2641: *a = 0;
2642: VecHIPCopyToGPU(v);
2643: *a = ((Vec_HIP*)v->spptr)->GPUarray;
2644: #endif
2645: return(0);
2646: }
2648: /*@C
2649: VecHIPRestoreArray - Restore a HIP device pointer previously acquired with VecHIPGetArray().
2651: This marks the host data as out of date. Subsequent access to the
2652: vector data on the host side with for instance VecGetArray() incurs a
2653: data transfer.
2655: Input Parameter:
2656: + v - the vector
2657: - a - the HIP device pointer. This pointer is invalid after
2658: VecHIPRestoreArray() returns.
2660: Fortran note:
2661: This function is not currently available from Fortran.
2663: Level: intermediate
2665: .seealso: VecHIPGetArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2666: @*/
2667: PETSC_EXTERN PetscErrorCode VecHIPRestoreArray(Vec v, PetscScalar **a)
2668: {
2673: #if defined(PETSC_HAVE_HIP)
2674: v->offloadmask = PETSC_OFFLOAD_GPU;
2675: #endif
2677: PetscObjectStateIncrease((PetscObject)v);
2678: return(0);
2679: }
2681: /*@C
2682: VecHIPGetArrayRead - Provides read access to the HIP buffer inside a vector.
2684: This function is analogous to VecGetArrayRead(): The pointer
2685: returned by this function points to a consistent view of the vector
2686: data. This may involve a copy operation of data from the host to the
2687: device if the data on the device is out of date. If the device
2688: memory hasn't been allocated previously it will be allocated as part
2689: of this function call. VecHIPGetArrayRead() assumes that the
2690: user will not modify the vector data. This is analgogous to
2691: intent(in) in Fortran.
2693: The HIP device pointer has to be released by calling
2694: VecHIPRestoreArrayRead(). If the data on the host side was
2695: previously up to date it will remain so, i.e. data on both the device
2696: and the host is up to date. Accessing data on the host side does not
2697: incur a device to host data transfer.
2699: Input Parameter:
2700: . v - the vector
2702: Output Parameter:
2703: . a - the HIP pointer.
2705: Fortran note:
2706: This function is not currently available from Fortran.
2708: Level: intermediate
2710: .seealso: VecHIPRestoreArrayRead(), VecHIPGetArray(), VecHIPGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2711: @*/
2712: PETSC_EXTERN PetscErrorCode VecHIPGetArrayRead(Vec v, const PetscScalar **a)
2713: {
2714: #if defined(PETSC_HAVE_HIP)
2716: #endif
2720: #if defined(PETSC_HAVE_HIP)
2721: *a = 0;
2722: VecHIPCopyToGPU(v);
2723: *a = ((Vec_HIP*)v->spptr)->GPUarray;
2724: #endif
2725: return(0);
2726: }
2728: /*@C
2729: VecHIPRestoreArrayRead - Restore a HIP device pointer previously acquired with VecHIPGetArrayRead().
2731: If the data on the host side was previously up to date it will remain
2732: so, i.e. data on both the device and the host is up to date.
2733: Accessing data on the host side e.g. with VecGetArray() does not
2734: incur a device to host data transfer.
2736: Input Parameter:
2737: + v - the vector
2738: - a - the HIP device pointer. This pointer is invalid after
2739: VecHIPRestoreArrayRead() returns.
2741: Fortran note:
2742: This function is not currently available from Fortran.
2744: Level: intermediate
2746: .seealso: VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecHIPGetArray(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2747: @*/
2748: PETSC_EXTERN PetscErrorCode VecHIPRestoreArrayRead(Vec v, const PetscScalar **a)
2749: {
2752: *a = NULL;
2753: return(0);
2754: }
2756: /*@C
2757: VecHIPGetArrayWrite - Provides write access to the HIP buffer inside a vector.
2759: The data pointed to by the device pointer is uninitialized. The user
2760: may not read from this data. Furthermore, the entire array needs to
2761: be filled by the user to obtain well-defined behaviour. The device
2762: memory will be allocated by this function if it hasn't been allocated
2763: previously. This is analogous to intent(out) in Fortran.
2765: The device pointer needs to be released with
2766: VecHIPRestoreArrayWrite(). When the pointer is released the
2767: host data of the vector is marked as out of data. Subsequent access
2768: of the host data with e.g. VecGetArray() incurs a device to host data
2769: transfer.
2772: Input Parameter:
2773: . v - the vector
2775: Output Parameter:
2776: . a - the HIP pointer
2778: Fortran note:
2779: This function is not currently available from Fortran.
2781: Level: advanced
2783: .seealso: VecHIPRestoreArrayWrite(), VecHIPGetArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2784: @*/
2785: PETSC_EXTERN PetscErrorCode VecHIPGetArrayWrite(Vec v, PetscScalar **a)
2786: {
2787: #if defined(PETSC_HAVE_HIP)
2789: #endif
2793: #if defined(PETSC_HAVE_HIP)
2794: *a = 0;
2795: VecHIPAllocateCheck(v);
2796: *a = ((Vec_HIP*)v->spptr)->GPUarray;
2797: #endif
2798: return(0);
2799: }
2801: /*@C
2802: VecHIPRestoreArrayWrite - Restore a HIP device pointer previously acquired with VecHIPGetArrayWrite().
2804: Data on the host will be marked as out of date. Subsequent access of
2805: the data on the host side e.g. with VecGetArray() will incur a device
2806: to host data transfer.
2808: Input Parameter:
2809: + v - the vector
2810: - a - the HIP device pointer. This pointer is invalid after
2811: VecHIPRestoreArrayWrite() returns.
2813: Fortran note:
2814: This function is not currently available from Fortran.
2816: Level: intermediate
2818: .seealso: VecHIPGetArrayWrite(), VecHIPGetArray(), VecHIPGetArrayRead(), VecHIPGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2819: @*/
2820: PETSC_EXTERN PetscErrorCode VecHIPRestoreArrayWrite(Vec v, PetscScalar **a)
2821: {
2826: #if defined(PETSC_HAVE_HIP)
2827: v->offloadmask = PETSC_OFFLOAD_GPU;
2828: #endif
2830: PetscObjectStateIncrease((PetscObject)v);
2831: return(0);
2832: }
2834: /*@C
2835: VecHIPPlaceArray - Allows one to replace the GPU array in a vector with a
2836: GPU array provided by the user. This is useful to avoid copying an
2837: array into a vector.
2839: Not Collective
2841: Input Parameters:
2842: + vec - the vector
2843: - array - the GPU array
2845: Notes:
2846: You can return to the original GPU array with a call to VecHIPResetArray()
2847: It is not possible to use VecHIPPlaceArray() and VecPlaceArray() at the
2848: same time on the same vector.
2850: Level: developer
2852: .seealso: VecPlaceArray(), VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecResetArray(), VecHIPResetArray(), VecHIPReplaceArray()
2854: @*/
2855: PetscErrorCode VecHIPPlaceArray(Vec vin,const PetscScalar a[])
2856: {
2861: #if defined(PETSC_HAVE_HIP)
2862: VecHIPCopyToGPU(vin);
2863: if (((Vec_Seq*)vin->data)->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"VecHIPPlaceArray()/VecPlaceArray() was already called on this vector, without a call to VecHIPResetArray()/VecResetArray()");
2864: ((Vec_Seq*)vin->data)->unplacedarray = (PetscScalar *) ((Vec_HIP*)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2865: ((Vec_HIP*)vin->spptr)->GPUarray = (PetscScalar*)a;
2866: vin->offloadmask = PETSC_OFFLOAD_GPU;
2867: #endif
2868: PetscObjectStateIncrease((PetscObject)vin);
2869: return(0);
2870: }
2872: /*@C
2873: VecHIPReplaceArray - Allows one to replace the GPU array in a vector
2874: with a GPU array provided by the user. This is useful to avoid copying
2875: a GPU array into a vector.
2877: Not Collective
2879: Input Parameters:
2880: + vec - the vector
2881: - array - the GPU array
2883: Notes:
2884: This permanently replaces the GPU array and frees the memory associated
2885: with the old GPU array.
2887: The memory passed in CANNOT be freed by the user. It will be freed
2888: when the vector is destroyed.
2890: Not supported from Fortran
2892: Level: developer
2894: .seealso: VecGetArray(), VecRestoreArray(), VecPlaceArray(), VecResetArray(), VecHIPResetArray(), VecHIPPlaceArray(), VecReplaceArray()
2896: @*/
2897: PetscErrorCode VecHIPReplaceArray(Vec vin,const PetscScalar a[])
2898: {
2899: #if defined(PETSC_HAVE_HIP)
2900: hipError_t err;
2901: #endif
2906: #if defined(PETSC_HAVE_HIP)
2907: err = hipFree(((Vec_HIP*)vin->spptr)->GPUarray);CHKERRHIP(err);
2908: ((Vec_HIP*)vin->spptr)->GPUarray = (PetscScalar*)a;
2909: vin->offloadmask = PETSC_OFFLOAD_GPU;
2910: #endif
2911: PetscObjectStateIncrease((PetscObject)vin);
2912: return(0);
2913: }
2915: /*@C
2916: VecHIPResetArray - Resets a vector to use its default memory. Call this
2917: after the use of VecHIPPlaceArray().
2919: Not Collective
2921: Input Parameters:
2922: . vec - the vector
2924: Level: developer
2926: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecPlaceArray(), VecResetArray(), VecHIPPlaceArray(), VecHIPReplaceArray()
2928: @*/
2929: PetscErrorCode VecHIPResetArray(Vec vin)
2930: {
2935: #if defined(PETSC_HAVE_HIP)
2936: VecHIPCopyToGPU(vin);
2937: ((Vec_HIP*)vin->spptr)->GPUarray = (PetscScalar *) ((Vec_Seq*)vin->data)->unplacedarray;
2938: ((Vec_Seq*)vin->data)->unplacedarray = 0;
2939: vin->offloadmask = PETSC_OFFLOAD_GPU;
2940: #endif
2941: PetscObjectStateIncrease((PetscObject)vin);
2942: return(0);
2943: }
2945: /*MC
2946: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
2947: and makes them accessible via a Fortran90 pointer.
2949: Synopsis:
2950: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
2952: Collective on Vec
2954: Input Parameters:
2955: + x - a vector to mimic
2956: - n - the number of vectors to obtain
2958: Output Parameters:
2959: + y - Fortran90 pointer to the array of vectors
2960: - ierr - error code
2962: Example of Usage:
2963: .vb
2964: #include <petsc/finclude/petscvec.h>
2965: use petscvec
2967: Vec x
2968: Vec, pointer :: y(:)
2969: ....
2970: call VecDuplicateVecsF90(x,2,y,ierr)
2971: call VecSet(y(2),alpha,ierr)
2972: call VecSet(y(2),alpha,ierr)
2973: ....
2974: call VecDestroyVecsF90(2,y,ierr)
2975: .ve
2977: Notes:
2978: Not yet supported for all F90 compilers
2980: Use VecDestroyVecsF90() to free the space.
2982: Level: beginner
2984: .seealso: VecDestroyVecsF90(), VecDuplicateVecs()
2986: M*/
2988: /*MC
2989: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2990: VecGetArrayF90().
2992: Synopsis:
2993: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2995: Logically Collective on Vec
2997: Input Parameters:
2998: + x - vector
2999: - xx_v - the Fortran90 pointer to the array
3001: Output Parameter:
3002: . ierr - error code
3004: Example of Usage:
3005: .vb
3006: #include <petsc/finclude/petscvec.h>
3007: use petscvec
3009: PetscScalar, pointer :: xx_v(:)
3010: ....
3011: call VecGetArrayF90(x,xx_v,ierr)
3012: xx_v(3) = a
3013: call VecRestoreArrayF90(x,xx_v,ierr)
3014: .ve
3016: Level: beginner
3018: .seealso: VecGetArrayF90(), VecGetArray(), VecRestoreArray(), UsingFortran, VecRestoreArrayReadF90()
3020: M*/
3022: /*MC
3023: VecDestroyVecsF90 - Frees a block of vectors obtained with VecDuplicateVecsF90().
3025: Synopsis:
3026: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
3028: Collective on Vec
3030: Input Parameters:
3031: + n - the number of vectors previously obtained
3032: - x - pointer to array of vector pointers
3034: Output Parameter:
3035: . ierr - error code
3037: Notes:
3038: Not yet supported for all F90 compilers
3040: Level: beginner
3042: .seealso: VecDestroyVecs(), VecDuplicateVecsF90()
3044: M*/
3046: /*MC
3047: VecGetArrayF90 - Accesses a vector array from Fortran90. For default PETSc
3048: vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
3049: this routine is implementation dependent. You MUST call VecRestoreArrayF90()
3050: when you no longer need access to the array.
3052: Synopsis:
3053: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
3055: Logically Collective on Vec
3057: Input Parameter:
3058: . x - vector
3060: Output Parameters:
3061: + xx_v - the Fortran90 pointer to the array
3062: - ierr - error code
3064: Example of Usage:
3065: .vb
3066: #include <petsc/finclude/petscvec.h>
3067: use petscvec
3069: PetscScalar, pointer :: xx_v(:)
3070: ....
3071: call VecGetArrayF90(x,xx_v,ierr)
3072: xx_v(3) = a
3073: call VecRestoreArrayF90(x,xx_v,ierr)
3074: .ve
3076: If you ONLY intend to read entries from the array and not change any entries you should use VecGetArrayReadF90().
3078: Level: beginner
3080: .seealso: VecRestoreArrayF90(), VecGetArray(), VecRestoreArray(), VecGetArrayReadF90(), UsingFortran
3082: M*/
3084: /*MC
3085: VecGetArrayReadF90 - Accesses a read only array from Fortran90. For default PETSc
3086: vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
3087: this routine is implementation dependent. You MUST call VecRestoreArrayReadF90()
3088: when you no longer need access to the array.
3090: Synopsis:
3091: VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
3093: Logically Collective on Vec
3095: Input Parameter:
3096: . x - vector
3098: Output Parameters:
3099: + xx_v - the Fortran90 pointer to the array
3100: - ierr - error code
3102: Example of Usage:
3103: .vb
3104: #include <petsc/finclude/petscvec.h>
3105: use petscvec
3107: PetscScalar, pointer :: xx_v(:)
3108: ....
3109: call VecGetArrayReadF90(x,xx_v,ierr)
3110: a = xx_v(3)
3111: call VecRestoreArrayReadF90(x,xx_v,ierr)
3112: .ve
3114: If you intend to write entries into the array you must use VecGetArrayF90().
3116: Level: beginner
3118: .seealso: VecRestoreArrayReadF90(), VecGetArray(), VecRestoreArray(), VecGetArrayRead(), VecRestoreArrayRead(), VecGetArrayF90(), UsingFortran
3120: M*/
3122: /*MC
3123: VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
3124: VecGetArrayReadF90().
3126: Synopsis:
3127: VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
3129: Logically Collective on Vec
3131: Input Parameters:
3132: + x - vector
3133: - xx_v - the Fortran90 pointer to the array
3135: Output Parameter:
3136: . ierr - error code
3138: Example of Usage:
3139: .vb
3140: #include <petsc/finclude/petscvec.h>
3141: use petscvec
3143: PetscScalar, pointer :: xx_v(:)
3144: ....
3145: call VecGetArrayReadF90(x,xx_v,ierr)
3146: a = xx_v(3)
3147: call VecRestoreArrayReadF90(x,xx_v,ierr)
3148: .ve
3150: Level: beginner
3152: .seealso: VecGetArrayReadF90(), VecGetArray(), VecRestoreArray(), VecGetArrayRead(), VecRestoreArrayRead(),UsingFortran, VecRestoreArrayF90()
3154: M*/
3156: /*@C
3157: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
3158: processor's portion of the vector data. You MUST call VecRestoreArray2d()
3159: when you no longer need access to the array.
3161: Logically Collective
3163: Input Parameter:
3164: + x - the vector
3165: . m - first dimension of two dimensional array
3166: . n - second dimension of two dimensional array
3167: . mstart - first index you will use in first coordinate direction (often 0)
3168: - nstart - first index in the second coordinate direction (often 0)
3170: Output Parameter:
3171: . a - location to put pointer to the array
3173: Level: developer
3175: Notes:
3176: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
3177: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3178: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3179: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
3181: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3183: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3184: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3185: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3186: @*/
3187: PetscErrorCode VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3188: {
3190: PetscInt i,N;
3191: PetscScalar *aa;
3197: VecGetLocalSize(x,&N);
3198: 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);
3199: VecGetArray(x,&aa);
3201: PetscMalloc1(m,a);
3202: for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
3203: *a -= mstart;
3204: return(0);
3205: }
3207: /*@C
3208: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
3209: processor's portion of the vector data. You MUST call VecRestoreArray2dWrite()
3210: when you no longer need access to the array.
3212: Logically Collective
3214: Input Parameter:
3215: + x - the vector
3216: . m - first dimension of two dimensional array
3217: . n - second dimension of two dimensional array
3218: . mstart - first index you will use in first coordinate direction (often 0)
3219: - nstart - first index in the second coordinate direction (often 0)
3221: Output Parameter:
3222: . a - location to put pointer to the array
3224: Level: developer
3226: Notes:
3227: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
3228: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3229: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3230: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
3232: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3234: Concepts: vector^accessing local values as 2d array
3236: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3237: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3238: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3239: @*/
3240: PetscErrorCode VecGetArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3241: {
3243: PetscInt i,N;
3244: PetscScalar *aa;
3250: VecGetLocalSize(x,&N);
3251: 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);
3252: VecGetArrayWrite(x,&aa);
3254: PetscMalloc1(m,a);
3255: for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
3256: *a -= mstart;
3257: return(0);
3258: }
3260: /*@C
3261: VecRestoreArray2d - Restores a vector after VecGetArray2d() has been called.
3263: Logically Collective
3265: Input Parameters:
3266: + x - the vector
3267: . m - first dimension of two dimensional array
3268: . n - second dimension of the two dimensional array
3269: . mstart - first index you will use in first coordinate direction (often 0)
3270: . nstart - first index in the second coordinate direction (often 0)
3271: - a - location of pointer to array obtained from VecGetArray2d()
3273: Level: developer
3275: Notes:
3276: For regular PETSc vectors this routine does not involve any copies. For
3277: any special vectors that do not store local vector data in a contiguous
3278: array, this routine will copy the data back into the underlying
3279: vector data structure from the array obtained with VecGetArray().
3281: This routine actually zeros out the a pointer.
3283: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3284: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3285: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3286: @*/
3287: PetscErrorCode VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3288: {
3290: void *dummy;
3296: dummy = (void*)(*a + mstart);
3297: PetscFree(dummy);
3298: VecRestoreArray(x,NULL);
3299: return(0);
3300: }
3302: /*@C
3303: VecRestoreArray2dWrite - Restores a vector after VecGetArray2dWrite() has been called.
3305: Logically Collective
3307: Input Parameters:
3308: + x - the vector
3309: . m - first dimension of two dimensional array
3310: . n - second dimension of the two dimensional array
3311: . mstart - first index you will use in first coordinate direction (often 0)
3312: . nstart - first index in the second coordinate direction (often 0)
3313: - a - location of pointer to array obtained from VecGetArray2d()
3315: Level: developer
3317: Notes:
3318: For regular PETSc vectors this routine does not involve any copies. For
3319: any special vectors that do not store local vector data in a contiguous
3320: array, this routine will copy the data back into the underlying
3321: vector data structure from the array obtained with VecGetArray().
3323: This routine actually zeros out the a pointer.
3325: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3326: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3327: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3328: @*/
3329: PetscErrorCode VecRestoreArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3330: {
3332: void *dummy;
3338: dummy = (void*)(*a + mstart);
3339: PetscFree(dummy);
3340: VecRestoreArrayWrite(x,NULL);
3341: return(0);
3342: }
3344: /*@C
3345: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
3346: processor's portion of the vector data. You MUST call VecRestoreArray1d()
3347: when you no longer need access to the array.
3349: Logically Collective
3351: Input Parameter:
3352: + x - the vector
3353: . m - first dimension of two dimensional array
3354: - mstart - first index you will use in first coordinate direction (often 0)
3356: Output Parameter:
3357: . a - location to put pointer to the array
3359: Level: developer
3361: Notes:
3362: For a vector obtained from DMCreateLocalVector() mstart are likely
3363: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3364: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
3366: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3368: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3369: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3370: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3371: @*/
3372: PetscErrorCode VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3373: {
3375: PetscInt N;
3381: VecGetLocalSize(x,&N);
3382: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3383: VecGetArray(x,a);
3384: *a -= mstart;
3385: return(0);
3386: }
3388: /*@C
3389: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
3390: processor's portion of the vector data. You MUST call VecRestoreArray1dWrite()
3391: when you no longer need access to the array.
3393: Logically Collective
3395: Input Parameter:
3396: + x - the vector
3397: . m - first dimension of two dimensional array
3398: - mstart - first index you will use in first coordinate direction (often 0)
3400: Output Parameter:
3401: . a - location to put pointer to the array
3403: Level: developer
3405: Notes:
3406: For a vector obtained from DMCreateLocalVector() mstart are likely
3407: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3408: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
3410: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3412: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3413: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3414: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3415: @*/
3416: PetscErrorCode VecGetArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3417: {
3419: PetscInt N;
3425: VecGetLocalSize(x,&N);
3426: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3427: VecGetArrayWrite(x,a);
3428: *a -= mstart;
3429: return(0);
3430: }
3432: /*@C
3433: VecRestoreArray1d - Restores a vector after VecGetArray1d() has been called.
3435: Logically Collective
3437: Input Parameters:
3438: + x - the vector
3439: . m - first dimension of two dimensional array
3440: . mstart - first index you will use in first coordinate direction (often 0)
3441: - a - location of pointer to array obtained from VecGetArray21()
3443: Level: developer
3445: Notes:
3446: For regular PETSc vectors this routine does not involve any copies. For
3447: any special vectors that do not store local vector data in a contiguous
3448: array, this routine will copy the data back into the underlying
3449: vector data structure from the array obtained with VecGetArray1d().
3451: This routine actually zeros out the a pointer.
3453: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3454: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3455: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3456: @*/
3457: PetscErrorCode VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3458: {
3464: VecRestoreArray(x,NULL);
3465: return(0);
3466: }
3468: /*@C
3469: VecRestoreArray1dWrite - Restores a vector after VecGetArray1dWrite() has been called.
3471: Logically Collective
3473: Input Parameters:
3474: + x - the vector
3475: . m - first dimension of two dimensional array
3476: . mstart - first index you will use in first coordinate direction (often 0)
3477: - a - location of pointer to array obtained from VecGetArray21()
3479: Level: developer
3481: Notes:
3482: For regular PETSc vectors this routine does not involve any copies. For
3483: any special vectors that do not store local vector data in a contiguous
3484: array, this routine will copy the data back into the underlying
3485: vector data structure from the array obtained with VecGetArray1d().
3487: This routine actually zeros out the a pointer.
3489: Concepts: vector^accessing local values as 1d array
3491: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3492: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3493: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3494: @*/
3495: PetscErrorCode VecRestoreArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3496: {
3502: VecRestoreArrayWrite(x,NULL);
3503: return(0);
3504: }
3506: /*@C
3507: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
3508: processor's portion of the vector data. You MUST call VecRestoreArray3d()
3509: when you no longer need access to the array.
3511: Logically Collective
3513: Input Parameter:
3514: + x - the vector
3515: . m - first dimension of three dimensional array
3516: . n - second dimension of three dimensional array
3517: . p - third dimension of three dimensional array
3518: . mstart - first index you will use in first coordinate direction (often 0)
3519: . nstart - first index in the second coordinate direction (often 0)
3520: - pstart - first index in the third coordinate direction (often 0)
3522: Output Parameter:
3523: . a - location to put pointer to the array
3525: Level: developer
3527: Notes:
3528: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3529: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3530: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3531: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
3533: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3535: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3536: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3537: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3538: @*/
3539: PetscErrorCode VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3540: {
3542: PetscInt i,N,j;
3543: PetscScalar *aa,**b;
3549: VecGetLocalSize(x,&N);
3550: 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);
3551: VecGetArray(x,&aa);
3553: PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3554: b = (PetscScalar**)((*a) + m);
3555: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3556: for (i=0; i<m; i++)
3557: for (j=0; j<n; j++)
3558: b[i*n+j] = aa + i*n*p + j*p - pstart;
3560: *a -= mstart;
3561: return(0);
3562: }
3564: /*@C
3565: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
3566: processor's portion of the vector data. You MUST call VecRestoreArray3dWrite()
3567: when you no longer need access to the array.
3569: Logically Collective
3571: Input Parameter:
3572: + x - the vector
3573: . m - first dimension of three dimensional array
3574: . n - second dimension of three dimensional array
3575: . p - third dimension of three dimensional array
3576: . mstart - first index you will use in first coordinate direction (often 0)
3577: . nstart - first index in the second coordinate direction (often 0)
3578: - pstart - first index in the third coordinate direction (often 0)
3580: Output Parameter:
3581: . a - location to put pointer to the array
3583: Level: developer
3585: Notes:
3586: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3587: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3588: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3589: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
3591: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3593: Concepts: vector^accessing local values as 3d array
3595: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3596: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3597: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3598: @*/
3599: PetscErrorCode VecGetArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3600: {
3602: PetscInt i,N,j;
3603: PetscScalar *aa,**b;
3609: VecGetLocalSize(x,&N);
3610: 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);
3611: VecGetArrayWrite(x,&aa);
3613: PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3614: b = (PetscScalar**)((*a) + m);
3615: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3616: for (i=0; i<m; i++)
3617: for (j=0; j<n; j++)
3618: b[i*n+j] = aa + i*n*p + j*p - pstart;
3620: *a -= mstart;
3621: return(0);
3622: }
3624: /*@C
3625: VecRestoreArray3d - Restores a vector after VecGetArray3d() has been called.
3627: Logically Collective
3629: Input Parameters:
3630: + x - the vector
3631: . m - first dimension of three dimensional array
3632: . n - second dimension of the three dimensional array
3633: . p - third dimension of the three dimensional array
3634: . mstart - first index you will use in first coordinate direction (often 0)
3635: . nstart - first index in the second coordinate direction (often 0)
3636: . pstart - first index in the third coordinate direction (often 0)
3637: - a - location of pointer to array obtained from VecGetArray3d()
3639: Level: developer
3641: Notes:
3642: For regular PETSc vectors this routine does not involve any copies. For
3643: any special vectors that do not store local vector data in a contiguous
3644: array, this routine will copy the data back into the underlying
3645: vector data structure from the array obtained with VecGetArray().
3647: This routine actually zeros out the a pointer.
3649: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3650: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3651: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3652: @*/
3653: PetscErrorCode VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3654: {
3656: void *dummy;
3662: dummy = (void*)(*a + mstart);
3663: PetscFree(dummy);
3664: VecRestoreArray(x,NULL);
3665: return(0);
3666: }
3668: /*@C
3669: VecRestoreArray3dWrite - Restores a vector after VecGetArray3dWrite() has been called.
3671: Logically Collective
3673: Input Parameters:
3674: + x - the vector
3675: . m - first dimension of three dimensional array
3676: . n - second dimension of the three dimensional array
3677: . p - third dimension of the three dimensional array
3678: . mstart - first index you will use in first coordinate direction (often 0)
3679: . nstart - first index in the second coordinate direction (often 0)
3680: . pstart - first index in the third coordinate direction (often 0)
3681: - a - location of pointer to array obtained from VecGetArray3d()
3683: Level: developer
3685: Notes:
3686: For regular PETSc vectors this routine does not involve any copies. For
3687: any special vectors that do not store local vector data in a contiguous
3688: array, this routine will copy the data back into the underlying
3689: vector data structure from the array obtained with VecGetArray().
3691: This routine actually zeros out the a pointer.
3693: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3694: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3695: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3696: @*/
3697: PetscErrorCode VecRestoreArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3698: {
3700: void *dummy;
3706: dummy = (void*)(*a + mstart);
3707: PetscFree(dummy);
3708: VecRestoreArrayWrite(x,NULL);
3709: return(0);
3710: }
3712: /*@C
3713: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3714: processor's portion of the vector data. You MUST call VecRestoreArray4d()
3715: when you no longer need access to the array.
3717: Logically Collective
3719: Input Parameter:
3720: + x - the vector
3721: . m - first dimension of four dimensional array
3722: . n - second dimension of four dimensional array
3723: . p - third dimension of four dimensional array
3724: . q - fourth dimension of four dimensional array
3725: . mstart - first index you will use in first coordinate direction (often 0)
3726: . nstart - first index in the second coordinate direction (often 0)
3727: . pstart - first index in the third coordinate direction (often 0)
3728: - qstart - first index in the fourth coordinate direction (often 0)
3730: Output Parameter:
3731: . a - location to put pointer to the array
3733: Level: beginner
3735: Notes:
3736: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3737: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3738: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3739: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
3741: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3743: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3744: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3745: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3746: @*/
3747: PetscErrorCode VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3748: {
3750: PetscInt i,N,j,k;
3751: PetscScalar *aa,***b,**c;
3757: VecGetLocalSize(x,&N);
3758: 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);
3759: VecGetArray(x,&aa);
3761: PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3762: b = (PetscScalar***)((*a) + m);
3763: c = (PetscScalar**)(b + m*n);
3764: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3765: for (i=0; i<m; i++)
3766: for (j=0; j<n; j++)
3767: b[i*n+j] = c + i*n*p + j*p - pstart;
3768: for (i=0; i<m; i++)
3769: for (j=0; j<n; j++)
3770: for (k=0; k<p; k++)
3771: c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3772: *a -= mstart;
3773: return(0);
3774: }
3776: /*@C
3777: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3778: processor's portion of the vector data. You MUST call VecRestoreArray4dWrite()
3779: when you no longer need access to the array.
3781: Logically Collective
3783: Input Parameter:
3784: + x - the vector
3785: . m - first dimension of four dimensional array
3786: . n - second dimension of four dimensional array
3787: . p - third dimension of four dimensional array
3788: . q - fourth dimension of four dimensional array
3789: . mstart - first index you will use in first coordinate direction (often 0)
3790: . nstart - first index in the second coordinate direction (often 0)
3791: . pstart - first index in the third coordinate direction (often 0)
3792: - qstart - first index in the fourth coordinate direction (often 0)
3794: Output Parameter:
3795: . a - location to put pointer to the array
3797: Level: beginner
3799: Notes:
3800: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3801: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3802: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3803: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
3805: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3807: Concepts: vector^accessing local values as 3d array
3809: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3810: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3811: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3812: @*/
3813: PetscErrorCode VecGetArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3814: {
3816: PetscInt i,N,j,k;
3817: PetscScalar *aa,***b,**c;
3823: VecGetLocalSize(x,&N);
3824: 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);
3825: VecGetArrayWrite(x,&aa);
3827: PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3828: b = (PetscScalar***)((*a) + m);
3829: c = (PetscScalar**)(b + m*n);
3830: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3831: for (i=0; i<m; i++)
3832: for (j=0; j<n; j++)
3833: b[i*n+j] = c + i*n*p + j*p - pstart;
3834: for (i=0; i<m; i++)
3835: for (j=0; j<n; j++)
3836: for (k=0; k<p; k++)
3837: c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3838: *a -= mstart;
3839: return(0);
3840: }
3842: /*@C
3843: VecRestoreArray4d - Restores a vector after VecGetArray3d() has been called.
3845: Logically Collective
3847: Input Parameters:
3848: + x - the vector
3849: . m - first dimension of four dimensional array
3850: . n - second dimension of the four dimensional array
3851: . p - third dimension of the four dimensional array
3852: . q - fourth dimension of the four dimensional array
3853: . mstart - first index you will use in first coordinate direction (often 0)
3854: . nstart - first index in the second coordinate direction (often 0)
3855: . pstart - first index in the third coordinate direction (often 0)
3856: . qstart - first index in the fourth coordinate direction (often 0)
3857: - a - location of pointer to array obtained from VecGetArray4d()
3859: Level: beginner
3861: Notes:
3862: For regular PETSc vectors this routine does not involve any copies. For
3863: any special vectors that do not store local vector data in a contiguous
3864: array, this routine will copy the data back into the underlying
3865: vector data structure from the array obtained with VecGetArray().
3867: This routine actually zeros out the a pointer.
3869: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3870: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3871: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3872: @*/
3873: PetscErrorCode VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3874: {
3876: void *dummy;
3882: dummy = (void*)(*a + mstart);
3883: PetscFree(dummy);
3884: VecRestoreArray(x,NULL);
3885: return(0);
3886: }
3888: /*@C
3889: VecRestoreArray4dWrite - Restores a vector after VecGetArray3dWrite() has been called.
3891: Logically Collective
3893: Input Parameters:
3894: + x - the vector
3895: . m - first dimension of four dimensional array
3896: . n - second dimension of the four dimensional array
3897: . p - third dimension of the four dimensional array
3898: . q - fourth dimension of the four dimensional array
3899: . mstart - first index you will use in first coordinate direction (often 0)
3900: . nstart - first index in the second coordinate direction (often 0)
3901: . pstart - first index in the third coordinate direction (often 0)
3902: . qstart - first index in the fourth coordinate direction (often 0)
3903: - a - location of pointer to array obtained from VecGetArray4d()
3905: Level: beginner
3907: Notes:
3908: For regular PETSc vectors this routine does not involve any copies. For
3909: any special vectors that do not store local vector data in a contiguous
3910: array, this routine will copy the data back into the underlying
3911: vector data structure from the array obtained with VecGetArray().
3913: This routine actually zeros out the a pointer.
3915: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3916: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3917: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3918: @*/
3919: PetscErrorCode VecRestoreArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3920: {
3922: void *dummy;
3928: dummy = (void*)(*a + mstart);
3929: PetscFree(dummy);
3930: VecRestoreArrayWrite(x,NULL);
3931: return(0);
3932: }
3934: /*@C
3935: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3936: processor's portion of the vector data. You MUST call VecRestoreArray2dRead()
3937: when you no longer need access to the array.
3939: Logically Collective
3941: Input Parameter:
3942: + x - the vector
3943: . m - first dimension of two dimensional array
3944: . n - second dimension of two dimensional array
3945: . mstart - first index you will use in first coordinate direction (often 0)
3946: - nstart - first index in the second coordinate direction (often 0)
3948: Output Parameter:
3949: . a - location to put pointer to the array
3951: Level: developer
3953: Notes:
3954: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
3955: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3956: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3957: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
3959: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3961: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3962: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3963: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3964: @*/
3965: PetscErrorCode VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3966: {
3967: PetscErrorCode ierr;
3968: PetscInt i,N;
3969: const PetscScalar *aa;
3975: VecGetLocalSize(x,&N);
3976: 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);
3977: VecGetArrayRead(x,&aa);
3979: PetscMalloc1(m,a);
3980: for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
3981: *a -= mstart;
3982: return(0);
3983: }
3985: /*@C
3986: VecRestoreArray2dRead - Restores a vector after VecGetArray2dRead() has been called.
3988: Logically Collective
3990: Input Parameters:
3991: + x - the vector
3992: . m - first dimension of two dimensional array
3993: . n - second dimension of the two dimensional array
3994: . mstart - first index you will use in first coordinate direction (often 0)
3995: . nstart - first index in the second coordinate direction (often 0)
3996: - a - location of pointer to array obtained from VecGetArray2d()
3998: Level: developer
4000: Notes:
4001: For regular PETSc vectors this routine does not involve any copies. For
4002: any special vectors that do not store local vector data in a contiguous
4003: array, this routine will copy the data back into the underlying
4004: vector data structure from the array obtained with VecGetArray().
4006: This routine actually zeros out the a pointer.
4008: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
4009: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
4010: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
4011: @*/
4012: PetscErrorCode VecRestoreArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
4013: {
4015: void *dummy;
4021: dummy = (void*)(*a + mstart);
4022: PetscFree(dummy);
4023: VecRestoreArrayRead(x,NULL);
4024: return(0);
4025: }
4027: /*@C
4028: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
4029: processor's portion of the vector data. You MUST call VecRestoreArray1dRead()
4030: when you no longer need access to the array.
4032: Logically Collective
4034: Input Parameter:
4035: + x - the vector
4036: . m - first dimension of two dimensional array
4037: - mstart - first index you will use in first coordinate direction (often 0)
4039: Output Parameter:
4040: . a - location to put pointer to the array
4042: Level: developer
4044: Notes:
4045: For a vector obtained from DMCreateLocalVector() mstart are likely
4046: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
4047: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
4049: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
4051: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
4052: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
4053: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
4054: @*/
4055: PetscErrorCode VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
4056: {
4058: PetscInt N;
4064: VecGetLocalSize(x,&N);
4065: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
4066: VecGetArrayRead(x,(const PetscScalar**)a);
4067: *a -= mstart;
4068: return(0);
4069: }
4071: /*@C
4072: VecRestoreArray1dRead - Restores a vector after VecGetArray1dRead() has been called.
4074: Logically Collective
4076: Input Parameters:
4077: + x - the vector
4078: . m - first dimension of two dimensional array
4079: . mstart - first index you will use in first coordinate direction (often 0)
4080: - a - location of pointer to array obtained from VecGetArray21()
4082: Level: developer
4084: Notes:
4085: For regular PETSc vectors this routine does not involve any copies. For
4086: any special vectors that do not store local vector data in a contiguous
4087: array, this routine will copy the data back into the underlying
4088: vector data structure from the array obtained with VecGetArray1dRead().
4090: This routine actually zeros out the a pointer.
4092: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
4093: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
4094: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
4095: @*/
4096: PetscErrorCode VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
4097: {
4103: VecRestoreArrayRead(x,NULL);
4104: return(0);
4105: }
4108: /*@C
4109: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
4110: processor's portion of the vector data. You MUST call VecRestoreArray3dRead()
4111: when you no longer need access to the array.
4113: Logically Collective
4115: Input Parameter:
4116: + x - the vector
4117: . m - first dimension of three dimensional array
4118: . n - second dimension of three dimensional array
4119: . p - third dimension of three dimensional array
4120: . mstart - first index you will use in first coordinate direction (often 0)
4121: . nstart - first index in the second coordinate direction (often 0)
4122: - pstart - first index in the third coordinate direction (often 0)
4124: Output Parameter:
4125: . a - location to put pointer to the array
4127: Level: developer
4129: Notes:
4130: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
4131: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
4132: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
4133: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3dRead().
4135: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
4137: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
4138: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
4139: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
4140: @*/
4141: PetscErrorCode VecGetArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
4142: {
4143: PetscErrorCode ierr;
4144: PetscInt i,N,j;
4145: const PetscScalar *aa;
4146: PetscScalar **b;
4152: VecGetLocalSize(x,&N);
4153: 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);
4154: VecGetArrayRead(x,&aa);
4156: PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
4157: b = (PetscScalar**)((*a) + m);
4158: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
4159: for (i=0; i<m; i++)
4160: for (j=0; j<n; j++)
4161: b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;
4163: *a -= mstart;
4164: return(0);
4165: }
4167: /*@C
4168: VecRestoreArray3dRead - Restores a vector after VecGetArray3dRead() has been called.
4170: Logically Collective
4172: Input Parameters:
4173: + x - the vector
4174: . m - first dimension of three dimensional array
4175: . n - second dimension of the three dimensional array
4176: . p - third dimension of the three dimensional array
4177: . mstart - first index you will use in first coordinate direction (often 0)
4178: . nstart - first index in the second coordinate direction (often 0)
4179: . pstart - first index in the third coordinate direction (often 0)
4180: - a - location of pointer to array obtained from VecGetArray3dRead()
4182: Level: developer
4184: Notes:
4185: For regular PETSc vectors this routine does not involve any copies. For
4186: any special vectors that do not store local vector data in a contiguous
4187: array, this routine will copy the data back into the underlying
4188: vector data structure from the array obtained with VecGetArray().
4190: This routine actually zeros out the a pointer.
4192: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
4193: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
4194: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
4195: @*/
4196: PetscErrorCode VecRestoreArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
4197: {
4199: void *dummy;
4205: dummy = (void*)(*a + mstart);
4206: PetscFree(dummy);
4207: VecRestoreArrayRead(x,NULL);
4208: return(0);
4209: }
4211: /*@C
4212: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
4213: processor's portion of the vector data. You MUST call VecRestoreArray4dRead()
4214: when you no longer need access to the array.
4216: Logically Collective
4218: Input Parameter:
4219: + x - the vector
4220: . m - first dimension of four dimensional array
4221: . n - second dimension of four dimensional array
4222: . p - third dimension of four dimensional array
4223: . q - fourth dimension of four dimensional array
4224: . mstart - first index you will use in first coordinate direction (often 0)
4225: . nstart - first index in the second coordinate direction (often 0)
4226: . pstart - first index in the third coordinate direction (often 0)
4227: - qstart - first index in the fourth coordinate direction (often 0)
4229: Output Parameter:
4230: . a - location to put pointer to the array
4232: Level: beginner
4234: Notes:
4235: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
4236: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
4237: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
4238: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
4240: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
4242: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
4243: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
4244: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
4245: @*/
4246: PetscErrorCode VecGetArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
4247: {
4248: PetscErrorCode ierr;
4249: PetscInt i,N,j,k;
4250: const PetscScalar *aa;
4251: PetscScalar ***b,**c;
4257: VecGetLocalSize(x,&N);
4258: 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);
4259: VecGetArrayRead(x,&aa);
4261: PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
4262: b = (PetscScalar***)((*a) + m);
4263: c = (PetscScalar**)(b + m*n);
4264: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
4265: for (i=0; i<m; i++)
4266: for (j=0; j<n; j++)
4267: b[i*n+j] = c + i*n*p + j*p - pstart;
4268: for (i=0; i<m; i++)
4269: for (j=0; j<n; j++)
4270: for (k=0; k<p; k++)
4271: c[i*n*p+j*p+k] = (PetscScalar*) aa + i*n*p*q + j*p*q + k*q - qstart;
4272: *a -= mstart;
4273: return(0);
4274: }
4276: /*@C
4277: VecRestoreArray4dRead - Restores a vector after VecGetArray3d() has been called.
4279: Logically Collective
4281: Input Parameters:
4282: + x - the vector
4283: . m - first dimension of four dimensional array
4284: . n - second dimension of the four dimensional array
4285: . p - third dimension of the four dimensional array
4286: . q - fourth dimension of the four dimensional array
4287: . mstart - first index you will use in first coordinate direction (often 0)
4288: . nstart - first index in the second coordinate direction (often 0)
4289: . pstart - first index in the third coordinate direction (often 0)
4290: . qstart - first index in the fourth coordinate direction (often 0)
4291: - a - location of pointer to array obtained from VecGetArray4dRead()
4293: Level: beginner
4295: Notes:
4296: For regular PETSc vectors this routine does not involve any copies. For
4297: any special vectors that do not store local vector data in a contiguous
4298: array, this routine will copy the data back into the underlying
4299: vector data structure from the array obtained with VecGetArray().
4301: This routine actually zeros out the a pointer.
4303: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
4304: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
4305: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
4306: @*/
4307: PetscErrorCode VecRestoreArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
4308: {
4310: void *dummy;
4316: dummy = (void*)(*a + mstart);
4317: PetscFree(dummy);
4318: VecRestoreArrayRead(x,NULL);
4319: return(0);
4320: }
4322: #if defined(PETSC_USE_DEBUG)
4324: /*@
4325: VecLockGet - Gets the current lock status of a vector
4327: Logically Collective on Vec
4329: Input Parameter:
4330: . x - the vector
4332: Output Parameter:
4333: . state - greater than zero indicates the vector is locked for read; less then zero indicates the vector is
4334: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
4336: Level: beginner
4338: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop()
4339: @*/
4340: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
4341: {
4344: *state = x->lock;
4345: return(0);
4346: }
4348: /*@
4349: VecLockReadPush - Pushes a read-only lock on a vector to prevent it from writing
4351: Logically Collective on Vec
4353: Input Parameter:
4354: . x - the vector
4356: Notes:
4357: If this is set then calls to VecGetArray() or VecSetValues() or any other routines that change the vectors values will fail.
4359: The call can be nested, i.e., called multiple times on the same vector, but each VecLockReadPush(x) has to have one matching
4360: VecLockReadPop(x), which removes the latest read-only lock.
4362: Level: beginner
4364: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPop(), VecLockGet()
4365: @*/
4366: PetscErrorCode VecLockReadPush(Vec x)
4367: {
4370: 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");
4371: x->lock++;
4372: return(0);
4373: }
4375: /*@
4376: VecLockReadPop - Pops a read-only lock from a vector
4378: Logically Collective on Vec
4380: Input Parameter:
4381: . x - the vector
4383: Level: beginner
4385: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockGet()
4386: @*/
4387: PetscErrorCode VecLockReadPop(Vec x)
4388: {
4391: x->lock--;
4392: if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector has been unlocked from read-only access too many times");
4393: return(0);
4394: }
4396: /*@C
4397: VecLockWriteSet_Private - Lock or unlock a vector for exclusive read/write access
4399: Logically Collective on Vec
4401: Input Parameter:
4402: + x - the vector
4403: - flg - PETSC_TRUE to lock the vector for writing; PETSC_FALSE to unlock it.
4405: Notes:
4406: The function is usefull in split-phase computations, which usually have a begin phase and an end phase.
4407: One can call VecLockWriteSet_Private(x,PETSC_TRUE) in the begin phase to lock a vector for exclusive
4408: access, and call VecLockWriteSet_Private(x,PETSC_FALSE) in the end phase to unlock the vector from exclusive
4409: access. In this way, one is ensured no other operations can access the vector in between. The code may like
4412: VecGetArray(x,&xdata); // begin phase
4413: VecLockWriteSet_Private(v,PETSC_TRUE);
4415: Other operations, which can not acceess x anymore (they can access xdata, of course)
4417: VecRestoreArray(x,&vdata); // end phase
4418: VecLockWriteSet_Private(v,PETSC_FALSE);
4420: The call can not be nested on the same vector, in other words, one can not call VecLockWriteSet_Private(x,PETSC_TRUE)
4421: again before calling VecLockWriteSet_Private(v,PETSC_FALSE).
4423: Level: beginner
4425: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop(), VecLockGet()
4426: @*/
4427: PetscErrorCode VecLockWriteSet_Private(Vec x,PetscBool flg)
4428: {
4431: if (flg) {
4432: 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");
4433: 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");
4434: else x->lock = -1;
4435: } else {
4436: 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");
4437: x->lock = 0;
4438: }
4439: return(0);
4440: }
4442: /*@
4443: VecLockPush - Pushes a read-only lock on a vector to prevent it from writing
4445: Level: deprecated
4447: .seealso: VecLockReadPush()
4448: @*/
4449: PetscErrorCode VecLockPush(Vec x)
4450: {
4453: VecLockReadPush(x);
4454: return(0);
4455: }
4457: /*@
4458: VecLockPop - Pops a read-only lock from a vector
4460: Level: deprecated
4462: .seealso: VecLockReadPop()
4463: @*/
4464: PetscErrorCode VecLockPop(Vec x)
4465: {
4468: VecLockReadPop(x);
4469: return(0);
4470: }
4472: #endif