Actual source code: rvector.c
petsc-master 2020-09-19
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/vecimpl.h>
6: #if defined(PETSC_HAVE_CUDA)
7: #include <../src/vec/vec/impls/dvecimpl.h>
8: #include <petsc/private/cudavecimpl.h>
9: #endif
10: static PetscInt VecGetSubVectorSavedStateId = -1;
12: PETSC_EXTERN PetscErrorCode VecValidValues(Vec vec,PetscInt argnum,PetscBool begin)
13: {
14: #if defined(PETSC_USE_DEBUG)
15: PetscErrorCode ierr;
16: PetscInt n,i;
17: const PetscScalar *x;
20: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL)
21: if ((vec->petscnative || vec->ops->getarray) && (vec->offloadmask & PETSC_OFFLOAD_CPU)) {
22: #else
23: if (vec->petscnative || vec->ops->getarray) {
24: #endif
25: VecGetLocalSize(vec,&n);
26: VecGetArrayRead(vec,&x);
27: for (i=0; i<n; i++) {
28: if (begin) {
29: 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);
30: } else {
31: if (PetscIsInfOrNanScalar(x[i])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FP,"Vec entry at local location %D is not-a-number or infinite at end of function: Parameter number %D",i,argnum);
32: }
33: }
34: VecRestoreArrayRead(vec,&x);
35: }
36: #else
38: #endif
39: return(0);
40: }
42: /*@
43: VecMaxPointwiseDivide - Computes the maximum of the componentwise division max = max_i abs(x_i/y_i).
45: Logically Collective on Vec
47: Input Parameters:
48: . x, y - the vectors
50: Output Parameter:
51: . max - the result
53: Level: advanced
55: Notes:
56: x and y may be the same vector
57: if a particular y_i is zero, it is treated as 1 in the above formula
59: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs()
60: @*/
61: PetscErrorCode VecMaxPointwiseDivide(Vec x,Vec y,PetscReal *max)
62: {
72: VecCheckSameSize(x,1,y,2);
73: (*x->ops->maxpointwisedivide)(x,y,max);
74: return(0);
75: }
77: /*@
78: VecDot - Computes the vector dot product.
80: Collective on Vec
82: Input Parameters:
83: . x, y - the vectors
85: Output Parameter:
86: . val - the dot product
88: Performance Issues:
89: $ per-processor memory bandwidth
90: $ interprocessor latency
91: $ work load inbalance that causes certain processes to arrive much earlier than others
93: Notes for Users of Complex Numbers:
94: For complex vectors, VecDot() computes
95: $ val = (x,y) = y^H x,
96: where y^H denotes the conjugate transpose of y. Note that this corresponds to the usual "mathematicians" complex
97: inner product where the SECOND argument gets the complex conjugate. Since the BLASdot() complex conjugates the first
98: first argument we call the BLASdot() with the arguments reversed.
100: Use VecTDot() for the indefinite form
101: $ val = (x,y) = y^T x,
102: where y^T denotes the transpose of y.
104: Level: intermediate
107: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDotRealPart()
108: @*/
109: PetscErrorCode VecDot(Vec x,Vec y,PetscScalar *val)
110: {
120: VecCheckSameSize(x,1,y,2);
122: PetscLogEventBegin(VEC_Dot,x,y,0,0);
123: (*x->ops->dot)(x,y,val);
124: PetscLogEventEnd(VEC_Dot,x,y,0,0);
125: return(0);
126: }
128: /*@
129: VecDotRealPart - Computes the real part of the vector dot product.
131: Collective on Vec
133: Input Parameters:
134: . x, y - the vectors
136: Output Parameter:
137: . val - the real part of the dot product;
139: Performance Issues:
140: $ per-processor memory bandwidth
141: $ interprocessor latency
142: $ work load inbalance that causes certain processes to arrive much earlier than others
144: Notes for Users of Complex Numbers:
145: See VecDot() for more details on the definition of the dot product for complex numbers
147: For real numbers this returns the same value as VecDot()
149: 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
150: the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)
152: Developer Note: This is not currently optimized to compute only the real part of the dot product.
154: Level: intermediate
157: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDot(), VecDotNorm2()
158: @*/
159: PetscErrorCode VecDotRealPart(Vec x,Vec y,PetscReal *val)
160: {
162: PetscScalar fdot;
165: VecDot(x,y,&fdot);
166: *val = PetscRealPart(fdot);
167: return(0);
168: }
170: /*@
171: VecNorm - Computes the vector norm.
173: Collective on Vec
175: Input Parameters:
176: + x - the vector
177: - type - one of NORM_1, NORM_2, NORM_INFINITY. Also available
178: NORM_1_AND_2, which computes both norms and stores them
179: in a two element array.
181: Output Parameter:
182: . val - the norm
184: Notes:
185: $ NORM_1 denotes sum_i |x_i|
186: $ NORM_2 denotes sqrt(sum_i |x_i|^2)
187: $ NORM_INFINITY denotes max_i |x_i|
189: For complex numbers NORM_1 will return the traditional 1 norm of the 2 norm of the complex numbers; that is the 1
190: norm of the absolutely values of the complex entries. In PETSc 3.6 and earlier releases it returned the 1 norm of
191: the 1 norm of the complex entries (what is returned by the BLAS routine asum()). Both are valid norms but most
192: people expect the former.
194: Level: intermediate
196: Performance Issues:
197: $ per-processor memory bandwidth
198: $ interprocessor latency
199: $ work load inbalance that causes certain processes to arrive much earlier than others
202: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNormAvailable(),
203: VecNormBegin(), VecNormEnd()
205: @*/
207: PetscErrorCode VecNorm(Vec x,NormType type,PetscReal *val)
208: {
209: PetscBool flg;
217: /*
218: * Cached data?
219: */
220: if (type!=NORM_1_AND_2) {
221: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,flg);
222: if (flg) return(0);
223: }
224: PetscLogEventBegin(VEC_Norm,x,0,0,0);
225: (*x->ops->norm)(x,type,val);
226: PetscLogEventEnd(VEC_Norm,x,0,0,0);
227: if (type!=NORM_1_AND_2) {
228: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[type],*val);
229: }
230: return(0);
231: }
233: /*@
234: VecNormAvailable - Returns the vector norm if it is already known.
236: Not Collective
238: Input Parameters:
239: + x - the vector
240: - type - one of NORM_1, NORM_2, NORM_INFINITY. Also available
241: NORM_1_AND_2, which computes both norms and stores them
242: in a two element array.
244: Output Parameter:
245: + available - PETSC_TRUE if the val returned is valid
246: - val - the norm
248: Notes:
249: $ NORM_1 denotes sum_i |x_i|
250: $ NORM_2 denotes sqrt(sum_i (x_i)^2)
251: $ NORM_INFINITY denotes max_i |x_i|
253: Level: intermediate
255: Performance Issues:
256: $ per-processor memory bandwidth
257: $ interprocessor latency
258: $ work load inbalance that causes certain processes to arrive much earlier than others
260: Compile Option:
261: PETSC_HAVE_SLOW_BLAS_NORM2 will cause a C (loop unrolled) version of the norm to be used, rather
262: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
263: (as opposed to vendor provided) because the FORTRAN BLAS NRM2() routine is very slow.
266: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNorm()
267: VecNormBegin(), VecNormEnd()
269: @*/
270: PetscErrorCode VecNormAvailable(Vec x,NormType type,PetscBool *available,PetscReal *val)
271: {
279: *available = PETSC_FALSE;
280: if (type!=NORM_1_AND_2) {
281: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,*available);
282: }
283: return(0);
284: }
286: /*@
287: VecNormalize - Normalizes a vector by 2-norm.
289: Collective on Vec
291: Input Parameters:
292: + x - the vector
294: Output Parameter:
295: . x - the normalized vector
296: - val - the vector norm before normalization
298: Level: intermediate
301: @*/
302: PetscErrorCode VecNormalize(Vec x,PetscReal *val)
303: {
305: PetscReal norm;
310: PetscLogEventBegin(VEC_Normalize,x,0,0,0);
311: VecNorm(x,NORM_2,&norm);
312: if (norm == 0.0) {
313: PetscInfo(x,"Vector of zero norm can not be normalized; Returning only the zero norm\n");
314: } else if (norm != 1.0) {
315: PetscScalar tmp = 1.0/norm;
316: VecScale(x,tmp);
317: }
318: if (val) *val = norm;
319: PetscLogEventEnd(VEC_Normalize,x,0,0,0);
320: return(0);
321: }
323: /*@C
324: VecMax - Determines the vector component with maximum real part and its location.
326: Collective on Vec
328: Input Parameter:
329: . x - the vector
331: Output Parameters:
332: + p - the location of val (pass NULL if you don't want this)
333: - val - the maximum component
335: Notes:
336: Returns the value PETSC_MIN_REAL and p = -1 if the vector is of length 0.
338: Returns the smallest index with the maximum value
339: Level: intermediate
342: .seealso: VecNorm(), VecMin()
343: @*/
344: PetscErrorCode VecMax(Vec x,PetscInt *p,PetscReal *val)
345: {
352: PetscLogEventBegin(VEC_Max,x,0,0,0);
353: (*x->ops->max)(x,p,val);
354: PetscLogEventEnd(VEC_Max,x,0,0,0);
355: return(0);
356: }
358: /*@C
359: VecMin - Determines the vector component with minimum real part and its location.
361: Collective on Vec
363: Input Parameters:
364: . x - the vector
366: Output Parameter:
367: + p - the location of val (pass NULL if you don't want this location)
368: - val - the minimum component
370: Level: intermediate
372: Notes:
373: Returns the value PETSC_MAX_REAL and p = -1 if the vector is of length 0.
375: This returns the smallest index with the minumum value
378: .seealso: VecMax()
379: @*/
380: PetscErrorCode VecMin(Vec x,PetscInt *p,PetscReal *val)
381: {
388: PetscLogEventBegin(VEC_Min,x,0,0,0);
389: (*x->ops->min)(x,p,val);
390: PetscLogEventEnd(VEC_Min,x,0,0,0);
391: return(0);
392: }
394: /*@
395: VecTDot - Computes an indefinite vector dot product. That is, this
396: routine does NOT use the complex conjugate.
398: Collective on Vec
400: Input Parameters:
401: . x, y - the vectors
403: Output Parameter:
404: . val - the dot product
406: Notes for Users of Complex Numbers:
407: For complex vectors, VecTDot() computes the indefinite form
408: $ val = (x,y) = y^T x,
409: where y^T denotes the transpose of y.
411: Use VecDot() for the inner product
412: $ val = (x,y) = y^H x,
413: where y^H denotes the conjugate transpose of y.
415: Level: intermediate
417: .seealso: VecDot(), VecMTDot()
418: @*/
419: PetscErrorCode VecTDot(Vec x,Vec y,PetscScalar *val)
420: {
430: VecCheckSameSize(x,1,y,2);
432: PetscLogEventBegin(VEC_TDot,x,y,0,0);
433: (*x->ops->tdot)(x,y,val);
434: PetscLogEventEnd(VEC_TDot,x,y,0,0);
435: return(0);
436: }
438: /*@
439: VecScale - Scales a vector.
441: Not collective on Vec
443: Input Parameters:
444: + x - the vector
445: - alpha - the scalar
447: Output Parameter:
448: . x - the scaled vector
450: Note:
451: For a vector with n components, VecScale() computes
452: $ x[i] = alpha * x[i], for i=1,...,n.
454: Level: intermediate
457: @*/
458: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
459: {
460: PetscReal norms[4] = {0.0,0.0,0.0, 0.0};
461: PetscBool flgs[4];
463: PetscInt i;
468: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
469: PetscLogEventBegin(VEC_Scale,x,0,0,0);
470: if (alpha != (PetscScalar)1.0) {
471: VecSetErrorIfLocked(x,1);
472: /* get current stashed norms */
473: for (i=0; i<4; i++) {
474: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
475: }
476: (*x->ops->scale)(x,alpha);
477: PetscObjectStateIncrease((PetscObject)x);
478: /* put the scaled stashed norms back into the Vec */
479: for (i=0; i<4; i++) {
480: if (flgs[i]) {
481: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],PetscAbsScalar(alpha)*norms[i]);
482: }
483: }
484: }
485: PetscLogEventEnd(VEC_Scale,x,0,0,0);
486: return(0);
487: }
489: /*@
490: VecSet - Sets all components of a vector to a single scalar value.
492: Logically Collective on Vec
494: Input Parameters:
495: + x - the vector
496: - alpha - the scalar
498: Output Parameter:
499: . x - the vector
501: Note:
502: For a vector of dimension n, VecSet() computes
503: $ x[i] = alpha, for i=1,...,n,
504: so that all vector entries then equal the identical
505: scalar value, alpha. Use the more general routine
506: VecSetValues() to set different vector entries.
508: You CANNOT call this after you have called VecSetValues() but before you call
509: VecAssemblyBegin/End().
511: Level: beginner
513: .seealso VecSetValues(), VecSetValuesBlocked(), VecSetRandom()
515: @*/
516: PetscErrorCode VecSet(Vec x,PetscScalar alpha)
517: {
518: PetscReal val;
524: 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()");
526: VecSetErrorIfLocked(x,1);
528: PetscLogEventBegin(VEC_Set,x,0,0,0);
529: (*x->ops->set)(x,alpha);
530: PetscLogEventEnd(VEC_Set,x,0,0,0);
531: PetscObjectStateIncrease((PetscObject)x);
533: /* norms can be simply set (if |alpha|*N not too large) */
534: val = PetscAbsScalar(alpha);
535: if (x->map->N == 0) {
536: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],0.0l);
537: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],0.0);
538: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],0.0);
539: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],0.0);
540: } else if (val > PETSC_MAX_REAL/x->map->N) {
541: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
542: } else {
543: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],x->map->N * val);
544: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
545: val = PetscSqrtReal((PetscReal)x->map->N) * val;
546: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],val);
547: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],val);
548: }
549: return(0);
550: }
553: /*@
554: VecAXPY - Computes y = alpha x + y.
556: Logically Collective on Vec
558: Input Parameters:
559: + alpha - the scalar
560: - x, y - the vectors
562: Output Parameter:
563: . y - output vector
565: Level: intermediate
567: Notes:
568: x and y MUST be different vectors
569: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
571: $ VecAXPY(y,alpha,x) y = alpha x + y
572: $ VecAYPX(y,beta,x) y = x + beta y
573: $ VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
574: $ VecWAXPY(w,alpha,x,y) w = alpha x + y
575: $ VecAXPBYPCZ(w,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
576: $ VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
579: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPBYPCZ(), VecAXPBY()
580: @*/
581: PetscErrorCode VecAXPY(Vec y,PetscScalar alpha,Vec x)
582: {
591: VecCheckSameSize(x,1,y,3);
592: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
594: if (alpha == (PetscScalar)0.0) return(0);
595: VecSetErrorIfLocked(y,1);
597: VecLockReadPush(x);
598: PetscLogEventBegin(VEC_AXPY,x,y,0,0);
599: (*y->ops->axpy)(y,alpha,x);
600: PetscLogEventEnd(VEC_AXPY,x,y,0,0);
601: VecLockReadPop(x);
602: PetscObjectStateIncrease((PetscObject)y);
603: return(0);
604: }
606: /*@
607: VecAXPBY - Computes y = alpha x + beta y.
609: Logically Collective on Vec
611: Input Parameters:
612: + alpha,beta - the scalars
613: - x, y - the vectors
615: Output Parameter:
616: . y - output vector
618: Level: intermediate
620: Notes:
621: x and y MUST be different vectors
622: The implementation is optimized for alpha and/or beta values of 0.0 and 1.0
625: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ()
626: @*/
627: PetscErrorCode VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
628: {
637: VecCheckSameSize(y,1,x,4);
638: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
641: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) return(0);
642: VecSetErrorIfLocked(y,1);
643: PetscLogEventBegin(VEC_AXPY,x,y,0,0);
644: (*y->ops->axpby)(y,alpha,beta,x);
645: PetscLogEventEnd(VEC_AXPY,x,y,0,0);
646: PetscObjectStateIncrease((PetscObject)y);
647: return(0);
648: }
650: /*@
651: VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z
653: Logically Collective on Vec
655: Input Parameters:
656: + alpha,beta, gamma - the scalars
657: - x, y, z - the vectors
659: Output Parameter:
660: . z - output vector
662: Level: intermediate
664: Notes:
665: x, y and z must be different vectors
666: The implementation is optimized for alpha of 1.0 and gamma of 1.0 or 0.0
669: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBY()
670: @*/
671: PetscErrorCode VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
672: {
684: VecCheckSameSize(x,1,y,5);
685: VecCheckSameSize(x,1,z,6);
686: if (x == y || x == z) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
687: if (y == z) SETERRQ(PetscObjectComm((PetscObject)y),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
691: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) return(0);
692: VecSetErrorIfLocked(z,1);
694: PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
695: (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
696: PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
697: PetscObjectStateIncrease((PetscObject)z);
698: return(0);
699: }
701: /*@
702: VecAYPX - Computes y = x + beta y.
704: Logically Collective on Vec
706: Input Parameters:
707: + beta - the scalar
708: - x, y - the vectors
710: Output Parameter:
711: . y - output vector
713: Level: intermediate
715: Notes:
716: x and y MUST be different vectors
717: The implementation is optimized for beta of -1.0, 0.0, and 1.0
720: .seealso: VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
721: @*/
722: PetscErrorCode VecAYPX(Vec y,PetscScalar beta,Vec x)
723: {
732: VecCheckSameSize(x,1,y,3);
733: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
735: VecSetErrorIfLocked(y,1);
737: PetscLogEventBegin(VEC_AYPX,x,y,0,0);
738: (*y->ops->aypx)(y,beta,x);
739: PetscLogEventEnd(VEC_AYPX,x,y,0,0);
740: PetscObjectStateIncrease((PetscObject)y);
741: return(0);
742: }
745: /*@
746: VecWAXPY - Computes w = alpha x + y.
748: Logically Collective on Vec
750: Input Parameters:
751: + alpha - the scalar
752: - x, y - the vectors
754: Output Parameter:
755: . w - the result
757: Level: intermediate
759: Notes:
760: w cannot be either x or y, but x and y can be the same
761: The implementation is optimzed for alpha of -1.0, 0.0, and 1.0
764: .seealso: VecAXPY(), VecAYPX(), VecAXPBY(), VecMAXPY(), VecAXPBYPCZ()
765: @*/
766: PetscErrorCode VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
767: {
779: VecCheckSameSize(x,3,y,4);
780: VecCheckSameSize(x,3,w,1);
781: if (w == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector y, suggest VecAXPY()");
782: if (w == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector x, suggest VecAYPX()");
784: VecSetErrorIfLocked(w,1);
786: PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
787: (*w->ops->waxpy)(w,alpha,x,y);
788: PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
789: PetscObjectStateIncrease((PetscObject)w);
790: return(0);
791: }
794: /*@C
795: VecSetValues - Inserts or adds values into certain locations of a vector.
797: Not Collective
799: Input Parameters:
800: + x - vector to insert in
801: . ni - number of elements to add
802: . ix - indices where to add
803: . y - array of values
804: - iora - either INSERT_VALUES or ADD_VALUES, where
805: ADD_VALUES adds values to any existing entries, and
806: INSERT_VALUES replaces existing entries with new values
808: Notes:
809: VecSetValues() sets x[ix[i]] = y[i], for i=0,...,ni-1.
811: Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
812: options cannot be mixed without intervening calls to the assembly
813: routines.
815: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
816: MUST be called after all calls to VecSetValues() have been completed.
818: VecSetValues() uses 0-based indices in Fortran as well as in C.
820: If you call VecSetOption(x, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE),
821: negative indices may be passed in ix. These rows are
822: simply ignored. This allows easily inserting element load matrices
823: with homogeneous Dirchlet boundary conditions that you don't want represented
824: in the vector.
826: Level: beginner
828: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
829: VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
830: @*/
831: PetscErrorCode VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
832: {
837: if (!ni) return(0);
842: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
843: (*x->ops->setvalues)(x,ni,ix,y,iora);
844: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
845: PetscObjectStateIncrease((PetscObject)x);
846: return(0);
847: }
849: /*@C
850: VecGetValues - Gets values from certain locations of a vector. Currently
851: can only get values on the same processor
853: Not Collective
855: Input Parameters:
856: + x - vector to get values from
857: . ni - number of elements to get
858: - ix - indices where to get them from (in global 1d numbering)
860: Output Parameter:
861: . y - array of values
863: Notes:
864: The user provides the allocated array y; it is NOT allocated in this routine
866: VecGetValues() gets y[i] = x[ix[i]], for i=0,...,ni-1.
868: VecAssemblyBegin() and VecAssemblyEnd() MUST be called before calling this
870: VecGetValues() uses 0-based indices in Fortran as well as in C.
872: If you call VecSetOption(x, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE),
873: negative indices may be passed in ix. These rows are
874: simply ignored.
876: Level: beginner
878: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues()
879: @*/
880: PetscErrorCode VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
881: {
886: if (!ni) return(0);
890: (*x->ops->getvalues)(x,ni,ix,y);
891: return(0);
892: }
894: /*@C
895: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
897: Not Collective
899: Input Parameters:
900: + x - vector to insert in
901: . ni - number of blocks to add
902: . ix - indices where to add in block count, rather than element count
903: . y - array of values
904: - iora - either INSERT_VALUES or ADD_VALUES, where
905: ADD_VALUES adds values to any existing entries, and
906: INSERT_VALUES replaces existing entries with new values
908: Notes:
909: VecSetValuesBlocked() sets x[bs*ix[i]+j] = y[bs*i+j],
910: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
912: Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
913: options cannot be mixed without intervening calls to the assembly
914: routines.
916: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
917: MUST be called after all calls to VecSetValuesBlocked() have been completed.
919: VecSetValuesBlocked() uses 0-based indices in Fortran as well as in C.
921: Negative indices may be passed in ix, these rows are
922: simply ignored. This allows easily inserting element load matrices
923: with homogeneous Dirchlet boundary conditions that you don't want represented
924: in the vector.
926: Level: intermediate
928: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
929: VecSetValues()
930: @*/
931: PetscErrorCode VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
932: {
937: if (!ni) return(0);
942: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
943: (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
944: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
945: PetscObjectStateIncrease((PetscObject)x);
946: return(0);
947: }
950: /*@C
951: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
952: using a local ordering of the nodes.
954: Not Collective
956: Input Parameters:
957: + x - vector to insert in
958: . ni - number of elements to add
959: . ix - indices where to add
960: . y - array of values
961: - iora - either INSERT_VALUES or ADD_VALUES, where
962: ADD_VALUES adds values to any existing entries, and
963: INSERT_VALUES replaces existing entries with new values
965: Level: intermediate
967: Notes:
968: VecSetValuesLocal() sets x[ix[i]] = y[i], for i=0,...,ni-1.
970: Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
971: options cannot be mixed without intervening calls to the assembly
972: routines.
974: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
975: MUST be called after all calls to VecSetValuesLocal() have been completed.
977: VecSetValuesLocal() uses 0-based indices in Fortran as well as in C.
979: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetLocalToGlobalMapping(),
980: VecSetValuesBlockedLocal()
981: @*/
982: PetscErrorCode VecSetValuesLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
983: {
985: PetscInt lixp[128],*lix = lixp;
989: if (!ni) return(0);
994: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
995: if (!x->ops->setvalueslocal) {
996: if (!x->map->mapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMapping()");
997: if (ni > 128) {
998: PetscMalloc1(ni,&lix);
999: }
1000: ISLocalToGlobalMappingApply(x->map->mapping,ni,(PetscInt*)ix,lix);
1001: (*x->ops->setvalues)(x,ni,lix,y,iora);
1002: if (ni > 128) {
1003: PetscFree(lix);
1004: }
1005: } else {
1006: (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
1007: }
1008: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1009: PetscObjectStateIncrease((PetscObject)x);
1010: return(0);
1011: }
1013: /*@
1014: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1015: using a local ordering of the nodes.
1017: Not Collective
1019: Input Parameters:
1020: + x - vector to insert in
1021: . ni - number of blocks to add
1022: . ix - indices where to add in block count, not element count
1023: . y - array of values
1024: - iora - either INSERT_VALUES or ADD_VALUES, where
1025: ADD_VALUES adds values to any existing entries, and
1026: INSERT_VALUES replaces existing entries with new values
1028: Level: intermediate
1030: Notes:
1031: VecSetValuesBlockedLocal() sets x[bs*ix[i]+j] = y[bs*i+j],
1032: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with VecSetBlockSize().
1034: Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1035: options cannot be mixed without intervening calls to the assembly
1036: routines.
1038: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
1039: MUST be called after all calls to VecSetValuesBlockedLocal() have been completed.
1041: VecSetValuesBlockedLocal() uses 0-based indices in Fortran as well as in C.
1044: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
1045: VecSetLocalToGlobalMapping()
1046: @*/
1047: PetscErrorCode VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1048: {
1050: PetscInt lixp[128],*lix = lixp;
1054: if (!ni) return(0);
1058: if (ni > 128) {
1059: PetscMalloc1(ni,&lix);
1060: }
1062: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1063: ISLocalToGlobalMappingApplyBlock(x->map->mapping,ni,(PetscInt*)ix,lix);
1064: (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1065: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1066: if (ni > 128) {
1067: PetscFree(lix);
1068: }
1069: PetscObjectStateIncrease((PetscObject)x);
1070: return(0);
1071: }
1073: /*@
1074: VecMTDot - Computes indefinite vector multiple dot products.
1075: That is, it does NOT use the complex conjugate.
1077: Collective on Vec
1079: Input Parameters:
1080: + x - one vector
1081: . nv - number of vectors
1082: - y - array of vectors. Note that vectors are pointers
1084: Output Parameter:
1085: . val - array of the dot products
1087: Notes for Users of Complex Numbers:
1088: For complex vectors, VecMTDot() computes the indefinite form
1089: $ val = (x,y) = y^T x,
1090: where y^T denotes the transpose of y.
1092: Use VecMDot() for the inner product
1093: $ val = (x,y) = y^H x,
1094: where y^H denotes the conjugate transpose of y.
1096: Level: intermediate
1099: .seealso: VecMDot(), VecTDot()
1100: @*/
1101: PetscErrorCode VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1102: {
1108: if (!nv) return(0);
1115: VecCheckSameSize(x,1,*y,3);
1117: PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1118: (*x->ops->mtdot)(x,nv,y,val);
1119: PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1120: return(0);
1121: }
1123: /*@
1124: VecMDot - Computes vector multiple dot products.
1126: Collective on Vec
1128: Input Parameters:
1129: + x - one vector
1130: . nv - number of vectors
1131: - y - array of vectors.
1133: Output Parameter:
1134: . val - array of the dot products (does not allocate the array)
1136: Notes for Users of Complex Numbers:
1137: For complex vectors, VecMDot() computes
1138: $ val = (x,y) = y^H x,
1139: where y^H denotes the conjugate transpose of y.
1141: Use VecMTDot() for the indefinite form
1142: $ val = (x,y) = y^T x,
1143: where y^T denotes the transpose of y.
1145: Level: intermediate
1148: .seealso: VecMTDot(), VecDot()
1149: @*/
1150: PetscErrorCode VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1151: {
1157: if (!nv) return(0);
1158: if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1165: VecCheckSameSize(x,1,*y,3);
1167: PetscLogEventBegin(VEC_MDot,x,*y,0,0);
1168: (*x->ops->mdot)(x,nv,y,val);
1169: PetscLogEventEnd(VEC_MDot,x,*y,0,0);
1170: return(0);
1171: }
1173: /*@
1174: VecMAXPY - Computes y = y + sum alpha[i] x[i]
1176: Logically Collective on Vec
1178: Input Parameters:
1179: + nv - number of scalars and x-vectors
1180: . alpha - array of scalars
1181: . y - one vector
1182: - x - array of vectors
1184: Level: intermediate
1186: Notes:
1187: y cannot be any of the x vectors
1189: .seealso: VecAYPX(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
1190: @*/
1191: PetscErrorCode VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1192: {
1194: PetscInt i;
1195: PetscBool nonzero;
1200: if (!nv) return(0);
1201: if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1208: VecCheckSameSize(y,1,*x,4);
1210: for (i=0, nonzero = PETSC_FALSE; i<nv && !nonzero; i++) nonzero = (PetscBool)(nonzero || alpha[i] != (PetscScalar)0.0);
1211: if (!nonzero) return(0);
1212: VecSetErrorIfLocked(y,1);
1213: PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1214: (*y->ops->maxpy)(y,nv,alpha,x);
1215: PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1216: PetscObjectStateIncrease((PetscObject)y);
1217: return(0);
1218: }
1220: /*@
1221: VecGetSubVector - Gets a vector representing part of another vector
1223: Collective on IS
1225: Input Arguments:
1226: + X - vector from which to extract a subvector
1227: - is - index set representing portion of X to extract
1229: Output Arguments:
1230: . Y - subvector corresponding to is
1232: Level: advanced
1234: Notes:
1235: The subvector Y should be returned with VecRestoreSubVector().
1237: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1238: modifying the subvector. Other non-overlapping subvectors can still be obtained from X using this function.
1240: .seealso: MatCreateSubMatrix()
1241: @*/
1242: PetscErrorCode VecGetSubVector(Vec X,IS is,Vec *Y)
1243: {
1244: PetscErrorCode ierr;
1245: Vec Z;
1251: if (X->ops->getsubvector) {
1252: (*X->ops->getsubvector)(X,is,&Z);
1253: } else { /* Default implementation currently does no caching */
1254: PetscInt gstart,gend,start;
1255: PetscBool contiguous,gcontiguous;
1257: VecGetOwnershipRange(X,&gstart,&gend);
1258: ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1259: MPIU_Allreduce(&contiguous,&gcontiguous,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1260: if (gcontiguous) { /* We can do a no-copy implementation */
1261: const PetscScalar *x;
1262: PetscInt n,N,bs;
1263: PetscInt state = 0;
1264: #if defined(PETSC_HAVE_CUDA)
1265: PetscBool iscuda;
1266: #endif
1268: ISGetSize(is,&N);
1269: ISGetLocalSize(is,&n);
1270: VecGetBlockSize(X,&bs);
1271: if (n%bs || bs == 1) bs = -1; /* Do not decide block size if we do not have to */
1272: #if defined(PETSC_HAVE_CUDA)
1273: PetscObjectTypeCompareAny((PetscObject)X,&iscuda,VECSEQCUDA,VECMPICUDA,"");
1274: if (iscuda) {
1275: const PetscScalar *x_d;
1276: PetscMPIInt size;
1277: PetscOffloadMask flg;
1279: VecCUDAGetArrays_Private(X,&x,&x_d,&flg);
1280: if (flg == PETSC_OFFLOAD_UNALLOCATED) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for PETSC_OFFLOAD_UNALLOCATED");
1281: if (n && !x && !x_d) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Missing vector data");
1282: if (x) x += start;
1283: if (x_d) x_d += start;
1284: MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1285: if (size == 1) {
1286: VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X),bs,n,x,x_d,&Z);
1287: } else {
1288: VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X),bs,n,N,x,x_d,&Z);
1289: }
1290: Z->offloadmask = flg;
1291: } else {
1292: #else
1293: {
1294: #endif
1295: VecGetArrayRead(X,&x);
1296: VecCreate(PetscObjectComm((PetscObject)X),&Z);
1297: VecSetType(Z,((PetscObject)X)->type_name);
1298: VecSetSizes(Z,n,N);
1299: VecSetBlockSize(Z,bs);
1300: VecPlaceArray(Z,x+start);
1301: VecRestoreArrayRead(X,&x);
1302: }
1304: /* this is relevant only in debug mode */
1305: VecLockGet(X,&state);
1306: if (state) {
1307: VecLockReadPush(Z);
1308: }
1309: Z->ops->placearray = NULL;
1310: Z->ops->replacearray = NULL;
1311: } else { /* Have to create a scatter and do a copy */
1312: VecScatter scatter;
1313: PetscInt n,N;
1315: ISGetLocalSize(is,&n);
1316: ISGetSize(is,&N);
1317: VecCreate(PetscObjectComm((PetscObject)is),&Z);
1318: VecSetSizes(Z,n,N);
1319: VecSetType(Z,((PetscObject)X)->type_name);
1320: VecScatterCreate(X,is,Z,NULL,&scatter);
1321: VecScatterBegin(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1322: VecScatterEnd(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1323: PetscObjectCompose((PetscObject)Z,"VecGetSubVector_Scatter",(PetscObject)scatter);
1324: VecScatterDestroy(&scatter);
1325: }
1326: }
1327: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1328: if (VecGetSubVectorSavedStateId < 0) {PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);}
1329: PetscObjectComposedDataSetInt((PetscObject)Z,VecGetSubVectorSavedStateId,1);
1330: *Y = Z;
1331: return(0);
1332: }
1334: /*@
1335: VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()
1337: Collective on IS
1339: Input Arguments:
1340: + X - vector from which subvector was obtained
1341: . is - index set representing the subset of X
1342: - Y - subvector being restored
1344: Level: advanced
1346: .seealso: VecGetSubVector()
1347: @*/
1348: PetscErrorCode VecRestoreSubVector(Vec X,IS is,Vec *Y)
1349: {
1357: if (X->ops->restoresubvector) {
1358: (*X->ops->restoresubvector)(X,is,Y);
1359: } else {
1360: PETSC_UNUSED PetscObjectState dummystate = 0;
1361: PetscBool valid;
1363: PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,dummystate,valid);
1364: if (!valid) {
1365: VecScatter scatter;
1366: PetscInt state;
1368: VecLockGet(X,&state);
1369: if (state != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vec X is locked for read-only or read/write access");
1371: PetscObjectQuery((PetscObject)*Y,"VecGetSubVector_Scatter",(PetscObject*)&scatter);
1372: if (scatter) {
1373: VecScatterBegin(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1374: VecScatterEnd(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1375: } else {
1376: #if defined(PETSC_HAVE_CUDA)
1377: PetscBool iscuda;
1379: PetscObjectTypeCompareAny((PetscObject)*Y,&iscuda,VECSEQCUDA,VECMPICUDA,"");
1380: if (iscuda) {
1381: PetscOffloadMask ymask = (*Y)->offloadmask;
1383: /* The offloadmask of X dictates where to move memory
1384: If X GPU data is valid, then move Y data on GPU if needed
1385: Otherwise, move back to the CPU */
1386: switch (X->offloadmask) {
1387: case PETSC_OFFLOAD_BOTH:
1388: if (ymask == PETSC_OFFLOAD_CPU) {
1389: VecCUDAResetArray(*Y);
1390: } else if (ymask == PETSC_OFFLOAD_GPU) {
1391: X->offloadmask = PETSC_OFFLOAD_GPU;
1392: }
1393: break;
1394: case PETSC_OFFLOAD_GPU:
1395: if (ymask == PETSC_OFFLOAD_CPU) {
1396: VecCUDAResetArray(*Y);
1397: }
1398: break;
1399: case PETSC_OFFLOAD_CPU:
1400: if (ymask == PETSC_OFFLOAD_GPU) {
1401: VecResetArray(*Y);
1402: }
1403: break;
1404: case PETSC_OFFLOAD_UNALLOCATED:
1405: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"This should not happen");
1406: break;
1407: }
1408: } else {
1409: #else
1410: {
1411: #endif
1412: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1413: VecResetArray(*Y);
1414: }
1415: PetscObjectStateIncrease((PetscObject)X);
1416: }
1417: }
1418: VecDestroy(Y);
1419: }
1420: return(0);
1421: }
1423: /*@
1424: VecGetLocalVectorRead - Maps the local portion of a vector into a
1425: vector. You must call VecRestoreLocalVectorRead() when the local
1426: vector is no longer needed.
1428: Not collective.
1430: Input parameter:
1431: . v - The vector for which the local vector is desired.
1433: Output parameter:
1434: . w - Upon exit this contains the local vector.
1436: Level: beginner
1438: Notes:
1439: This function is similar to VecGetArrayRead() which maps the local
1440: portion into a raw pointer. VecGetLocalVectorRead() is usually
1441: almost as efficient as VecGetArrayRead() but in certain circumstances
1442: VecGetLocalVectorRead() can be much more efficient than
1443: VecGetArrayRead(). This is because the construction of a contiguous
1444: array representing the vector data required by VecGetArrayRead() can
1445: be an expensive operation for certain vector types. For example, for
1446: GPU vectors VecGetArrayRead() requires that the data between device
1447: and host is synchronized.
1449: Unlike VecGetLocalVector(), this routine is not collective and
1450: preserves cached information.
1452: .seealso: VecRestoreLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1453: @*/
1454: PetscErrorCode VecGetLocalVectorRead(Vec v,Vec w)
1455: {
1457: PetscScalar *a;
1462: VecCheckSameLocalSize(v,1,w,2);
1463: if (v->ops->getlocalvectorread) {
1464: (*v->ops->getlocalvectorread)(v,w);
1465: } else {
1466: VecGetArrayRead(v,(const PetscScalar**)&a);
1467: VecPlaceArray(w,a);
1468: }
1469: return(0);
1470: }
1472: /*@
1473: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1474: previously mapped into a vector using VecGetLocalVectorRead().
1476: Not collective.
1478: Input parameter:
1479: + v - The local portion of this vector was previously mapped into w using VecGetLocalVectorRead().
1480: - w - The vector into which the local portion of v was mapped.
1482: Level: beginner
1484: .seealso: VecGetLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1485: @*/
1486: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec w)
1487: {
1489: PetscScalar *a;
1494: if (v->ops->restorelocalvectorread) {
1495: (*v->ops->restorelocalvectorread)(v,w);
1496: } else {
1497: VecGetArrayRead(w,(const PetscScalar**)&a);
1498: VecRestoreArrayRead(v,(const PetscScalar**)&a);
1499: VecResetArray(w);
1500: }
1501: return(0);
1502: }
1504: /*@
1505: VecGetLocalVector - Maps the local portion of a vector into a
1506: vector.
1508: Collective on v, not collective on w.
1510: Input parameter:
1511: . v - The vector for which the local vector is desired.
1513: Output parameter:
1514: . w - Upon exit this contains the local vector.
1516: Level: beginner
1518: Notes:
1519: This function is similar to VecGetArray() which maps the local
1520: portion into a raw pointer. VecGetLocalVector() is usually about as
1521: efficient as VecGetArray() but in certain circumstances
1522: VecGetLocalVector() can be much more efficient than VecGetArray().
1523: This is because the construction of a contiguous array representing
1524: the vector data required by VecGetArray() can be an expensive
1525: operation for certain vector types. For example, for GPU vectors
1526: VecGetArray() requires that the data between device and host is
1527: synchronized.
1529: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1530: @*/
1531: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1532: {
1534: PetscScalar *a;
1539: VecCheckSameLocalSize(v,1,w,2);
1540: if (v->ops->getlocalvector) {
1541: (*v->ops->getlocalvector)(v,w);
1542: } else {
1543: VecGetArray(v,&a);
1544: VecPlaceArray(w,a);
1545: }
1546: return(0);
1547: }
1549: /*@
1550: VecRestoreLocalVector - Unmaps the local portion of a vector
1551: previously mapped into a vector using VecGetLocalVector().
1553: Logically collective.
1555: Input parameter:
1556: + v - The local portion of this vector was previously mapped into w using VecGetLocalVector().
1557: - w - The vector into which the local portion of v was mapped.
1559: Level: beginner
1561: .seealso: VecGetLocalVector(), VecGetLocalVectorRead(), VecRestoreLocalVectorRead(), LocalVectorRead(), VecGetArrayRead(), VecGetArray()
1562: @*/
1563: PetscErrorCode VecRestoreLocalVector(Vec v,Vec w)
1564: {
1566: PetscScalar *a;
1571: if (v->ops->restorelocalvector) {
1572: (*v->ops->restorelocalvector)(v,w);
1573: } else {
1574: VecGetArray(w,&a);
1575: VecRestoreArray(v,&a);
1576: VecResetArray(w);
1577: }
1578: return(0);
1579: }
1581: /*@C
1582: VecGetArray - Returns a pointer to a contiguous array that contains this
1583: processor's portion of the vector data. For the standard PETSc
1584: vectors, VecGetArray() returns a pointer to the local data array and
1585: does not use any copies. If the underlying vector data is not stored
1586: in a contiguous array this routine will copy the data to a contiguous
1587: array and return a pointer to that. You MUST call VecRestoreArray()
1588: when you no longer need access to the array.
1590: Logically Collective on Vec
1592: Input Parameter:
1593: . x - the vector
1595: Output Parameter:
1596: . a - location to put pointer to the array
1598: Fortran Note:
1599: This routine is used differently from Fortran 77
1600: $ Vec x
1601: $ PetscScalar x_array(1)
1602: $ PetscOffset i_x
1603: $ PetscErrorCode ierr
1604: $ call VecGetArray(x,x_array,i_x,ierr)
1605: $
1606: $ Access first local entry in vector with
1607: $ value = x_array(i_x + 1)
1608: $
1609: $ ...... other code
1610: $ call VecRestoreArray(x,x_array,i_x,ierr)
1611: For Fortran 90 see VecGetArrayF90()
1613: See the Fortran chapter of the users manual and
1614: petsc/src/snes/tutorials/ex5f.F for details.
1616: Level: beginner
1618: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1619: VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1620: @*/
1621: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1622: {
1624: #if defined(PETSC_HAVE_VIENNACL)
1625: PetscBool is_viennacltype = PETSC_FALSE;
1626: #endif
1630: VecSetErrorIfLocked(x,1);
1631: if (x->petscnative) {
1632: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1633: if (x->offloadmask == PETSC_OFFLOAD_GPU) {
1634: #if defined(PETSC_HAVE_VIENNACL)
1635: PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1636: if (is_viennacltype) {
1637: VecViennaCLCopyFromGPU(x);
1638: } else
1639: #endif
1640: {
1641: #if defined(PETSC_HAVE_CUDA)
1642: VecCUDACopyFromGPU(x);
1643: #endif
1644: }
1645: } else if (x->offloadmask == PETSC_OFFLOAD_UNALLOCATED) {
1646: #if defined(PETSC_HAVE_VIENNACL)
1647: PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1648: if (is_viennacltype) {
1649: VecViennaCLAllocateCheckHost(x);
1650: } else
1651: #endif
1652: {
1653: #if defined(PETSC_HAVE_CUDA)
1654: VecCUDAAllocateCheckHost(x);
1655: #endif
1656: }
1657: }
1658: #endif
1659: *a = *((PetscScalar**)x->data);
1660: } else {
1661: if (x->ops->getarray) {
1662: (*x->ops->getarray)(x,a);
1663: } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array for vector type \"%s\"",((PetscObject)x)->type_name);
1664: }
1665: return(0);
1666: }
1668: /*@C
1669: VecGetArrayInPlace - Like VecGetArray(), but if this is a CUDA vector and it is currently offloaded to GPU,
1670: the returned pointer will be a GPU pointer to the GPU memory that contains this processor's portion of the
1671: vector data. Otherwise, it functions as VecGetArray().
1673: Logically Collective on Vec
1675: Input Parameter:
1676: . x - the vector
1678: Output Parameter:
1679: . a - location to put pointer to the array
1681: Level: beginner
1683: .seealso: VecRestoreArrayInPlace(), VecRestoreArrayInPlace(), VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(),
1684: VecPlaceArray(), VecGetArray2d(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1685: @*/
1686: PetscErrorCode VecGetArrayInPlace(Vec x,PetscScalar **a)
1687: {
1692: VecSetErrorIfLocked(x,1);
1694: #if defined(PETSC_HAVE_CUDA)
1695: if (x->petscnative && (x->offloadmask & PETSC_OFFLOAD_GPU)) { /* Prefer working on GPU when offloadmask is PETSC_OFFLOAD_BOTH */
1696: PetscBool is_cudatype = PETSC_FALSE;
1697: PetscObjectTypeCompareAny((PetscObject)x,&is_cudatype,VECSEQCUDA,VECMPICUDA,VECCUDA,"");
1698: if (is_cudatype) {
1699: VecCUDAGetArray(x,a);
1700: x->offloadmask = PETSC_OFFLOAD_GPU; /* Change the mask once GPU gets write access, don't wait until restore array */
1701: return(0);
1702: }
1703: }
1704: #endif
1705: VecGetArray(x,a);
1706: return(0);
1707: }
1709: /*@C
1710: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contains this
1711: processor's portion of the vector data. The values in this array are NOT valid, the routine calling this
1712: routine is responsible for putting values into the array; any values it does not set will be invalid
1714: Logically Collective on Vec
1716: Input Parameter:
1717: . x - the vector
1719: Output Parameter:
1720: . a - location to put pointer to the array
1722: Level: intermediate
1724: This is for vectors associate with GPUs, the vector is not copied up before giving access. If you need correct
1725: values in the array use VecGetArray()
1727: Concepts: vector^accessing local values
1729: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1730: VecGetArrayPair(), VecRestoreArrayPair(), VecGetArray(), VecRestoreArrayWrite()
1731: @*/
1732: PetscErrorCode VecGetArrayWrite(Vec x,PetscScalar **a)
1733: {
1738: VecSetErrorIfLocked(x,1);
1739: if (!x->ops->getarraywrite) {
1740: VecGetArray(x,a);
1741: } else {
1742: (*x->ops->getarraywrite)(x,a);
1743: }
1744: return(0);
1745: }
1747: /*@C
1748: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
1750: Not Collective
1752: Input Parameters:
1753: . x - the vector
1755: Output Parameter:
1756: . a - the array
1758: Level: beginner
1760: Notes:
1761: The array must be returned using a matching call to VecRestoreArrayRead().
1763: Unlike VecGetArray(), this routine is not collective and preserves cached information like vector norms.
1765: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
1766: implementations may require a copy, but must such implementations should cache the contiguous representation so that
1767: only one copy is performed when this routine is called multiple times in sequence.
1769: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1770: @*/
1771: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1772: {
1774: #if defined(PETSC_HAVE_VIENNACL)
1775: PetscBool is_viennacltype = PETSC_FALSE;
1776: #endif
1780: if (x->petscnative) {
1781: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1782: if (x->offloadmask == PETSC_OFFLOAD_GPU) {
1783: #if defined(PETSC_HAVE_VIENNACL)
1784: PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1785: if (is_viennacltype) {
1786: VecViennaCLCopyFromGPU(x);
1787: } else
1788: #endif
1789: {
1790: #if defined(PETSC_HAVE_CUDA)
1791: VecCUDACopyFromGPU(x);
1792: #endif
1793: }
1794: }
1795: #endif
1796: *a = *((PetscScalar **)x->data);
1797: } else if (x->ops->getarrayread) {
1798: (*x->ops->getarrayread)(x,a);
1799: } else {
1800: (*x->ops->getarray)(x,(PetscScalar**)a);
1801: }
1802: return(0);
1803: }
1805: /*@C
1806: VecGetArrayReadInPlace - Like VecGetArrayRead(), but if this is a CUDA vector and it is currently offloaded to GPU,
1807: the returned pointer will be a GPU pointer to the GPU memory that contains this processor's portion of the
1808: vector data. Otherwise, it functions as VecGetArrayRead().
1810: Not Collective
1812: Input Parameters:
1813: . x - the vector
1815: Output Parameter:
1816: . a - the array
1818: Level: beginner
1820: Notes:
1821: The array must be returned using a matching call to VecRestoreArrayReadInPlace().
1824: .seealso: VecRestoreArrayReadInPlace(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayInPlace()
1825: @*/
1826: PetscErrorCode VecGetArrayReadInPlace(Vec x,const PetscScalar **a)
1827: {
1832: #if defined(PETSC_HAVE_CUDA)
1833: if (x->petscnative && x->offloadmask & PETSC_OFFLOAD_GPU) {
1834: PetscBool is_cudatype = PETSC_FALSE;
1835: PetscObjectTypeCompareAny((PetscObject)x,&is_cudatype,VECSEQCUDA,VECMPICUDA,VECCUDA,"");
1836: if (is_cudatype) {
1837: VecCUDAGetArrayRead(x,a);
1838: return(0);
1839: }
1840: }
1841: #endif
1842: VecGetArrayRead(x,a);
1843: return(0);
1844: }
1846: /*@C
1847: VecGetArrays - Returns a pointer to the arrays in a set of vectors
1848: that were created by a call to VecDuplicateVecs(). You MUST call
1849: VecRestoreArrays() when you no longer need access to the array.
1851: Logically Collective on Vec
1853: Input Parameter:
1854: + x - the vectors
1855: - n - the number of vectors
1857: Output Parameter:
1858: . a - location to put pointer to the array
1860: Fortran Note:
1861: This routine is not supported in Fortran.
1863: Level: intermediate
1865: .seealso: VecGetArray(), VecRestoreArrays()
1866: @*/
1867: PetscErrorCode VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1868: {
1870: PetscInt i;
1871: PetscScalar **q;
1877: if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1878: PetscMalloc1(n,&q);
1879: for (i=0; i<n; ++i) {
1880: VecGetArray(x[i],&q[i]);
1881: }
1882: *a = q;
1883: return(0);
1884: }
1886: /*@C
1887: VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1888: has been called.
1890: Logically Collective on Vec
1892: Input Parameters:
1893: + x - the vector
1894: . n - the number of vectors
1895: - a - location of pointer to arrays obtained from VecGetArrays()
1897: Notes:
1898: For regular PETSc vectors this routine does not involve any copies. For
1899: any special vectors that do not store local vector data in a contiguous
1900: array, this routine will copy the data back into the underlying
1901: vector data structure from the arrays obtained with VecGetArrays().
1903: Fortran Note:
1904: This routine is not supported in Fortran.
1906: Level: intermediate
1908: .seealso: VecGetArrays(), VecRestoreArray()
1909: @*/
1910: PetscErrorCode VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1911: {
1913: PetscInt i;
1914: PetscScalar **q = *a;
1921: for (i=0; i<n; ++i) {
1922: VecRestoreArray(x[i],&q[i]);
1923: }
1924: PetscFree(q);
1925: return(0);
1926: }
1928: /*@C
1929: VecRestoreArray - Restores a vector after VecGetArray() has been called.
1931: Logically Collective on Vec
1933: Input Parameters:
1934: + x - the vector
1935: - a - location of pointer to array obtained from VecGetArray()
1937: Level: beginner
1939: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1940: VecGetArrayPair(), VecRestoreArrayPair()
1941: @*/
1942: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1943: {
1948: if (x->petscnative) {
1949: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1950: x->offloadmask = PETSC_OFFLOAD_CPU;
1951: #endif
1952: } else {
1953: (*x->ops->restorearray)(x,a);
1954: }
1955: if (a) *a = NULL;
1956: PetscObjectStateIncrease((PetscObject)x);
1957: return(0);
1958: }
1960: /*@C
1961: VecRestoreArrayInPlace - Restores a vector after VecGetArrayInPlace() has been called.
1963: Logically Collective on Vec
1965: Input Parameters:
1966: + x - the vector
1967: - a - location of pointer to array obtained from VecGetArrayInPlace()
1969: Level: beginner
1971: .seealso: VecGetArrayInPlace(), VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(),
1972: VecPlaceArray(), VecRestoreArray2d(), VecGetArrayPair(), VecRestoreArrayPair()
1973: @*/
1974: PetscErrorCode VecRestoreArrayInPlace(Vec x,PetscScalar **a)
1975: {
1980: #if defined(PETSC_HAVE_CUDA)
1981: if (x->petscnative && x->offloadmask == PETSC_OFFLOAD_GPU) {
1982: PetscBool is_cudatype = PETSC_FALSE;
1983: PetscObjectTypeCompareAny((PetscObject)x,&is_cudatype,VECSEQCUDA,VECMPICUDA,VECCUDA,"");
1984: if (is_cudatype) {
1985: VecCUDARestoreArray(x,a);
1986: return(0);
1987: }
1988: }
1989: #endif
1990: VecRestoreArray(x,a);
1991: return(0);
1992: }
1995: /*@C
1996: VecRestoreArrayWrite - Restores a vector after VecGetArrayWrite() has been called.
1998: Logically Collective on Vec
2000: Input Parameters:
2001: + x - the vector
2002: - a - location of pointer to array obtained from VecGetArray()
2004: Level: beginner
2006: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
2007: VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite()
2008: @*/
2009: PetscErrorCode VecRestoreArrayWrite(Vec x,PetscScalar **a)
2010: {
2015: if (x->petscnative) {
2016: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
2017: x->offloadmask = PETSC_OFFLOAD_CPU;
2018: #endif
2019: } else {
2020: if (x->ops->restorearraywrite) {
2021: (*x->ops->restorearraywrite)(x,a);
2022: } else {
2023: (*x->ops->restorearray)(x,a);
2024: }
2025: }
2026: if (a) *a = NULL;
2027: PetscObjectStateIncrease((PetscObject)x);
2028: return(0);
2029: }
2031: /*@C
2032: VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()
2034: Not Collective
2036: Input Parameters:
2037: + vec - the vector
2038: - array - the array
2040: Level: beginner
2042: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
2043: @*/
2044: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
2045: {
2050: if (x->petscnative) {
2051: /* nothing */
2052: } else if (x->ops->restorearrayread) {
2053: (*x->ops->restorearrayread)(x,a);
2054: } else {
2055: (*x->ops->restorearray)(x,(PetscScalar**)a);
2056: }
2057: if (a) *a = NULL;
2058: return(0);
2059: }
2061: /*@C
2062: VecRestoreArrayReadInPlace - Restore array obtained with VecGetArrayReadInPlace()
2064: Not Collective
2066: Input Parameters:
2067: + vec - the vector
2068: - array - the array
2070: Level: beginner
2072: .seealso: VecGetArrayReadInPlace(), VecRestoreArrayInPlace(), VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
2073: @*/
2074: PetscErrorCode VecRestoreArrayReadInPlace(Vec x,const PetscScalar **a)
2075: {
2079: VecRestoreArrayRead(x,a);
2080: return(0);
2081: }
2083: /*@
2084: VecPlaceArray - Allows one to replace the array in a vector with an
2085: array provided by the user. This is useful to avoid copying an array
2086: into a vector.
2088: Not Collective
2090: Input Parameters:
2091: + vec - the vector
2092: - array - the array
2094: Notes:
2095: You can return to the original array with a call to VecResetArray()
2097: Level: developer
2099: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecResetArray()
2101: @*/
2102: PetscErrorCode VecPlaceArray(Vec vec,const PetscScalar array[])
2103: {
2110: if (vec->ops->placearray) {
2111: (*vec->ops->placearray)(vec,array);
2112: } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
2113: PetscObjectStateIncrease((PetscObject)vec);
2114: return(0);
2115: }
2117: /*@C
2118: VecReplaceArray - Allows one to replace the array in a vector with an
2119: array provided by the user. This is useful to avoid copying an array
2120: into a vector.
2122: Not Collective
2124: Input Parameters:
2125: + vec - the vector
2126: - array - the array
2128: Notes:
2129: This permanently replaces the array and frees the memory associated
2130: with the old array.
2132: The memory passed in MUST be obtained with PetscMalloc() and CANNOT be
2133: freed by the user. It will be freed when the vector is destroyed.
2135: Not supported from Fortran
2137: Level: developer
2139: .seealso: VecGetArray(), VecRestoreArray(), VecPlaceArray(), VecResetArray()
2141: @*/
2142: PetscErrorCode VecReplaceArray(Vec vec,const PetscScalar array[])
2143: {
2149: if (vec->ops->replacearray) {
2150: (*vec->ops->replacearray)(vec,array);
2151: } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
2152: PetscObjectStateIncrease((PetscObject)vec);
2153: return(0);
2154: }
2157: /*@C
2158: VecCUDAGetArray - Provides access to the CUDA buffer inside a vector.
2160: This function has semantics similar to VecGetArray(): the pointer
2161: returned by this function points to a consistent view of the vector
2162: data. This may involve a copy operation of data from the host to the
2163: device if the data on the device is out of date. If the device
2164: memory hasn't been allocated previously it will be allocated as part
2165: of this function call. VecCUDAGetArray() assumes that
2166: the user will modify the vector data. This is similar to
2167: intent(inout) in fortran.
2169: The CUDA device pointer has to be released by calling
2170: VecCUDARestoreArray(). Upon restoring the vector data
2171: the data on the host will be marked as out of date. A subsequent
2172: access of the host data will thus incur a data transfer from the
2173: device to the host.
2176: Input Parameter:
2177: . v - the vector
2179: Output Parameter:
2180: . a - the CUDA device pointer
2182: Fortran note:
2183: This function is not currently available from Fortran.
2185: Level: intermediate
2187: .seealso: VecCUDARestoreArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2188: @*/
2189: PETSC_EXTERN PetscErrorCode VecCUDAGetArray(Vec v, PetscScalar **a)
2190: {
2191: #if defined(PETSC_HAVE_CUDA)
2193: #endif
2197: #if defined(PETSC_HAVE_CUDA)
2198: *a = 0;
2199: VecCUDACopyToGPU(v);
2200: *a = ((Vec_CUDA*)v->spptr)->GPUarray;
2201: #endif
2202: return(0);
2203: }
2205: /*@C
2206: VecCUDARestoreArray - Restore a CUDA device pointer previously acquired with VecCUDAGetArray().
2208: This marks the host data as out of date. Subsequent access to the
2209: vector data on the host side with for instance VecGetArray() incurs a
2210: data transfer.
2212: Input Parameter:
2213: + v - the vector
2214: - a - the CUDA device pointer. This pointer is invalid after
2215: VecCUDARestoreArray() returns.
2217: Fortran note:
2218: This function is not currently available from Fortran.
2220: Level: intermediate
2222: .seealso: VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2223: @*/
2224: PETSC_EXTERN PetscErrorCode VecCUDARestoreArray(Vec v, PetscScalar **a)
2225: {
2230: #if defined(PETSC_HAVE_CUDA)
2231: v->offloadmask = PETSC_OFFLOAD_GPU;
2232: #endif
2234: PetscObjectStateIncrease((PetscObject)v);
2235: return(0);
2236: }
2238: /*@C
2239: VecCUDAGetArrayRead - Provides read access to the CUDA buffer inside a vector.
2241: This function is analogous to VecGetArrayRead(): The pointer
2242: returned by this function points to a consistent view of the vector
2243: data. This may involve a copy operation of data from the host to the
2244: device if the data on the device is out of date. If the device
2245: memory hasn't been allocated previously it will be allocated as part
2246: of this function call. VecCUDAGetArrayRead() assumes that the
2247: user will not modify the vector data. This is analgogous to
2248: intent(in) in Fortran.
2250: The CUDA device pointer has to be released by calling
2251: VecCUDARestoreArrayRead(). If the data on the host side was
2252: previously up to date it will remain so, i.e. data on both the device
2253: and the host is up to date. Accessing data on the host side does not
2254: incur a device to host data transfer.
2256: Input Parameter:
2257: . v - the vector
2259: Output Parameter:
2260: . a - the CUDA pointer.
2262: Fortran note:
2263: This function is not currently available from Fortran.
2265: Level: intermediate
2267: .seealso: VecCUDARestoreArrayRead(), VecCUDAGetArray(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2268: @*/
2269: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayRead(Vec v, const PetscScalar **a)
2270: {
2271: #if defined(PETSC_HAVE_CUDA)
2273: #endif
2277: #if defined(PETSC_HAVE_CUDA)
2278: *a = 0;
2279: VecCUDACopyToGPU(v);
2280: *a = ((Vec_CUDA*)v->spptr)->GPUarray;
2281: #endif
2282: return(0);
2283: }
2285: /*@C
2286: VecCUDARestoreArrayRead - Restore a CUDA device pointer previously acquired with VecCUDAGetArrayRead().
2288: If the data on the host side was previously up to date it will remain
2289: so, i.e. data on both the device and the host is up to date.
2290: Accessing data on the host side e.g. with VecGetArray() does not
2291: incur a device to host data transfer.
2293: Input Parameter:
2294: + v - the vector
2295: - a - the CUDA device pointer. This pointer is invalid after
2296: VecCUDARestoreArrayRead() returns.
2298: Fortran note:
2299: This function is not currently available from Fortran.
2301: Level: intermediate
2303: .seealso: VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2304: @*/
2305: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayRead(Vec v, const PetscScalar **a)
2306: {
2309: *a = NULL;
2310: return(0);
2311: }
2313: /*@C
2314: VecCUDAGetArrayWrite - Provides write access to the CUDA buffer inside a vector.
2316: The data pointed to by the device pointer is uninitialized. The user
2317: may not read from this data. Furthermore, the entire array needs to
2318: be filled by the user to obtain well-defined behaviour. The device
2319: memory will be allocated by this function if it hasn't been allocated
2320: previously. This is analogous to intent(out) in Fortran.
2322: The device pointer needs to be released with
2323: VecCUDARestoreArrayWrite(). When the pointer is released the
2324: host data of the vector is marked as out of data. Subsequent access
2325: of the host data with e.g. VecGetArray() incurs a device to host data
2326: transfer.
2329: Input Parameter:
2330: . v - the vector
2332: Output Parameter:
2333: . a - the CUDA pointer
2335: Fortran note:
2336: This function is not currently available from Fortran.
2338: Level: advanced
2340: .seealso: VecCUDARestoreArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecGetArrayRead()
2341: @*/
2342: PETSC_EXTERN PetscErrorCode VecCUDAGetArrayWrite(Vec v, PetscScalar **a)
2343: {
2344: #if defined(PETSC_HAVE_CUDA)
2346: #endif
2350: #if defined(PETSC_HAVE_CUDA)
2351: *a = 0;
2352: VecCUDAAllocateCheck(v);
2353: *a = ((Vec_CUDA*)v->spptr)->GPUarray;
2354: #endif
2355: return(0);
2356: }
2358: /*@C
2359: VecCUDARestoreArrayWrite - Restore a CUDA device pointer previously acquired with VecCUDAGetArrayWrite().
2361: Data on the host will be marked as out of date. Subsequent access of
2362: the data on the host side e.g. with VecGetArray() will incur a device
2363: to host data transfer.
2365: Input Parameter:
2366: + v - the vector
2367: - a - the CUDA device pointer. This pointer is invalid after
2368: VecCUDARestoreArrayWrite() returns.
2370: Fortran note:
2371: This function is not currently available from Fortran.
2373: Level: intermediate
2375: .seealso: VecCUDAGetArrayWrite(), VecCUDAGetArray(), VecCUDAGetArrayRead(), VecCUDAGetArrayWrite(), VecGetArray(), VecRestoreArray(), VecGetArrayRead()
2376: @*/
2377: PETSC_EXTERN PetscErrorCode VecCUDARestoreArrayWrite(Vec v, PetscScalar **a)
2378: {
2383: #if defined(PETSC_HAVE_CUDA)
2384: v->offloadmask = PETSC_OFFLOAD_GPU;
2385: #endif
2387: PetscObjectStateIncrease((PetscObject)v);
2388: return(0);
2389: }
2391: /*@C
2392: VecCUDAPlaceArray - Allows one to replace the GPU array in a vector with a
2393: GPU array provided by the user. This is useful to avoid copying an
2394: array into a vector.
2396: Not Collective
2398: Input Parameters:
2399: + vec - the vector
2400: - array - the GPU array
2402: Notes:
2403: You can return to the original GPU array with a call to VecCUDAResetArray()
2404: It is not possible to use VecCUDAPlaceArray() and VecPlaceArray() at the
2405: same time on the same vector.
2407: Level: developer
2409: .seealso: VecPlaceArray(), VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecResetArray(), VecCUDAResetArray(), VecCUDAReplaceArray()
2411: @*/
2412: PetscErrorCode VecCUDAPlaceArray(Vec vin,const PetscScalar a[])
2413: {
2418: #if defined(PETSC_HAVE_CUDA)
2419: VecCUDACopyToGPU(vin);
2420: 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()");
2421: ((Vec_Seq*)vin->data)->unplacedarray = (PetscScalar *) ((Vec_CUDA*)vin->spptr)->GPUarray; /* save previous GPU array so reset can bring it back */
2422: ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2423: vin->offloadmask = PETSC_OFFLOAD_GPU;
2424: #endif
2425: PetscObjectStateIncrease((PetscObject)vin);
2426: return(0);
2427: }
2429: /*@C
2430: VecCUDAReplaceArray - Allows one to replace the GPU array in a vector
2431: with a GPU array provided by the user. This is useful to avoid copying
2432: a GPU array into a vector.
2434: Not Collective
2436: Input Parameters:
2437: + vec - the vector
2438: - array - the GPU array
2440: Notes:
2441: This permanently replaces the GPU array and frees the memory associated
2442: with the old GPU array.
2444: The memory passed in CANNOT be freed by the user. It will be freed
2445: when the vector is destroyed.
2447: Not supported from Fortran
2449: Level: developer
2451: .seealso: VecGetArray(), VecRestoreArray(), VecPlaceArray(), VecResetArray(), VecCUDAResetArray(), VecCUDAPlaceArray(), VecReplaceArray()
2453: @*/
2454: PetscErrorCode VecCUDAReplaceArray(Vec vin,const PetscScalar a[])
2455: {
2456: #if defined(PETSC_HAVE_CUDA)
2457: cudaError_t err;
2458: #endif
2463: #if defined(PETSC_HAVE_CUDA)
2464: err = cudaFree(((Vec_CUDA*)vin->spptr)->GPUarray);CHKERRCUDA(err);
2465: ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar*)a;
2466: vin->offloadmask = PETSC_OFFLOAD_GPU;
2467: #endif
2468: PetscObjectStateIncrease((PetscObject)vin);
2469: return(0);
2470: }
2472: /*@C
2473: VecCUDAResetArray - Resets a vector to use its default memory. Call this
2474: after the use of VecCUDAPlaceArray().
2476: Not Collective
2478: Input Parameters:
2479: . vec - the vector
2481: Level: developer
2483: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecPlaceArray(), VecResetArray(), VecCUDAPlaceArray(), VecCUDAReplaceArray()
2485: @*/
2486: PetscErrorCode VecCUDAResetArray(Vec vin)
2487: {
2492: #if defined(PETSC_HAVE_CUDA)
2493: VecCUDACopyToGPU(vin);
2494: ((Vec_CUDA*)vin->spptr)->GPUarray = (PetscScalar *) ((Vec_Seq*)vin->data)->unplacedarray;
2495: ((Vec_Seq*)vin->data)->unplacedarray = 0;
2496: vin->offloadmask = PETSC_OFFLOAD_GPU;
2497: #endif
2498: PetscObjectStateIncrease((PetscObject)vin);
2499: return(0);
2500: }
2502: /*MC
2503: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
2504: and makes them accessible via a Fortran90 pointer.
2506: Synopsis:
2507: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
2509: Collective on Vec
2511: Input Parameters:
2512: + x - a vector to mimic
2513: - n - the number of vectors to obtain
2515: Output Parameters:
2516: + y - Fortran90 pointer to the array of vectors
2517: - ierr - error code
2519: Example of Usage:
2520: .vb
2521: #include <petsc/finclude/petscvec.h>
2522: use petscvec
2524: Vec x
2525: Vec, pointer :: y(:)
2526: ....
2527: call VecDuplicateVecsF90(x,2,y,ierr)
2528: call VecSet(y(2),alpha,ierr)
2529: call VecSet(y(2),alpha,ierr)
2530: ....
2531: call VecDestroyVecsF90(2,y,ierr)
2532: .ve
2534: Notes:
2535: Not yet supported for all F90 compilers
2537: Use VecDestroyVecsF90() to free the space.
2539: Level: beginner
2541: .seealso: VecDestroyVecsF90(), VecDuplicateVecs()
2543: M*/
2545: /*MC
2546: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2547: VecGetArrayF90().
2549: Synopsis:
2550: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2552: Logically Collective on Vec
2554: Input Parameters:
2555: + x - vector
2556: - xx_v - the Fortran90 pointer to the array
2558: Output Parameter:
2559: . ierr - error code
2561: Example of Usage:
2562: .vb
2563: #include <petsc/finclude/petscvec.h>
2564: use petscvec
2566: PetscScalar, pointer :: xx_v(:)
2567: ....
2568: call VecGetArrayF90(x,xx_v,ierr)
2569: xx_v(3) = a
2570: call VecRestoreArrayF90(x,xx_v,ierr)
2571: .ve
2573: Level: beginner
2575: .seealso: VecGetArrayF90(), VecGetArray(), VecRestoreArray(), UsingFortran, VecRestoreArrayReadF90()
2577: M*/
2579: /*MC
2580: VecDestroyVecsF90 - Frees a block of vectors obtained with VecDuplicateVecsF90().
2582: Synopsis:
2583: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
2585: Collective on Vec
2587: Input Parameters:
2588: + n - the number of vectors previously obtained
2589: - x - pointer to array of vector pointers
2591: Output Parameter:
2592: . ierr - error code
2594: Notes:
2595: Not yet supported for all F90 compilers
2597: Level: beginner
2599: .seealso: VecDestroyVecs(), VecDuplicateVecsF90()
2601: M*/
2603: /*MC
2604: VecGetArrayF90 - Accesses a vector array from Fortran90. For default PETSc
2605: vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
2606: this routine is implementation dependent. You MUST call VecRestoreArrayF90()
2607: when you no longer need access to the array.
2609: Synopsis:
2610: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2612: Logically Collective on Vec
2614: Input Parameter:
2615: . x - vector
2617: Output Parameters:
2618: + xx_v - the Fortran90 pointer to the array
2619: - ierr - error code
2621: Example of Usage:
2622: .vb
2623: #include <petsc/finclude/petscvec.h>
2624: use petscvec
2626: PetscScalar, pointer :: xx_v(:)
2627: ....
2628: call VecGetArrayF90(x,xx_v,ierr)
2629: xx_v(3) = a
2630: call VecRestoreArrayF90(x,xx_v,ierr)
2631: .ve
2633: If you ONLY intend to read entries from the array and not change any entries you should use VecGetArrayReadF90().
2635: Level: beginner
2637: .seealso: VecRestoreArrayF90(), VecGetArray(), VecRestoreArray(), VecGetArrayReadF90(), UsingFortran
2639: M*/
2641: /*MC
2642: VecGetArrayReadF90 - Accesses a read only array from Fortran90. For default PETSc
2643: vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
2644: this routine is implementation dependent. You MUST call VecRestoreArrayReadF90()
2645: when you no longer need access to the array.
2647: Synopsis:
2648: VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2650: Logically Collective on Vec
2652: Input Parameter:
2653: . x - vector
2655: Output Parameters:
2656: + xx_v - the Fortran90 pointer to the array
2657: - ierr - error code
2659: Example of Usage:
2660: .vb
2661: #include <petsc/finclude/petscvec.h>
2662: use petscvec
2664: PetscScalar, pointer :: xx_v(:)
2665: ....
2666: call VecGetArrayReadF90(x,xx_v,ierr)
2667: a = xx_v(3)
2668: call VecRestoreArrayReadF90(x,xx_v,ierr)
2669: .ve
2671: If you intend to write entries into the array you must use VecGetArrayF90().
2673: Level: beginner
2675: .seealso: VecRestoreArrayReadF90(), VecGetArray(), VecRestoreArray(), VecGetArrayRead(), VecRestoreArrayRead(), VecGetArrayF90(), UsingFortran
2677: M*/
2679: /*MC
2680: VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2681: VecGetArrayReadF90().
2683: Synopsis:
2684: VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2686: Logically Collective on Vec
2688: Input Parameters:
2689: + x - vector
2690: - xx_v - the Fortran90 pointer to the array
2692: Output Parameter:
2693: . ierr - error code
2695: Example of Usage:
2696: .vb
2697: #include <petsc/finclude/petscvec.h>
2698: use petscvec
2700: PetscScalar, pointer :: xx_v(:)
2701: ....
2702: call VecGetArrayReadF90(x,xx_v,ierr)
2703: a = xx_v(3)
2704: call VecRestoreArrayReadF90(x,xx_v,ierr)
2705: .ve
2707: Level: beginner
2709: .seealso: VecGetArrayReadF90(), VecGetArray(), VecRestoreArray(), VecGetArrayRead(), VecRestoreArrayRead(),UsingFortran, VecRestoreArrayF90()
2711: M*/
2713: /*@C
2714: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2715: processor's portion of the vector data. You MUST call VecRestoreArray2d()
2716: when you no longer need access to the array.
2718: Logically Collective
2720: Input Parameter:
2721: + x - the vector
2722: . m - first dimension of two dimensional array
2723: . n - second dimension of two dimensional array
2724: . mstart - first index you will use in first coordinate direction (often 0)
2725: - nstart - first index in the second coordinate direction (often 0)
2727: Output Parameter:
2728: . a - location to put pointer to the array
2730: Level: developer
2732: Notes:
2733: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2734: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2735: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2736: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
2738: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2740: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2741: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2742: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2743: @*/
2744: PetscErrorCode VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2745: {
2747: PetscInt i,N;
2748: PetscScalar *aa;
2754: VecGetLocalSize(x,&N);
2755: 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);
2756: VecGetArray(x,&aa);
2758: PetscMalloc1(m,a);
2759: for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2760: *a -= mstart;
2761: return(0);
2762: }
2764: /*@C
2765: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2766: processor's portion of the vector data. You MUST call VecRestoreArray2dWrite()
2767: when you no longer need access to the array.
2769: Logically Collective
2771: Input Parameter:
2772: + x - the vector
2773: . m - first dimension of two dimensional array
2774: . n - second dimension of two dimensional array
2775: . mstart - first index you will use in first coordinate direction (often 0)
2776: - nstart - first index in the second coordinate direction (often 0)
2778: Output Parameter:
2779: . a - location to put pointer to the array
2781: Level: developer
2783: Notes:
2784: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2785: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2786: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2787: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
2789: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2791: Concepts: vector^accessing local values as 2d array
2793: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2794: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2795: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2796: @*/
2797: PetscErrorCode VecGetArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2798: {
2800: PetscInt i,N;
2801: PetscScalar *aa;
2807: VecGetLocalSize(x,&N);
2808: 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);
2809: VecGetArrayWrite(x,&aa);
2811: PetscMalloc1(m,a);
2812: for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2813: *a -= mstart;
2814: return(0);
2815: }
2817: /*@C
2818: VecRestoreArray2d - Restores a vector after VecGetArray2d() has been called.
2820: Logically Collective
2822: Input Parameters:
2823: + x - the vector
2824: . m - first dimension of two dimensional array
2825: . n - second dimension of the two dimensional array
2826: . mstart - first index you will use in first coordinate direction (often 0)
2827: . nstart - first index in the second coordinate direction (often 0)
2828: - a - location of pointer to array obtained from VecGetArray2d()
2830: Level: developer
2832: Notes:
2833: For regular PETSc vectors this routine does not involve any copies. For
2834: any special vectors that do not store local vector data in a contiguous
2835: array, this routine will copy the data back into the underlying
2836: vector data structure from the array obtained with VecGetArray().
2838: This routine actually zeros out the a pointer.
2840: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2841: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2842: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2843: @*/
2844: PetscErrorCode VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2845: {
2847: void *dummy;
2853: dummy = (void*)(*a + mstart);
2854: PetscFree(dummy);
2855: VecRestoreArray(x,NULL);
2856: return(0);
2857: }
2859: /*@C
2860: VecRestoreArray2dWrite - Restores a vector after VecGetArray2dWrite() has been called.
2862: Logically Collective
2864: Input Parameters:
2865: + x - the vector
2866: . m - first dimension of two dimensional array
2867: . n - second dimension of the two dimensional array
2868: . mstart - first index you will use in first coordinate direction (often 0)
2869: . nstart - first index in the second coordinate direction (often 0)
2870: - a - location of pointer to array obtained from VecGetArray2d()
2872: Level: developer
2874: Notes:
2875: For regular PETSc vectors this routine does not involve any copies. For
2876: any special vectors that do not store local vector data in a contiguous
2877: array, this routine will copy the data back into the underlying
2878: vector data structure from the array obtained with VecGetArray().
2880: This routine actually zeros out the a pointer.
2882: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2883: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2884: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2885: @*/
2886: PetscErrorCode VecRestoreArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2887: {
2889: void *dummy;
2895: dummy = (void*)(*a + mstart);
2896: PetscFree(dummy);
2897: VecRestoreArrayWrite(x,NULL);
2898: return(0);
2899: }
2901: /*@C
2902: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2903: processor's portion of the vector data. You MUST call VecRestoreArray1d()
2904: when you no longer need access to the array.
2906: Logically Collective
2908: Input Parameter:
2909: + x - the vector
2910: . m - first dimension of two dimensional array
2911: - mstart - first index you will use in first coordinate direction (often 0)
2913: Output Parameter:
2914: . a - location to put pointer to the array
2916: Level: developer
2918: Notes:
2919: For a vector obtained from DMCreateLocalVector() mstart are likely
2920: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2921: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
2923: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2925: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2926: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2927: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2928: @*/
2929: PetscErrorCode VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2930: {
2932: PetscInt N;
2938: VecGetLocalSize(x,&N);
2939: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2940: VecGetArray(x,a);
2941: *a -= mstart;
2942: return(0);
2943: }
2945: /*@C
2946: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
2947: processor's portion of the vector data. You MUST call VecRestoreArray1dWrite()
2948: when you no longer need access to the array.
2950: Logically Collective
2952: Input Parameter:
2953: + x - the vector
2954: . m - first dimension of two dimensional array
2955: - mstart - first index you will use in first coordinate direction (often 0)
2957: Output Parameter:
2958: . a - location to put pointer to the array
2960: Level: developer
2962: Notes:
2963: For a vector obtained from DMCreateLocalVector() mstart are likely
2964: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2965: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
2967: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2969: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2970: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2971: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2972: @*/
2973: PetscErrorCode VecGetArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2974: {
2976: PetscInt N;
2982: VecGetLocalSize(x,&N);
2983: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2984: VecGetArrayWrite(x,a);
2985: *a -= mstart;
2986: return(0);
2987: }
2989: /*@C
2990: VecRestoreArray1d - Restores a vector after VecGetArray1d() has been called.
2992: Logically Collective
2994: Input Parameters:
2995: + x - the vector
2996: . m - first dimension of two dimensional array
2997: . mstart - first index you will use in first coordinate direction (often 0)
2998: - a - location of pointer to array obtained from VecGetArray21()
3000: Level: developer
3002: Notes:
3003: For regular PETSc vectors this routine does not involve any copies. For
3004: any special vectors that do not store local vector data in a contiguous
3005: array, this routine will copy the data back into the underlying
3006: vector data structure from the array obtained with VecGetArray1d().
3008: This routine actually zeros out the a pointer.
3010: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3011: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3012: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3013: @*/
3014: PetscErrorCode VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3015: {
3021: VecRestoreArray(x,NULL);
3022: return(0);
3023: }
3025: /*@C
3026: VecRestoreArray1dWrite - Restores a vector after VecGetArray1dWrite() has been called.
3028: Logically Collective
3030: Input Parameters:
3031: + x - the vector
3032: . m - first dimension of two dimensional array
3033: . mstart - first index you will use in first coordinate direction (often 0)
3034: - a - location of pointer to array obtained from VecGetArray21()
3036: Level: developer
3038: Notes:
3039: For regular PETSc vectors this routine does not involve any copies. For
3040: any special vectors that do not store local vector data in a contiguous
3041: array, this routine will copy the data back into the underlying
3042: vector data structure from the array obtained with VecGetArray1d().
3044: This routine actually zeros out the a pointer.
3046: Concepts: vector^accessing local values as 1d array
3048: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3049: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3050: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3051: @*/
3052: PetscErrorCode VecRestoreArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3053: {
3059: VecRestoreArrayWrite(x,NULL);
3060: return(0);
3061: }
3063: /*@C
3064: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
3065: processor's portion of the vector data. You MUST call VecRestoreArray3d()
3066: when you no longer need access to the array.
3068: Logically Collective
3070: Input Parameter:
3071: + x - the vector
3072: . m - first dimension of three dimensional array
3073: . n - second dimension of three dimensional array
3074: . p - third dimension of three dimensional array
3075: . mstart - first index you will use in first coordinate direction (often 0)
3076: . nstart - first index in the second coordinate direction (often 0)
3077: - pstart - first index in the third coordinate direction (often 0)
3079: Output Parameter:
3080: . a - location to put pointer to the array
3082: Level: developer
3084: Notes:
3085: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3086: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3087: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3088: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
3090: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3092: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3093: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3094: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3095: @*/
3096: PetscErrorCode VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3097: {
3099: PetscInt i,N,j;
3100: PetscScalar *aa,**b;
3106: VecGetLocalSize(x,&N);
3107: 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);
3108: VecGetArray(x,&aa);
3110: PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3111: b = (PetscScalar**)((*a) + m);
3112: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3113: for (i=0; i<m; i++)
3114: for (j=0; j<n; j++)
3115: b[i*n+j] = aa + i*n*p + j*p - pstart;
3117: *a -= mstart;
3118: return(0);
3119: }
3121: /*@C
3122: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
3123: processor's portion of the vector data. You MUST call VecRestoreArray3dWrite()
3124: when you no longer need access to the array.
3126: Logically Collective
3128: Input Parameter:
3129: + x - the vector
3130: . m - first dimension of three dimensional array
3131: . n - second dimension of three dimensional array
3132: . p - third dimension of three dimensional array
3133: . mstart - first index you will use in first coordinate direction (often 0)
3134: . nstart - first index in the second coordinate direction (often 0)
3135: - pstart - first index in the third coordinate direction (often 0)
3137: Output Parameter:
3138: . a - location to put pointer to the array
3140: Level: developer
3142: Notes:
3143: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3144: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3145: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3146: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
3148: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3150: Concepts: vector^accessing local values as 3d array
3152: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3153: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3154: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3155: @*/
3156: PetscErrorCode VecGetArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3157: {
3159: PetscInt i,N,j;
3160: PetscScalar *aa,**b;
3166: VecGetLocalSize(x,&N);
3167: 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);
3168: VecGetArrayWrite(x,&aa);
3170: PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3171: b = (PetscScalar**)((*a) + m);
3172: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3173: for (i=0; i<m; i++)
3174: for (j=0; j<n; j++)
3175: b[i*n+j] = aa + i*n*p + j*p - pstart;
3177: *a -= mstart;
3178: return(0);
3179: }
3181: /*@C
3182: VecRestoreArray3d - Restores a vector after VecGetArray3d() has been called.
3184: Logically Collective
3186: Input Parameters:
3187: + x - the vector
3188: . m - first dimension of three dimensional array
3189: . n - second dimension of the three dimensional array
3190: . p - third dimension of the three dimensional array
3191: . mstart - first index you will use in first coordinate direction (often 0)
3192: . nstart - first index in the second coordinate direction (often 0)
3193: . pstart - first index in the third coordinate direction (often 0)
3194: - a - location of pointer to array obtained from VecGetArray3d()
3196: Level: developer
3198: Notes:
3199: For regular PETSc vectors this routine does not involve any copies. For
3200: any special vectors that do not store local vector data in a contiguous
3201: array, this routine will copy the data back into the underlying
3202: vector data structure from the array obtained with VecGetArray().
3204: This routine actually zeros out the a pointer.
3206: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3207: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3208: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3209: @*/
3210: PetscErrorCode VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3211: {
3213: void *dummy;
3219: dummy = (void*)(*a + mstart);
3220: PetscFree(dummy);
3221: VecRestoreArray(x,NULL);
3222: return(0);
3223: }
3225: /*@C
3226: VecRestoreArray3dWrite - Restores a vector after VecGetArray3dWrite() has been called.
3228: Logically Collective
3230: Input Parameters:
3231: + x - the vector
3232: . m - first dimension of three dimensional array
3233: . n - second dimension of the three dimensional array
3234: . p - third dimension of the three dimensional array
3235: . mstart - first index you will use in first coordinate direction (often 0)
3236: . nstart - first index in the second coordinate direction (often 0)
3237: . pstart - first index in the third coordinate direction (often 0)
3238: - a - location of pointer to array obtained from VecGetArray3d()
3240: Level: developer
3242: Notes:
3243: For regular PETSc vectors this routine does not involve any copies. For
3244: any special vectors that do not store local vector data in a contiguous
3245: array, this routine will copy the data back into the underlying
3246: vector data structure from the array obtained with VecGetArray().
3248: This routine actually zeros out the a pointer.
3250: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3251: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3252: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3253: @*/
3254: PetscErrorCode VecRestoreArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3255: {
3257: void *dummy;
3263: dummy = (void*)(*a + mstart);
3264: PetscFree(dummy);
3265: VecRestoreArrayWrite(x,NULL);
3266: return(0);
3267: }
3269: /*@C
3270: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3271: processor's portion of the vector data. You MUST call VecRestoreArray4d()
3272: when you no longer need access to the array.
3274: Logically Collective
3276: Input Parameter:
3277: + x - the vector
3278: . m - first dimension of four dimensional array
3279: . n - second dimension of four dimensional array
3280: . p - third dimension of four dimensional array
3281: . q - fourth dimension of four dimensional array
3282: . mstart - first index you will use in first coordinate direction (often 0)
3283: . nstart - first index in the second coordinate direction (often 0)
3284: . pstart - first index in the third coordinate direction (often 0)
3285: - qstart - first index in the fourth coordinate direction (often 0)
3287: Output Parameter:
3288: . a - location to put pointer to the array
3290: Level: beginner
3292: Notes:
3293: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3294: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3295: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3296: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
3298: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3300: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3301: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3302: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3303: @*/
3304: PetscErrorCode VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3305: {
3307: PetscInt i,N,j,k;
3308: PetscScalar *aa,***b,**c;
3314: VecGetLocalSize(x,&N);
3315: 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);
3316: VecGetArray(x,&aa);
3318: PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3319: b = (PetscScalar***)((*a) + m);
3320: c = (PetscScalar**)(b + m*n);
3321: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3322: for (i=0; i<m; i++)
3323: for (j=0; j<n; j++)
3324: b[i*n+j] = c + i*n*p + j*p - pstart;
3325: for (i=0; i<m; i++)
3326: for (j=0; j<n; j++)
3327: for (k=0; k<p; k++)
3328: c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3329: *a -= mstart;
3330: return(0);
3331: }
3333: /*@C
3334: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3335: processor's portion of the vector data. You MUST call VecRestoreArray4dWrite()
3336: when you no longer need access to the array.
3338: Logically Collective
3340: Input Parameter:
3341: + x - the vector
3342: . m - first dimension of four dimensional array
3343: . n - second dimension of four dimensional array
3344: . p - third dimension of four dimensional array
3345: . q - fourth dimension of four dimensional array
3346: . mstart - first index you will use in first coordinate direction (often 0)
3347: . nstart - first index in the second coordinate direction (often 0)
3348: . pstart - first index in the third coordinate direction (often 0)
3349: - qstart - first index in the fourth coordinate direction (often 0)
3351: Output Parameter:
3352: . a - location to put pointer to the array
3354: Level: beginner
3356: Notes:
3357: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3358: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3359: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3360: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
3362: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3364: Concepts: vector^accessing local values as 3d array
3366: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3367: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3368: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3369: @*/
3370: PetscErrorCode VecGetArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3371: {
3373: PetscInt i,N,j,k;
3374: PetscScalar *aa,***b,**c;
3380: VecGetLocalSize(x,&N);
3381: 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);
3382: VecGetArrayWrite(x,&aa);
3384: PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3385: b = (PetscScalar***)((*a) + m);
3386: c = (PetscScalar**)(b + m*n);
3387: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3388: for (i=0; i<m; i++)
3389: for (j=0; j<n; j++)
3390: b[i*n+j] = c + i*n*p + j*p - pstart;
3391: for (i=0; i<m; i++)
3392: for (j=0; j<n; j++)
3393: for (k=0; k<p; k++)
3394: c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
3395: *a -= mstart;
3396: return(0);
3397: }
3399: /*@C
3400: VecRestoreArray4d - Restores a vector after VecGetArray3d() has been called.
3402: Logically Collective
3404: Input Parameters:
3405: + x - the vector
3406: . m - first dimension of four dimensional array
3407: . n - second dimension of the four dimensional array
3408: . p - third dimension of the four dimensional array
3409: . q - fourth dimension of the four dimensional array
3410: . mstart - first index you will use in first coordinate direction (often 0)
3411: . nstart - first index in the second coordinate direction (often 0)
3412: . pstart - first index in the third coordinate direction (often 0)
3413: . qstart - first index in the fourth coordinate direction (often 0)
3414: - a - location of pointer to array obtained from VecGetArray4d()
3416: Level: beginner
3418: Notes:
3419: For regular PETSc vectors this routine does not involve any copies. For
3420: any special vectors that do not store local vector data in a contiguous
3421: array, this routine will copy the data back into the underlying
3422: vector data structure from the array obtained with VecGetArray().
3424: This routine actually zeros out the a pointer.
3426: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3427: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3428: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3429: @*/
3430: PetscErrorCode VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3431: {
3433: void *dummy;
3439: dummy = (void*)(*a + mstart);
3440: PetscFree(dummy);
3441: VecRestoreArray(x,NULL);
3442: return(0);
3443: }
3445: /*@C
3446: VecRestoreArray4dWrite - Restores a vector after VecGetArray3dWrite() has been called.
3448: Logically Collective
3450: Input Parameters:
3451: + x - the vector
3452: . m - first dimension of four dimensional array
3453: . n - second dimension of the four dimensional array
3454: . p - third dimension of the four dimensional array
3455: . q - fourth dimension of the four dimensional array
3456: . mstart - first index you will use in first coordinate direction (often 0)
3457: . nstart - first index in the second coordinate direction (often 0)
3458: . pstart - first index in the third coordinate direction (often 0)
3459: . qstart - first index in the fourth coordinate direction (often 0)
3460: - a - location of pointer to array obtained from VecGetArray4d()
3462: Level: beginner
3464: Notes:
3465: For regular PETSc vectors this routine does not involve any copies. For
3466: any special vectors that do not store local vector data in a contiguous
3467: array, this routine will copy the data back into the underlying
3468: vector data structure from the array obtained with VecGetArray().
3470: This routine actually zeros out the a pointer.
3472: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3473: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3474: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3475: @*/
3476: PetscErrorCode VecRestoreArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3477: {
3479: void *dummy;
3485: dummy = (void*)(*a + mstart);
3486: PetscFree(dummy);
3487: VecRestoreArrayWrite(x,NULL);
3488: return(0);
3489: }
3491: /*@C
3492: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3493: processor's portion of the vector data. You MUST call VecRestoreArray2dRead()
3494: when you no longer need access to the array.
3496: Logically Collective
3498: Input Parameter:
3499: + x - the vector
3500: . m - first dimension of two dimensional array
3501: . n - second dimension of two dimensional array
3502: . mstart - first index you will use in first coordinate direction (often 0)
3503: - nstart - first index in the second coordinate direction (often 0)
3505: Output Parameter:
3506: . a - location to put pointer to the array
3508: Level: developer
3510: Notes:
3511: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
3512: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3513: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3514: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
3516: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3518: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3519: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3520: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3521: @*/
3522: PetscErrorCode VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3523: {
3524: PetscErrorCode ierr;
3525: PetscInt i,N;
3526: const PetscScalar *aa;
3532: VecGetLocalSize(x,&N);
3533: 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);
3534: VecGetArrayRead(x,&aa);
3536: PetscMalloc1(m,a);
3537: for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
3538: *a -= mstart;
3539: return(0);
3540: }
3542: /*@C
3543: VecRestoreArray2dRead - Restores a vector after VecGetArray2dRead() has been called.
3545: Logically Collective
3547: Input Parameters:
3548: + x - the vector
3549: . m - first dimension of two dimensional array
3550: . n - second dimension of the two dimensional array
3551: . mstart - first index you will use in first coordinate direction (often 0)
3552: . nstart - first index in the second coordinate direction (often 0)
3553: - a - location of pointer to array obtained from VecGetArray2d()
3555: Level: developer
3557: Notes:
3558: For regular PETSc vectors this routine does not involve any copies. For
3559: any special vectors that do not store local vector data in a contiguous
3560: array, this routine will copy the data back into the underlying
3561: vector data structure from the array obtained with VecGetArray().
3563: This routine actually zeros out the a pointer.
3565: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3566: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3567: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3568: @*/
3569: PetscErrorCode VecRestoreArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
3570: {
3572: void *dummy;
3578: dummy = (void*)(*a + mstart);
3579: PetscFree(dummy);
3580: VecRestoreArrayRead(x,NULL);
3581: return(0);
3582: }
3584: /*@C
3585: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3586: processor's portion of the vector data. You MUST call VecRestoreArray1dRead()
3587: when you no longer need access to the array.
3589: Logically Collective
3591: Input Parameter:
3592: + x - the vector
3593: . m - first dimension of two dimensional array
3594: - mstart - first index you will use in first coordinate direction (often 0)
3596: Output Parameter:
3597: . a - location to put pointer to the array
3599: Level: developer
3601: Notes:
3602: For a vector obtained from DMCreateLocalVector() mstart are likely
3603: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3604: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
3606: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3608: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3609: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3610: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3611: @*/
3612: PetscErrorCode VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3613: {
3615: PetscInt N;
3621: VecGetLocalSize(x,&N);
3622: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3623: VecGetArrayRead(x,(const PetscScalar**)a);
3624: *a -= mstart;
3625: return(0);
3626: }
3628: /*@C
3629: VecRestoreArray1dRead - Restores a vector after VecGetArray1dRead() has been called.
3631: Logically Collective
3633: Input Parameters:
3634: + x - the vector
3635: . m - first dimension of two dimensional array
3636: . mstart - first index you will use in first coordinate direction (often 0)
3637: - a - location of pointer to array obtained from VecGetArray21()
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 VecGetArray1dRead().
3647: This routine actually zeros out the a pointer.
3649: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3650: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3651: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3652: @*/
3653: PetscErrorCode VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3654: {
3660: VecRestoreArrayRead(x,NULL);
3661: return(0);
3662: }
3665: /*@C
3666: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3667: processor's portion of the vector data. You MUST call VecRestoreArray3dRead()
3668: when you no longer need access to the array.
3670: Logically Collective
3672: Input Parameter:
3673: + x - the vector
3674: . m - first dimension of three dimensional array
3675: . n - second dimension of three dimensional array
3676: . p - third dimension of three dimensional array
3677: . mstart - first index you will use in first coordinate direction (often 0)
3678: . nstart - first index in the second coordinate direction (often 0)
3679: - pstart - first index in the third coordinate direction (often 0)
3681: Output Parameter:
3682: . a - location to put pointer to the array
3684: Level: developer
3686: Notes:
3687: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3688: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3689: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3690: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3dRead().
3692: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3694: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3695: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3696: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3697: @*/
3698: PetscErrorCode VecGetArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3699: {
3700: PetscErrorCode ierr;
3701: PetscInt i,N,j;
3702: const PetscScalar *aa;
3703: PetscScalar **b;
3709: VecGetLocalSize(x,&N);
3710: 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);
3711: VecGetArrayRead(x,&aa);
3713: PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3714: b = (PetscScalar**)((*a) + m);
3715: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3716: for (i=0; i<m; i++)
3717: for (j=0; j<n; j++)
3718: b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;
3720: *a -= mstart;
3721: return(0);
3722: }
3724: /*@C
3725: VecRestoreArray3dRead - Restores a vector after VecGetArray3dRead() has been called.
3727: Logically Collective
3729: Input Parameters:
3730: + x - the vector
3731: . m - first dimension of three dimensional array
3732: . n - second dimension of the three dimensional array
3733: . p - third dimension of the three dimensional array
3734: . mstart - first index you will use in first coordinate direction (often 0)
3735: . nstart - first index in the second coordinate direction (often 0)
3736: . pstart - first index in the third coordinate direction (often 0)
3737: - a - location of pointer to array obtained from VecGetArray3dRead()
3739: Level: developer
3741: Notes:
3742: For regular PETSc vectors this routine does not involve any copies. For
3743: any special vectors that do not store local vector data in a contiguous
3744: array, this routine will copy the data back into the underlying
3745: vector data structure from the array obtained with VecGetArray().
3747: This routine actually zeros out the a pointer.
3749: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3750: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3751: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3752: @*/
3753: PetscErrorCode VecRestoreArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3754: {
3756: void *dummy;
3762: dummy = (void*)(*a + mstart);
3763: PetscFree(dummy);
3764: VecRestoreArrayRead(x,NULL);
3765: return(0);
3766: }
3768: /*@C
3769: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3770: processor's portion of the vector data. You MUST call VecRestoreArray4dRead()
3771: when you no longer need access to the array.
3773: Logically Collective
3775: Input Parameter:
3776: + x - the vector
3777: . m - first dimension of four dimensional array
3778: . n - second dimension of four dimensional array
3779: . p - third dimension of four dimensional array
3780: . q - fourth dimension of four dimensional array
3781: . mstart - first index you will use in first coordinate direction (often 0)
3782: . nstart - first index in the second coordinate direction (often 0)
3783: . pstart - first index in the third coordinate direction (often 0)
3784: - qstart - first index in the fourth coordinate direction (often 0)
3786: Output Parameter:
3787: . a - location to put pointer to the array
3789: Level: beginner
3791: Notes:
3792: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3793: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3794: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3795: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
3797: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3799: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3800: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3801: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3802: @*/
3803: PetscErrorCode VecGetArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3804: {
3805: PetscErrorCode ierr;
3806: PetscInt i,N,j,k;
3807: const PetscScalar *aa;
3808: PetscScalar ***b,**c;
3814: VecGetLocalSize(x,&N);
3815: 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);
3816: VecGetArrayRead(x,&aa);
3818: PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3819: b = (PetscScalar***)((*a) + m);
3820: c = (PetscScalar**)(b + m*n);
3821: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3822: for (i=0; i<m; i++)
3823: for (j=0; j<n; j++)
3824: b[i*n+j] = c + i*n*p + j*p - pstart;
3825: for (i=0; i<m; i++)
3826: for (j=0; j<n; j++)
3827: for (k=0; k<p; k++)
3828: c[i*n*p+j*p+k] = (PetscScalar*) aa + i*n*p*q + j*p*q + k*q - qstart;
3829: *a -= mstart;
3830: return(0);
3831: }
3833: /*@C
3834: VecRestoreArray4dRead - Restores a vector after VecGetArray3d() has been called.
3836: Logically Collective
3838: Input Parameters:
3839: + x - the vector
3840: . m - first dimension of four dimensional array
3841: . n - second dimension of the four dimensional array
3842: . p - third dimension of the four dimensional array
3843: . q - fourth dimension of the four dimensional array
3844: . mstart - first index you will use in first coordinate direction (often 0)
3845: . nstart - first index in the second coordinate direction (often 0)
3846: . pstart - first index in the third coordinate direction (often 0)
3847: . qstart - first index in the fourth coordinate direction (often 0)
3848: - a - location of pointer to array obtained from VecGetArray4dRead()
3850: Level: beginner
3852: Notes:
3853: For regular PETSc vectors this routine does not involve any copies. For
3854: any special vectors that do not store local vector data in a contiguous
3855: array, this routine will copy the data back into the underlying
3856: vector data structure from the array obtained with VecGetArray().
3858: This routine actually zeros out the a pointer.
3860: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3861: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3862: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3863: @*/
3864: PetscErrorCode VecRestoreArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3865: {
3867: void *dummy;
3873: dummy = (void*)(*a + mstart);
3874: PetscFree(dummy);
3875: VecRestoreArrayRead(x,NULL);
3876: return(0);
3877: }
3879: #if defined(PETSC_USE_DEBUG)
3881: /*@
3882: VecLockGet - Gets the current lock status of a vector
3884: Logically Collective on Vec
3886: Input Parameter:
3887: . x - the vector
3889: Output Parameter:
3890: . state - greater than zero indicates the vector is locked for read; less then zero indicates the vector is
3891: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3893: Level: beginner
3895: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop()
3896: @*/
3897: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
3898: {
3901: *state = x->lock;
3902: return(0);
3903: }
3905: /*@
3906: VecLockReadPush - Pushes a read-only lock on a vector to prevent it from writing
3908: Logically Collective on Vec
3910: Input Parameter:
3911: . x - the vector
3913: Notes:
3914: If this is set then calls to VecGetArray() or VecSetValues() or any other routines that change the vectors values will fail.
3916: The call can be nested, i.e., called multiple times on the same vector, but each VecLockReadPush(x) has to have one matching
3917: VecLockReadPop(x), which removes the latest read-only lock.
3919: Level: beginner
3921: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPop(), VecLockGet()
3922: @*/
3923: PetscErrorCode VecLockReadPush(Vec x)
3924: {
3927: 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");
3928: x->lock++;
3929: return(0);
3930: }
3932: /*@
3933: VecLockReadPop - Pops a read-only lock from a vector
3935: Logically Collective on Vec
3937: Input Parameter:
3938: . x - the vector
3940: Level: beginner
3942: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockGet()
3943: @*/
3944: PetscErrorCode VecLockReadPop(Vec x)
3945: {
3948: x->lock--;
3949: if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector has been unlocked from read-only access too many times");
3950: return(0);
3951: }
3953: /*@C
3954: VecLockWriteSet_Private - Lock or unlock a vector for exclusive read/write access
3956: Logically Collective on Vec
3958: Input Parameter:
3959: + x - the vector
3960: - flg - PETSC_TRUE to lock the vector for writing; PETSC_FALSE to unlock it.
3962: Notes:
3963: The function is usefull in split-phase computations, which usually have a begin phase and an end phase.
3964: One can call VecLockWriteSet_Private(x,PETSC_TRUE) in the begin phase to lock a vector for exclusive
3965: access, and call VecLockWriteSet_Private(x,PETSC_FALSE) in the end phase to unlock the vector from exclusive
3966: access. In this way, one is ensured no other operations can access the vector in between. The code may like
3969: VecGetArray(x,&xdata); // begin phase
3970: VecLockWriteSet_Private(v,PETSC_TRUE);
3972: Other operations, which can not acceess x anymore (they can access xdata, of course)
3974: VecRestoreArray(x,&vdata); // end phase
3975: VecLockWriteSet_Private(v,PETSC_FALSE);
3977: The call can not be nested on the same vector, in other words, one can not call VecLockWriteSet_Private(x,PETSC_TRUE)
3978: again before calling VecLockWriteSet_Private(v,PETSC_FALSE).
3980: Level: beginner
3982: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop(), VecLockGet()
3983: @*/
3984: PetscErrorCode VecLockWriteSet_Private(Vec x,PetscBool flg)
3985: {
3988: if (flg) {
3989: 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");
3990: 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");
3991: else x->lock = -1;
3992: } else {
3993: 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");
3994: x->lock = 0;
3995: }
3996: return(0);
3997: }
3999: /*@
4000: VecLockPush - Pushes a read-only lock on a vector to prevent it from writing
4002: Level: deprecated
4004: .seealso: VecLockReadPush()
4005: @*/
4006: PetscErrorCode VecLockPush(Vec x)
4007: {
4010: VecLockReadPush(x);
4011: return(0);
4012: }
4014: /*@
4015: VecLockPop - Pops a read-only lock from a vector
4017: Level: deprecated
4019: .seealso: VecLockReadPop()
4020: @*/
4021: PetscErrorCode VecLockPop(Vec x)
4022: {
4025: VecLockReadPop(x);
4026: return(0);
4027: }
4029: #endif