Actual source code: rvector.c
petsc-master 2019-08-18
2: /*
3: Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
4: These are the vector functions the user calls.
5: */
6: #include <petsc/private/vecimpl.h>
7: static PetscInt VecGetSubVectorSavedStateId = -1;
9: PETSC_EXTERN PetscErrorCode VecValidValues(Vec vec,PetscInt argnum,PetscBool begin)
10: {
11: #if defined(PETSC_USE_DEBUG)
12: PetscErrorCode ierr;
13: PetscInt n,i;
14: const PetscScalar *x;
17: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL)
18: if ((vec->petscnative || vec->ops->getarray) && (vec->valid_GPU_array == PETSC_OFFLOAD_CPU || vec->valid_GPU_array == PETSC_OFFLOAD_BOTH)) {
19: #else
20: if (vec->petscnative || vec->ops->getarray) {
21: #endif
22: VecGetLocalSize(vec,&n);
23: VecGetArrayRead(vec,&x);
24: for (i=0; i<n; i++) {
25: if (begin) {
26: if (PetscIsInfOrNanScalar(x[i])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FP,"Vec entry at local location %D is not-a-number or infinite at beginning of function: Parameter number %D",i,argnum);
27: } else {
28: if (PetscIsInfOrNanScalar(x[i])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FP,"Vec entry at local location %D is not-a-number or infinite at end of function: Parameter number %D",i,argnum);
29: }
30: }
31: VecRestoreArrayRead(vec,&x);
32: }
33: return(0);
34: #else
35: return 0;
36: #endif
37: }
39: /*@
40: VecMaxPointwiseDivide - Computes the maximum of the componentwise division max = max_i abs(x_i/y_i).
42: Logically Collective on Vec
44: Input Parameters:
45: . x, y - the vectors
47: Output Parameter:
48: . max - the result
50: Level: advanced
52: Notes:
53: x and y may be the same vector
54: if a particular y_i is zero, it is treated as 1 in the above formula
56: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs()
57: @*/
58: PetscErrorCode VecMaxPointwiseDivide(Vec x,Vec y,PetscReal *max)
59: {
69: VecCheckSameSize(x,1,y,2);
70: (*x->ops->maxpointwisedivide)(x,y,max);
71: return(0);
72: }
74: /*@
75: VecDot - Computes the vector dot product.
77: Collective on Vec
79: Input Parameters:
80: . x, y - the vectors
82: Output Parameter:
83: . val - the dot product
85: Performance Issues:
86: $ per-processor memory bandwidth
87: $ interprocessor latency
88: $ work load inbalance that causes certain processes to arrive much earlier than others
90: Notes for Users of Complex Numbers:
91: For complex vectors, VecDot() computes
92: $ val = (x,y) = y^H x,
93: where y^H denotes the conjugate transpose of y. Note that this corresponds to the usual "mathematicians" complex
94: inner product where the SECOND argument gets the complex conjugate. Since the BLASdot() complex conjugates the first
95: first argument we call the BLASdot() with the arguments reversed.
97: Use VecTDot() for the indefinite form
98: $ val = (x,y) = y^T x,
99: where y^T denotes the transpose of y.
101: Level: intermediate
104: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDotRealPart()
105: @*/
106: PetscErrorCode VecDot(Vec x,Vec y,PetscScalar *val)
107: {
117: VecCheckSameSize(x,1,y,2);
119: PetscLogEventBegin(VEC_Dot,x,y,0,0);
120: (*x->ops->dot)(x,y,val);
121: PetscLogEventEnd(VEC_Dot,x,y,0,0);
122: return(0);
123: }
125: /*@
126: VecDotRealPart - Computes the real part of the vector dot product.
128: Collective on Vec
130: Input Parameters:
131: . x, y - the vectors
133: Output Parameter:
134: . val - the real part of the dot product;
136: Performance Issues:
137: $ per-processor memory bandwidth
138: $ interprocessor latency
139: $ work load inbalance that causes certain processes to arrive much earlier than others
141: Notes for Users of Complex Numbers:
142: See VecDot() for more details on the definition of the dot product for complex numbers
144: For real numbers this returns the same value as VecDot()
146: 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
147: the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)
149: Developer Note: This is not currently optimized to compute only the real part of the dot product.
151: Level: intermediate
154: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDot(), VecDotNorm2()
155: @*/
156: PetscErrorCode VecDotRealPart(Vec x,Vec y,PetscReal *val)
157: {
159: PetscScalar fdot;
162: VecDot(x,y,&fdot);
163: *val = PetscRealPart(fdot);
164: return(0);
165: }
167: /*@
168: VecNorm - Computes the vector norm.
170: Collective on Vec
172: Input Parameters:
173: + x - the vector
174: - type - one of NORM_1, NORM_2, NORM_INFINITY. Also available
175: NORM_1_AND_2, which computes both norms and stores them
176: in a two element array.
178: Output Parameter:
179: . val - the norm
181: Notes:
182: $ NORM_1 denotes sum_i |x_i|
183: $ NORM_2 denotes sqrt(sum_i |x_i|^2)
184: $ NORM_INFINITY denotes max_i |x_i|
186: For complex numbers NORM_1 will return the traditional 1 norm of the 2 norm of the complex numbers; that is the 1
187: norm of the absolutely values of the complex entries. In PETSc 3.6 and earlier releases it returned the 1 norm of
188: the 1 norm of the complex entries (what is returned by the BLAS routine asum()). Both are valid norms but most
189: people expect the former.
191: Level: intermediate
193: Performance Issues:
194: $ per-processor memory bandwidth
195: $ interprocessor latency
196: $ work load inbalance that causes certain processes to arrive much earlier than others
199: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNormAvailable(),
200: VecNormBegin(), VecNormEnd()
202: @*/
204: PetscErrorCode VecNorm(Vec x,NormType type,PetscReal *val)
205: {
206: PetscBool flg;
214: /*
215: * Cached data?
216: */
217: if (type!=NORM_1_AND_2) {
218: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,flg);
219: if (flg) return(0);
220: }
221: PetscLogEventBegin(VEC_Norm,x,0,0,0);
222: (*x->ops->norm)(x,type,val);
223: PetscLogEventEnd(VEC_Norm,x,0,0,0);
224: if (type!=NORM_1_AND_2) {
225: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[type],*val);
226: }
227: return(0);
228: }
230: /*@
231: VecNormAvailable - Returns the vector norm if it is already known.
233: Not Collective
235: Input Parameters:
236: + x - the vector
237: - type - one of NORM_1, NORM_2, NORM_INFINITY. Also available
238: NORM_1_AND_2, which computes both norms and stores them
239: in a two element array.
241: Output Parameter:
242: + available - PETSC_TRUE if the val returned is valid
243: - val - the norm
245: Notes:
246: $ NORM_1 denotes sum_i |x_i|
247: $ NORM_2 denotes sqrt(sum_i (x_i)^2)
248: $ NORM_INFINITY denotes max_i |x_i|
250: Level: intermediate
252: Performance Issues:
253: $ per-processor memory bandwidth
254: $ interprocessor latency
255: $ work load inbalance that causes certain processes to arrive much earlier than others
257: Compile Option:
258: PETSC_HAVE_SLOW_BLAS_NORM2 will cause a C (loop unrolled) version of the norm to be used, rather
259: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
260: (as opposed to vendor provided) because the FORTRAN BLAS NRM2() routine is very slow.
263: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNorm()
264: VecNormBegin(), VecNormEnd()
266: @*/
267: PetscErrorCode VecNormAvailable(Vec x,NormType type,PetscBool *available,PetscReal *val)
268: {
276: *available = PETSC_FALSE;
277: if (type!=NORM_1_AND_2) {
278: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,*available);
279: }
280: return(0);
281: }
283: /*@
284: VecNormalize - Normalizes a vector by 2-norm.
286: Collective on Vec
288: Input Parameters:
289: + x - the vector
291: Output Parameter:
292: . x - the normalized vector
293: - val - the vector norm before normalization
295: Level: intermediate
298: @*/
299: PetscErrorCode VecNormalize(Vec x,PetscReal *val)
300: {
302: PetscReal norm;
307: PetscLogEventBegin(VEC_Normalize,x,0,0,0);
308: VecNorm(x,NORM_2,&norm);
309: if (norm == 0.0) {
310: PetscInfo(x,"Vector of zero norm can not be normalized; Returning only the zero norm\n");
311: } else if (norm != 1.0) {
312: PetscScalar tmp = 1.0/norm;
313: VecScale(x,tmp);
314: }
315: if (val) *val = norm;
316: PetscLogEventEnd(VEC_Normalize,x,0,0,0);
317: return(0);
318: }
320: /*@C
321: VecMax - Determines the vector component with maximum real part and its location.
323: Collective on Vec
325: Input Parameter:
326: . x - the vector
328: Output Parameters:
329: + p - the location of val (pass NULL if you don't want this)
330: - val - the maximum component
332: Notes:
333: Returns the value PETSC_MIN_REAL and p = -1 if the vector is of length 0.
335: Returns the smallest index with the maximum value
336: Level: intermediate
339: .seealso: VecNorm(), VecMin()
340: @*/
341: PetscErrorCode VecMax(Vec x,PetscInt *p,PetscReal *val)
342: {
349: PetscLogEventBegin(VEC_Max,x,0,0,0);
350: (*x->ops->max)(x,p,val);
351: PetscLogEventEnd(VEC_Max,x,0,0,0);
352: return(0);
353: }
355: /*@C
356: VecMin - Determines the vector component with minimum real part and its location.
358: Collective on Vec
360: Input Parameters:
361: . x - the vector
363: Output Parameter:
364: + p - the location of val (pass NULL if you don't want this location)
365: - val - the minimum component
367: Level: intermediate
369: Notes:
370: Returns the value PETSC_MAX_REAL and p = -1 if the vector is of length 0.
372: This returns the smallest index with the minumum value
375: .seealso: VecMax()
376: @*/
377: PetscErrorCode VecMin(Vec x,PetscInt *p,PetscReal *val)
378: {
385: PetscLogEventBegin(VEC_Min,x,0,0,0);
386: (*x->ops->min)(x,p,val);
387: PetscLogEventEnd(VEC_Min,x,0,0,0);
388: return(0);
389: }
391: /*@
392: VecTDot - Computes an indefinite vector dot product. That is, this
393: routine does NOT use the complex conjugate.
395: Collective on Vec
397: Input Parameters:
398: . x, y - the vectors
400: Output Parameter:
401: . val - the dot product
403: Notes for Users of Complex Numbers:
404: For complex vectors, VecTDot() computes the indefinite form
405: $ val = (x,y) = y^T x,
406: where y^T denotes the transpose of y.
408: Use VecDot() for the inner product
409: $ val = (x,y) = y^H x,
410: where y^H denotes the conjugate transpose of y.
412: Level: intermediate
414: .seealso: VecDot(), VecMTDot()
415: @*/
416: PetscErrorCode VecTDot(Vec x,Vec y,PetscScalar *val)
417: {
427: VecCheckSameSize(x,1,y,2);
429: PetscLogEventBegin(VEC_TDot,x,y,0,0);
430: (*x->ops->tdot)(x,y,val);
431: PetscLogEventEnd(VEC_TDot,x,y,0,0);
432: return(0);
433: }
435: /*@
436: VecScale - Scales a vector.
438: Not collective on Vec
440: Input Parameters:
441: + x - the vector
442: - alpha - the scalar
444: Output Parameter:
445: . x - the scaled vector
447: Note:
448: For a vector with n components, VecScale() computes
449: $ x[i] = alpha * x[i], for i=1,...,n.
451: Level: intermediate
454: @*/
455: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
456: {
457: PetscReal norms[4] = {0.0,0.0,0.0, 0.0};
458: PetscBool flgs[4];
460: PetscInt i;
465: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
466: PetscLogEventBegin(VEC_Scale,x,0,0,0);
467: if (alpha != (PetscScalar)1.0) {
468: /* get current stashed norms */
469: for (i=0; i<4; i++) {
470: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
471: }
472: (*x->ops->scale)(x,alpha);
473: PetscObjectStateIncrease((PetscObject)x);
474: /* put the scaled stashed norms back into the Vec */
475: for (i=0; i<4; i++) {
476: if (flgs[i]) {
477: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],PetscAbsScalar(alpha)*norms[i]);
478: }
479: }
480: }
481: PetscLogEventEnd(VEC_Scale,x,0,0,0);
482: return(0);
483: }
485: /*@
486: VecSet - Sets all components of a vector to a single scalar value.
488: Logically Collective on Vec
490: Input Parameters:
491: + x - the vector
492: - alpha - the scalar
494: Output Parameter:
495: . x - the vector
497: Note:
498: For a vector of dimension n, VecSet() computes
499: $ x[i] = alpha, for i=1,...,n,
500: so that all vector entries then equal the identical
501: scalar value, alpha. Use the more general routine
502: VecSetValues() to set different vector entries.
504: You CANNOT call this after you have called VecSetValues() but before you call
505: VecAssemblyBegin/End().
507: Level: beginner
509: .seealso VecSetValues(), VecSetValuesBlocked(), VecSetRandom()
511: @*/
512: PetscErrorCode VecSet(Vec x,PetscScalar alpha)
513: {
514: PetscReal val;
520: 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()");
522: VecSetErrorIfLocked(x,1);
524: PetscLogEventBegin(VEC_Set,x,0,0,0);
525: (*x->ops->set)(x,alpha);
526: PetscLogEventEnd(VEC_Set,x,0,0,0);
527: PetscObjectStateIncrease((PetscObject)x);
529: /* norms can be simply set (if |alpha|*N not too large) */
530: val = PetscAbsScalar(alpha);
531: if (x->map->N == 0) {
532: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],0.0l);
533: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],0.0);
534: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],0.0);
535: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],0.0);
536: } else if (val > PETSC_MAX_REAL/x->map->N) {
537: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
538: } else {
539: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],x->map->N * val);
540: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
541: val = PetscSqrtReal((PetscReal)x->map->N) * val;
542: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],val);
543: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],val);
544: }
545: return(0);
546: }
549: /*@
550: VecAXPY - Computes y = alpha x + y.
552: Logically Collective on Vec
554: Input Parameters:
555: + alpha - the scalar
556: - x, y - the vectors
558: Output Parameter:
559: . y - output vector
561: Level: intermediate
563: Notes:
564: x and y MUST be different vectors
565: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
567: $ VecAXPY(y,alpha,x) y = alpha x + y
568: $ VecAYPX(y,beta,x) y = x + beta y
569: $ VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
570: $ VecWAXPY(w,alpha,x,y) w = alpha x + y
571: $ VecAXPBYPCZ(w,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
572: $ VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
575: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPBYPCZ(), VecAXPBY()
576: @*/
577: PetscErrorCode VecAXPY(Vec y,PetscScalar alpha,Vec x)
578: {
587: VecCheckSameSize(x,1,y,3);
588: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
590: VecSetErrorIfLocked(y,1);
592: VecLockReadPush(x);
593: PetscLogEventBegin(VEC_AXPY,x,y,0,0);
594: (*y->ops->axpy)(y,alpha,x);
595: PetscLogEventEnd(VEC_AXPY,x,y,0,0);
596: VecLockReadPop(x);
597: PetscObjectStateIncrease((PetscObject)y);
598: return(0);
599: }
601: /*@
602: VecAXPBY - Computes y = alpha x + beta y.
604: Logically Collective on Vec
606: Input Parameters:
607: + alpha,beta - the scalars
608: - x, y - the vectors
610: Output Parameter:
611: . y - output vector
613: Level: intermediate
615: Notes:
616: x and y MUST be different vectors
617: The implementation is optimized for alpha and/or beta values of 0.0 and 1.0
620: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ()
621: @*/
622: PetscErrorCode VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
623: {
632: VecCheckSameSize(y,1,x,4);
633: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
636: PetscLogEventBegin(VEC_AXPY,x,y,0,0);
637: (*y->ops->axpby)(y,alpha,beta,x);
638: PetscLogEventEnd(VEC_AXPY,x,y,0,0);
639: PetscObjectStateIncrease((PetscObject)y);
640: return(0);
641: }
643: /*@
644: VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z
646: Logically Collective on Vec
648: Input Parameters:
649: + alpha,beta, gamma - the scalars
650: - x, y, z - the vectors
652: Output Parameter:
653: . z - output vector
655: Level: intermediate
657: Notes:
658: x, y and z must be different vectors
659: The implementation is optimized for alpha of 1.0 and gamma of 1.0 or 0.0
662: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBY()
663: @*/
664: PetscErrorCode VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
665: {
677: VecCheckSameSize(x,1,y,5);
678: VecCheckSameSize(x,1,z,6);
679: if (x == y || x == z) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
680: if (y == z) SETERRQ(PetscObjectComm((PetscObject)y),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
685: PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
686: (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
687: PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
688: PetscObjectStateIncrease((PetscObject)z);
689: return(0);
690: }
692: /*@
693: VecAYPX - Computes y = x + beta y.
695: Logically Collective on Vec
697: Input Parameters:
698: + beta - the scalar
699: - x, y - the vectors
701: Output Parameter:
702: . y - output vector
704: Level: intermediate
706: Notes:
707: x and y MUST be different vectors
708: The implementation is optimized for beta of -1.0, 0.0, and 1.0
711: .seealso: VecMAXPY(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
712: @*/
713: PetscErrorCode VecAYPX(Vec y,PetscScalar beta,Vec x)
714: {
723: VecCheckSameSize(x,1,y,3);
724: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
727: PetscLogEventBegin(VEC_AYPX,x,y,0,0);
728: (*y->ops->aypx)(y,beta,x);
729: PetscLogEventEnd(VEC_AYPX,x,y,0,0);
730: PetscObjectStateIncrease((PetscObject)y);
731: return(0);
732: }
735: /*@
736: VecWAXPY - Computes w = alpha x + y.
738: Logically Collective on Vec
740: Input Parameters:
741: + alpha - the scalar
742: - x, y - the vectors
744: Output Parameter:
745: . w - the result
747: Level: intermediate
749: Notes:
750: w cannot be either x or y, but x and y can be the same
751: The implementation is optimzed for alpha of -1.0, 0.0, and 1.0
754: .seealso: VecAXPY(), VecAYPX(), VecAXPBY(), VecMAXPY(), VecAXPBYPCZ()
755: @*/
756: PetscErrorCode VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
757: {
769: VecCheckSameSize(x,3,y,4);
770: VecCheckSameSize(x,3,w,1);
771: if (w == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector y, suggest VecAXPY()");
772: if (w == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector x, suggest VecAYPX()");
775: PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
776: (*w->ops->waxpy)(w,alpha,x,y);
777: PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
778: PetscObjectStateIncrease((PetscObject)w);
779: return(0);
780: }
783: /*@C
784: VecSetValues - Inserts or adds values into certain locations of a vector.
786: Not Collective
788: Input Parameters:
789: + x - vector to insert in
790: . ni - number of elements to add
791: . ix - indices where to add
792: . y - array of values
793: - iora - either INSERT_VALUES or ADD_VALUES, where
794: ADD_VALUES adds values to any existing entries, and
795: INSERT_VALUES replaces existing entries with new values
797: Notes:
798: VecSetValues() sets x[ix[i]] = y[i], for i=0,...,ni-1.
800: Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
801: options cannot be mixed without intervening calls to the assembly
802: routines.
804: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
805: MUST be called after all calls to VecSetValues() have been completed.
807: VecSetValues() uses 0-based indices in Fortran as well as in C.
809: If you call VecSetOption(x, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE),
810: negative indices may be passed in ix. These rows are
811: simply ignored. This allows easily inserting element load matrices
812: with homogeneous Dirchlet boundary conditions that you don't want represented
813: in the vector.
815: Level: beginner
817: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
818: VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
819: @*/
820: PetscErrorCode VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
821: {
826: if (!ni) return(0);
830: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
831: (*x->ops->setvalues)(x,ni,ix,y,iora);
832: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
833: PetscObjectStateIncrease((PetscObject)x);
834: return(0);
835: }
837: /*@
838: VecGetValues - Gets values from certain locations of a vector. Currently
839: can only get values on the same processor
841: Not Collective
843: Input Parameters:
844: + x - vector to get values from
845: . ni - number of elements to get
846: - ix - indices where to get them from (in global 1d numbering)
848: Output Parameter:
849: . y - array of values
851: Notes:
852: The user provides the allocated array y; it is NOT allocated in this routine
854: VecGetValues() gets y[i] = x[ix[i]], for i=0,...,ni-1.
856: VecAssemblyBegin() and VecAssemblyEnd() MUST be called before calling this
858: VecGetValues() uses 0-based indices in Fortran as well as in C.
860: If you call VecSetOption(x, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE),
861: negative indices may be passed in ix. These rows are
862: simply ignored.
864: Level: beginner
866: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues()
867: @*/
868: PetscErrorCode VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
869: {
874: if (!ni) return(0);
878: (*x->ops->getvalues)(x,ni,ix,y);
879: return(0);
880: }
882: /*@C
883: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
885: Not Collective
887: Input Parameters:
888: + x - vector to insert in
889: . ni - number of blocks to add
890: . ix - indices where to add in block count, rather than element count
891: . y - array of values
892: - iora - either INSERT_VALUES or ADD_VALUES, where
893: ADD_VALUES adds values to any existing entries, and
894: INSERT_VALUES replaces existing entries with new values
896: Notes:
897: VecSetValuesBlocked() sets x[bs*ix[i]+j] = y[bs*i+j],
898: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
900: Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
901: options cannot be mixed without intervening calls to the assembly
902: routines.
904: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
905: MUST be called after all calls to VecSetValuesBlocked() have been completed.
907: VecSetValuesBlocked() uses 0-based indices in Fortran as well as in C.
909: Negative indices may be passed in ix, these rows are
910: simply ignored. This allows easily inserting element load matrices
911: with homogeneous Dirchlet boundary conditions that you don't want represented
912: in the vector.
914: Level: intermediate
916: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
917: VecSetValues()
918: @*/
919: PetscErrorCode VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
920: {
925: if (!ni) return(0);
929: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
930: (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
931: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
932: PetscObjectStateIncrease((PetscObject)x);
933: return(0);
934: }
937: /*@C
938: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
939: using a local ordering of the nodes.
941: Not Collective
943: Input Parameters:
944: + x - vector to insert in
945: . ni - number of elements to add
946: . ix - indices where to add
947: . y - array of values
948: - iora - either INSERT_VALUES or ADD_VALUES, where
949: ADD_VALUES adds values to any existing entries, and
950: INSERT_VALUES replaces existing entries with new values
952: Level: intermediate
954: Notes:
955: VecSetValuesLocal() sets x[ix[i]] = y[i], for i=0,...,ni-1.
957: Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
958: options cannot be mixed without intervening calls to the assembly
959: routines.
961: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
962: MUST be called after all calls to VecSetValuesLocal() have been completed.
964: VecSetValuesLocal() uses 0-based indices in Fortran as well as in C.
966: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetLocalToGlobalMapping(),
967: VecSetValuesBlockedLocal()
968: @*/
969: PetscErrorCode VecSetValuesLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
970: {
972: PetscInt lixp[128],*lix = lixp;
976: if (!ni) return(0);
981: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
982: if (!x->ops->setvalueslocal) {
983: if (!x->map->mapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMapping()");
984: if (ni > 128) {
985: PetscMalloc1(ni,&lix);
986: }
987: ISLocalToGlobalMappingApply(x->map->mapping,ni,(PetscInt*)ix,lix);
988: (*x->ops->setvalues)(x,ni,lix,y,iora);
989: if (ni > 128) {
990: PetscFree(lix);
991: }
992: } else {
993: (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
994: }
995: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
996: PetscObjectStateIncrease((PetscObject)x);
997: return(0);
998: }
1000: /*@
1001: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1002: using a local ordering of the nodes.
1004: Not Collective
1006: Input Parameters:
1007: + x - vector to insert in
1008: . ni - number of blocks to add
1009: . ix - indices where to add in block count, not element count
1010: . y - array of values
1011: - iora - either INSERT_VALUES or ADD_VALUES, where
1012: ADD_VALUES adds values to any existing entries, and
1013: INSERT_VALUES replaces existing entries with new values
1015: Level: intermediate
1017: Notes:
1018: VecSetValuesBlockedLocal() sets x[bs*ix[i]+j] = y[bs*i+j],
1019: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with VecSetBlockSize().
1021: Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1022: options cannot be mixed without intervening calls to the assembly
1023: routines.
1025: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
1026: MUST be called after all calls to VecSetValuesBlockedLocal() have been completed.
1028: VecSetValuesBlockedLocal() uses 0-based indices in Fortran as well as in C.
1031: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
1032: VecSetLocalToGlobalMapping()
1033: @*/
1034: PetscErrorCode VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1035: {
1037: PetscInt lixp[128],*lix = lixp;
1041: if (!ni) return(0);
1045: if (ni > 128) {
1046: PetscMalloc1(ni,&lix);
1047: }
1049: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1050: ISLocalToGlobalMappingApplyBlock(x->map->mapping,ni,(PetscInt*)ix,lix);
1051: (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1052: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1053: if (ni > 128) {
1054: PetscFree(lix);
1055: }
1056: PetscObjectStateIncrease((PetscObject)x);
1057: return(0);
1058: }
1060: /*@
1061: VecMTDot - Computes indefinite vector multiple dot products.
1062: That is, it does NOT use the complex conjugate.
1064: Collective on Vec
1066: Input Parameters:
1067: + x - one vector
1068: . nv - number of vectors
1069: - y - array of vectors. Note that vectors are pointers
1071: Output Parameter:
1072: . val - array of the dot products
1074: Notes for Users of Complex Numbers:
1075: For complex vectors, VecMTDot() computes the indefinite form
1076: $ val = (x,y) = y^T x,
1077: where y^T denotes the transpose of y.
1079: Use VecMDot() for the inner product
1080: $ val = (x,y) = y^H x,
1081: where y^H denotes the conjugate transpose of y.
1083: Level: intermediate
1086: .seealso: VecMDot(), VecTDot()
1087: @*/
1088: PetscErrorCode VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1089: {
1095: if (!nv) return(0);
1102: VecCheckSameSize(x,1,*y,3);
1104: PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1105: (*x->ops->mtdot)(x,nv,y,val);
1106: PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1107: return(0);
1108: }
1110: /*@
1111: VecMDot - Computes vector multiple dot products.
1113: Collective on Vec
1115: Input Parameters:
1116: + x - one vector
1117: . nv - number of vectors
1118: - y - array of vectors.
1120: Output Parameter:
1121: . val - array of the dot products (does not allocate the array)
1123: Notes for Users of Complex Numbers:
1124: For complex vectors, VecMDot() computes
1125: $ val = (x,y) = y^H x,
1126: where y^H denotes the conjugate transpose of y.
1128: Use VecMTDot() for the indefinite form
1129: $ val = (x,y) = y^T x,
1130: where y^T denotes the transpose of y.
1132: Level: intermediate
1135: .seealso: VecMTDot(), VecDot()
1136: @*/
1137: PetscErrorCode VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1138: {
1144: if (!nv) return(0);
1145: if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1152: VecCheckSameSize(x,1,*y,3);
1154: PetscLogEventBegin(VEC_MDot,x,*y,0,0);
1155: (*x->ops->mdot)(x,nv,y,val);
1156: PetscLogEventEnd(VEC_MDot,x,*y,0,0);
1157: return(0);
1158: }
1160: /*@
1161: VecMAXPY - Computes y = y + sum alpha[i] x[i]
1163: Logically Collective on Vec
1165: Input Parameters:
1166: + nv - number of scalars and x-vectors
1167: . alpha - array of scalars
1168: . y - one vector
1169: - x - array of vectors
1171: Level: intermediate
1173: Notes:
1174: y cannot be any of the x vectors
1176: .seealso: VecAYPX(), VecWAXPY(), VecAXPY(), VecAXPBYPCZ(), VecAXPBY()
1177: @*/
1178: PetscErrorCode VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1179: {
1181: PetscInt i;
1186: if (!nv) return(0);
1187: if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1194: VecCheckSameSize(y,1,*x,4);
1197: PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1198: (*y->ops->maxpy)(y,nv,alpha,x);
1199: PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1200: PetscObjectStateIncrease((PetscObject)y);
1201: return(0);
1202: }
1204: /*@
1205: VecGetSubVector - Gets a vector representing part of another vector
1207: Collective on IS
1209: Input Arguments:
1210: + X - vector from which to extract a subvector
1211: - is - index set representing portion of X to extract
1213: Output Arguments:
1214: . Y - subvector corresponding to is
1216: Level: advanced
1218: Notes:
1219: The subvector Y should be returned with VecRestoreSubVector().
1221: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1222: modifying the subvector. Other non-overlapping subvectors can still be obtained from X using this function.
1224: .seealso: MatCreateSubMatrix()
1225: @*/
1226: PetscErrorCode VecGetSubVector(Vec X,IS is,Vec *Y)
1227: {
1228: PetscErrorCode ierr;
1229: Vec Z;
1235: if (X->ops->getsubvector) {
1236: (*X->ops->getsubvector)(X,is,&Z);
1237: } else { /* Default implementation currently does no caching */
1238: PetscInt gstart,gend,start;
1239: PetscBool contiguous,gcontiguous;
1240: VecGetOwnershipRange(X,&gstart,&gend);
1241: ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1242: MPIU_Allreduce(&contiguous,&gcontiguous,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1243: if (gcontiguous) { /* We can do a no-copy implementation */
1244: PetscInt n,N,bs;
1245: PetscInt state;
1247: ISGetSize(is,&N);
1248: ISGetLocalSize(is,&n);
1249: VecGetBlockSize(X,&bs);
1250: if (n%bs || bs == 1 || !n) bs = -1; /* Do not decide block size if we do not have to */
1251: VecLockGet(X,&state);
1252: if (state) {
1253: const PetscScalar *x;
1254: VecGetArrayRead(X,&x);
1255: VecCreate(PetscObjectComm((PetscObject)X),&Z);
1256: VecSetType(Z,((PetscObject)X)->type_name);
1257: VecSetSizes(Z,n,N);
1258: VecSetBlockSize(Z,bs);
1259: VecPlaceArray(Z,(PetscScalar*)x+start);
1260: VecLockReadPush(Z);
1261: VecRestoreArrayRead(X,&x);
1262: } else {
1263: PetscScalar *x;
1264: VecGetArray(X,&x);
1265: VecCreate(PetscObjectComm((PetscObject)X),&Z);
1266: VecSetType(Z,((PetscObject)X)->type_name);
1267: VecSetSizes(Z,n,N);
1268: VecSetBlockSize(Z,bs);
1269: VecPlaceArray(Z,(PetscScalar*)x+start);
1270: VecRestoreArray(X,&x);
1271: }
1272: } else { /* Have to create a scatter and do a copy */
1273: VecScatter scatter;
1274: PetscInt n,N;
1275: ISGetLocalSize(is,&n);
1276: ISGetSize(is,&N);
1277: VecCreate(PetscObjectComm((PetscObject)is),&Z);
1278: VecSetSizes(Z,n,N);
1279: VecSetType(Z,((PetscObject)X)->type_name);
1280: VecScatterCreate(X,is,Z,NULL,&scatter);
1281: VecScatterBegin(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1282: VecScatterEnd(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1283: PetscObjectCompose((PetscObject)Z,"VecGetSubVector_Scatter",(PetscObject)scatter);
1284: VecScatterDestroy(&scatter);
1285: }
1286: }
1287: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1288: if (VecGetSubVectorSavedStateId < 0) {PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);}
1289: PetscObjectComposedDataSetInt((PetscObject)Z,VecGetSubVectorSavedStateId,1);
1290: *Y = Z;
1291: return(0);
1292: }
1294: /*@
1295: VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()
1297: Collective on IS
1299: Input Arguments:
1300: + X - vector from which subvector was obtained
1301: . is - index set representing the subset of X
1302: - Y - subvector being restored
1304: Level: advanced
1306: .seealso: VecGetSubVector()
1307: @*/
1308: PetscErrorCode VecRestoreSubVector(Vec X,IS is,Vec *Y)
1309: {
1317: if (X->ops->restoresubvector) {
1318: (*X->ops->restoresubvector)(X,is,Y);
1319: } else {
1320: PETSC_UNUSED PetscObjectState dummystate = 0;
1321: PetscBool valid;
1322: PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,dummystate,valid);
1323: if (!valid) {
1324: VecScatter scatter;
1326: PetscObjectQuery((PetscObject)*Y,"VecGetSubVector_Scatter",(PetscObject*)&scatter);
1327: if (scatter) {
1328: VecScatterBegin(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1329: VecScatterEnd(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1330: }
1331: }
1332: VecDestroy(Y);
1333: }
1334: return(0);
1335: }
1337: /*@
1338: VecGetLocalVectorRead - Maps the local portion of a vector into a
1339: vector. You must call VecRestoreLocalVectorRead() when the local
1340: vector is no longer needed.
1342: Not collective.
1344: Input parameter:
1345: . v - The vector for which the local vector is desired.
1347: Output parameter:
1348: . w - Upon exit this contains the local vector.
1350: Level: beginner
1352: Notes:
1353: This function is similar to VecGetArrayRead() which maps the local
1354: portion into a raw pointer. VecGetLocalVectorRead() is usually
1355: almost as efficient as VecGetArrayRead() but in certain circumstances
1356: VecGetLocalVectorRead() can be much more efficient than
1357: VecGetArrayRead(). This is because the construction of a contiguous
1358: array representing the vector data required by VecGetArrayRead() can
1359: be an expensive operation for certain vector types. For example, for
1360: GPU vectors VecGetArrayRead() requires that the data between device
1361: and host is synchronized.
1363: Unlike VecGetLocalVector(), this routine is not collective and
1364: preserves cached information.
1366: .seealso: VecRestoreLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1367: @*/
1368: PetscErrorCode VecGetLocalVectorRead(Vec v,Vec w)
1369: {
1371: PetscScalar *a;
1376: VecCheckSameLocalSize(v,1,w,2);
1377: if (v->ops->getlocalvectorread) {
1378: (*v->ops->getlocalvectorread)(v,w);
1379: } else {
1380: VecGetArrayRead(v,(const PetscScalar**)&a);
1381: VecPlaceArray(w,a);
1382: }
1383: return(0);
1384: }
1386: /*@
1387: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1388: previously mapped into a vector using VecGetLocalVectorRead().
1390: Not collective.
1392: Input parameter:
1393: + v - The local portion of this vector was previously mapped into w using VecGetLocalVectorRead().
1394: - w - The vector into which the local portion of v was mapped.
1396: Level: beginner
1398: .seealso: VecGetLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1399: @*/
1400: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec w)
1401: {
1403: PetscScalar *a;
1408: if (v->ops->restorelocalvectorread) {
1409: (*v->ops->restorelocalvectorread)(v,w);
1410: } else {
1411: VecGetArrayRead(w,(const PetscScalar**)&a);
1412: VecRestoreArrayRead(v,(const PetscScalar**)&a);
1413: VecResetArray(w);
1414: }
1415: return(0);
1416: }
1418: /*@
1419: VecGetLocalVector - Maps the local portion of a vector into a
1420: vector.
1422: Collective on v, not collective on w.
1424: Input parameter:
1425: . v - The vector for which the local vector is desired.
1427: Output parameter:
1428: . w - Upon exit this contains the local vector.
1430: Level: beginner
1432: Notes:
1433: This function is similar to VecGetArray() which maps the local
1434: portion into a raw pointer. VecGetLocalVector() is usually about as
1435: efficient as VecGetArray() but in certain circumstances
1436: VecGetLocalVector() can be much more efficient than VecGetArray().
1437: This is because the construction of a contiguous array representing
1438: the vector data required by VecGetArray() can be an expensive
1439: operation for certain vector types. For example, for GPU vectors
1440: VecGetArray() requires that the data between device and host is
1441: synchronized.
1443: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1444: @*/
1445: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1446: {
1448: PetscScalar *a;
1453: VecCheckSameLocalSize(v,1,w,2);
1454: if (v->ops->getlocalvector) {
1455: (*v->ops->getlocalvector)(v,w);
1456: } else {
1457: VecGetArray(v,&a);
1458: VecPlaceArray(w,a);
1459: }
1460: return(0);
1461: }
1463: /*@
1464: VecRestoreLocalVector - Unmaps the local portion of a vector
1465: previously mapped into a vector using VecGetLocalVector().
1467: Logically collective.
1469: Input parameter:
1470: + v - The local portion of this vector was previously mapped into w using VecGetLocalVector().
1471: - w - The vector into which the local portion of v was mapped.
1473: Level: beginner
1475: .seealso: VecGetLocalVector(), VecGetLocalVectorRead(), VecRestoreLocalVectorRead(), LocalVectorRead(), VecGetArrayRead(), VecGetArray()
1476: @*/
1477: PetscErrorCode VecRestoreLocalVector(Vec v,Vec w)
1478: {
1480: PetscScalar *a;
1485: if (v->ops->restorelocalvector) {
1486: (*v->ops->restorelocalvector)(v,w);
1487: } else {
1488: VecGetArray(w,&a);
1489: VecRestoreArray(v,&a);
1490: VecResetArray(w);
1491: }
1492: return(0);
1493: }
1495: /*@C
1496: VecGetArray - Returns a pointer to a contiguous array that contains this
1497: processor's portion of the vector data. For the standard PETSc
1498: vectors, VecGetArray() returns a pointer to the local data array and
1499: does not use any copies. If the underlying vector data is not stored
1500: in a contiguous array this routine will copy the data to a contiguous
1501: array and return a pointer to that. You MUST call VecRestoreArray()
1502: when you no longer need access to the array.
1504: Logically Collective on Vec
1506: Input Parameter:
1507: . x - the vector
1509: Output Parameter:
1510: . a - location to put pointer to the array
1512: Fortran Note:
1513: This routine is used differently from Fortran 77
1514: $ Vec x
1515: $ PetscScalar x_array(1)
1516: $ PetscOffset i_x
1517: $ PetscErrorCode ierr
1518: $ call VecGetArray(x,x_array,i_x,ierr)
1519: $
1520: $ Access first local entry in vector with
1521: $ value = x_array(i_x + 1)
1522: $
1523: $ ...... other code
1524: $ call VecRestoreArray(x,x_array,i_x,ierr)
1525: For Fortran 90 see VecGetArrayF90()
1527: See the Fortran chapter of the users manual and
1528: petsc/src/snes/examples/tutorials/ex5f.F for details.
1530: Level: beginner
1532: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1533: VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite(), VecRestoreArrayWrite()
1534: @*/
1535: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1536: {
1538: #if defined(PETSC_HAVE_VIENNACL)
1539: PetscBool is_viennacltype = PETSC_FALSE;
1540: #endif
1544: VecSetErrorIfLocked(x,1);
1545: if (x->petscnative) {
1546: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1547: if (x->valid_GPU_array == PETSC_OFFLOAD_GPU) {
1548: #if defined(PETSC_HAVE_VIENNACL)
1549: PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1550: if (is_viennacltype) {
1551: VecViennaCLCopyFromGPU(x);
1552: } else
1553: #endif
1554: {
1555: #if defined(PETSC_HAVE_CUDA)
1556: VecCUDACopyFromGPU(x);
1557: #endif
1558: }
1559: } else if (x->valid_GPU_array == PETSC_OFFLOAD_UNALLOCATED) {
1560: #if defined(PETSC_HAVE_VIENNACL)
1561: PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1562: if (is_viennacltype) {
1563: VecViennaCLAllocateCheckHost(x);
1564: } else
1565: #endif
1566: {
1567: #if defined(PETSC_HAVE_CUDA)
1568: VecCUDAAllocateCheckHost(x);
1569: #endif
1570: }
1571: }
1572: #endif
1573: *a = *((PetscScalar**)x->data);
1574: } else {
1575: if (x->ops->getarray) {
1576: (*x->ops->getarray)(x,a);
1577: } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array for vector type \"%s\"",((PetscObject)x)->type_name);
1578: }
1579: return(0);
1580: }
1582: /*@C
1583: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contains this
1584: processor's portion of the vector data. The values in this array are NOT valid, the routine calling this
1585: routine is responsible for putting values into the array; any values it does not set will be invalid
1587: Logically Collective on Vec
1589: Input Parameter:
1590: . x - the vector
1592: Output Parameter:
1593: . a - location to put pointer to the array
1595: Level: intermediate
1597: This is for vectors associate with GPUs, the vector is not copied up before giving access. If you need correct
1598: values in the array use VecGetArray()
1600: Concepts: vector^accessing local values
1602: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1603: VecGetArrayPair(), VecRestoreArrayPair(), VecGetArray(), VecRestoreArrayWrite()
1604: @*/
1605: PetscErrorCode VecGetArrayWrite(Vec x,PetscScalar **a)
1606: {
1611: VecSetErrorIfLocked(x,1);
1612: if (!x->ops->getarraywrite) {
1613: VecGetArray(x,a);
1614: } else {
1615: (*x->ops->getarraywrite)(x,a);
1616: }
1617: return(0);
1618: }
1620: /*@C
1621: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
1623: Not Collective
1625: Input Parameters:
1626: . x - the vector
1628: Output Parameter:
1629: . a - the array
1631: Level: beginner
1633: Notes:
1634: The array must be returned using a matching call to VecRestoreArrayRead().
1636: Unlike VecGetArray(), this routine is not collective and preserves cached information like vector norms.
1638: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
1639: implementations may require a copy, but must such implementations should cache the contiguous representation so that
1640: only one copy is performed when this routine is called multiple times in sequence.
1642: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1643: @*/
1644: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1645: {
1647: #if defined(PETSC_HAVE_VIENNACL)
1648: PetscBool is_viennacltype = PETSC_FALSE;
1649: #endif
1653: if (x->petscnative) {
1654: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1655: if (x->valid_GPU_array == PETSC_OFFLOAD_GPU) {
1656: #if defined(PETSC_HAVE_VIENNACL)
1657: PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1658: if (is_viennacltype) {
1659: VecViennaCLCopyFromGPU(x);
1660: } else
1661: #endif
1662: {
1663: #if defined(PETSC_HAVE_CUDA)
1664: VecCUDACopyFromGPU(x);
1665: #endif
1666: }
1667: }
1668: #endif
1669: *a = *((PetscScalar **)x->data);
1670: } else if (x->ops->getarrayread) {
1671: (*x->ops->getarrayread)(x,a);
1672: } else {
1673: (*x->ops->getarray)(x,(PetscScalar**)a);
1674: }
1675: return(0);
1676: }
1678: /*@C
1679: VecGetArrays - Returns a pointer to the arrays in a set of vectors
1680: that were created by a call to VecDuplicateVecs(). You MUST call
1681: VecRestoreArrays() when you no longer need access to the array.
1683: Logically Collective on Vec
1685: Input Parameter:
1686: + x - the vectors
1687: - n - the number of vectors
1689: Output Parameter:
1690: . a - location to put pointer to the array
1692: Fortran Note:
1693: This routine is not supported in Fortran.
1695: Level: intermediate
1697: .seealso: VecGetArray(), VecRestoreArrays()
1698: @*/
1699: PetscErrorCode VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1700: {
1702: PetscInt i;
1703: PetscScalar **q;
1709: if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1710: PetscMalloc1(n,&q);
1711: for (i=0; i<n; ++i) {
1712: VecGetArray(x[i],&q[i]);
1713: }
1714: *a = q;
1715: return(0);
1716: }
1718: /*@C
1719: VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1720: has been called.
1722: Logically Collective on Vec
1724: Input Parameters:
1725: + x - the vector
1726: . n - the number of vectors
1727: - a - location of pointer to arrays obtained from VecGetArrays()
1729: Notes:
1730: For regular PETSc vectors this routine does not involve any copies. For
1731: any special vectors that do not store local vector data in a contiguous
1732: array, this routine will copy the data back into the underlying
1733: vector data structure from the arrays obtained with VecGetArrays().
1735: Fortran Note:
1736: This routine is not supported in Fortran.
1738: Level: intermediate
1740: .seealso: VecGetArrays(), VecRestoreArray()
1741: @*/
1742: PetscErrorCode VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1743: {
1745: PetscInt i;
1746: PetscScalar **q = *a;
1753: for (i=0; i<n; ++i) {
1754: VecRestoreArray(x[i],&q[i]);
1755: }
1756: PetscFree(q);
1757: return(0);
1758: }
1760: /*@C
1761: VecRestoreArray - Restores a vector after VecGetArray() has been called.
1763: Logically Collective on Vec
1765: Input Parameters:
1766: + x - the vector
1767: - a - location of pointer to array obtained from VecGetArray()
1769: Level: beginner
1771: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1772: VecGetArrayPair(), VecRestoreArrayPair()
1773: @*/
1774: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1775: {
1780: if (x->petscnative) {
1781: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1782: x->valid_GPU_array = PETSC_OFFLOAD_CPU;
1783: #endif
1784: } else {
1785: (*x->ops->restorearray)(x,a);
1786: }
1787: if (a) *a = NULL;
1788: PetscObjectStateIncrease((PetscObject)x);
1789: return(0);
1790: }
1792: /*@C
1793: VecRestoreArrayWrite - Restores a vector after VecGetArrayWrite() has been called.
1795: Logically Collective on Vec
1797: Input Parameters:
1798: + x - the vector
1799: - a - location of pointer to array obtained from VecGetArray()
1801: Level: beginner
1803: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1804: VecGetArrayPair(), VecRestoreArrayPair(), VecGetArrayWrite()
1805: @*/
1806: PetscErrorCode VecRestoreArrayWrite(Vec x,PetscScalar **a)
1807: {
1812: if (x->petscnative) {
1813: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1814: x->valid_GPU_array = PETSC_OFFLOAD_CPU;
1815: #endif
1816: } else {
1817: if (x->ops->restorearraywrite) {
1818: (*x->ops->restorearraywrite)(x,a);
1819: } else {
1820: (*x->ops->restorearray)(x,a);
1821: }
1822: }
1823: if (a) *a = NULL;
1824: PetscObjectStateIncrease((PetscObject)x);
1825: return(0);
1826: }
1828: /*@C
1829: VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()
1831: Not Collective
1833: Input Parameters:
1834: + vec - the vector
1835: - array - the array
1837: Level: beginner
1839: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1840: @*/
1841: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1842: {
1847: if (x->petscnative) {
1848: /* nothing */
1849: } else if (x->ops->restorearrayread) {
1850: (*x->ops->restorearrayread)(x,a);
1851: } else {
1852: (*x->ops->restorearray)(x,(PetscScalar**)a);
1853: }
1854: if (a) *a = NULL;
1855: return(0);
1856: }
1858: /*@
1859: VecPlaceArray - Allows one to replace the array in a vector with an
1860: array provided by the user. This is useful to avoid copying an array
1861: into a vector.
1863: Not Collective
1865: Input Parameters:
1866: + vec - the vector
1867: - array - the array
1869: Notes:
1870: You can return to the original array with a call to VecResetArray()
1872: Level: developer
1874: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecResetArray()
1876: @*/
1877: PetscErrorCode VecPlaceArray(Vec vec,const PetscScalar array[])
1878: {
1885: if (vec->ops->placearray) {
1886: (*vec->ops->placearray)(vec,array);
1887: } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
1888: PetscObjectStateIncrease((PetscObject)vec);
1889: return(0);
1890: }
1892: /*@C
1893: VecReplaceArray - Allows one to replace the array in a vector with an
1894: array provided by the user. This is useful to avoid copying an array
1895: into a vector.
1897: Not Collective
1899: Input Parameters:
1900: + vec - the vector
1901: - array - the array
1903: Notes:
1904: This permanently replaces the array and frees the memory associated
1905: with the old array.
1907: The memory passed in MUST be obtained with PetscMalloc() and CANNOT be
1908: freed by the user. It will be freed when the vector is destroy.
1910: Not supported from Fortran
1912: Level: developer
1914: .seealso: VecGetArray(), VecRestoreArray(), VecPlaceArray(), VecResetArray()
1916: @*/
1917: PetscErrorCode VecReplaceArray(Vec vec,const PetscScalar array[])
1918: {
1924: if (vec->ops->replacearray) {
1925: (*vec->ops->replacearray)(vec,array);
1926: } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
1927: PetscObjectStateIncrease((PetscObject)vec);
1928: return(0);
1929: }
1931: /*MC
1932: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
1933: and makes them accessible via a Fortran90 pointer.
1935: Synopsis:
1936: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
1938: Collective on Vec
1940: Input Parameters:
1941: + x - a vector to mimic
1942: - n - the number of vectors to obtain
1944: Output Parameters:
1945: + y - Fortran90 pointer to the array of vectors
1946: - ierr - error code
1948: Example of Usage:
1949: .vb
1950: #include <petsc/finclude/petscvec.h>
1951: use petscvec
1953: Vec x
1954: Vec, pointer :: y(:)
1955: ....
1956: call VecDuplicateVecsF90(x,2,y,ierr)
1957: call VecSet(y(2),alpha,ierr)
1958: call VecSet(y(2),alpha,ierr)
1959: ....
1960: call VecDestroyVecsF90(2,y,ierr)
1961: .ve
1963: Notes:
1964: Not yet supported for all F90 compilers
1966: Use VecDestroyVecsF90() to free the space.
1968: Level: beginner
1970: .seealso: VecDestroyVecsF90(), VecDuplicateVecs()
1972: M*/
1974: /*MC
1975: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
1976: VecGetArrayF90().
1978: Synopsis:
1979: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
1981: Logically Collective on Vec
1983: Input Parameters:
1984: + x - vector
1985: - xx_v - the Fortran90 pointer to the array
1987: Output Parameter:
1988: . ierr - error code
1990: Example of Usage:
1991: .vb
1992: #include <petsc/finclude/petscvec.h>
1993: use petscvec
1995: PetscScalar, pointer :: xx_v(:)
1996: ....
1997: call VecGetArrayF90(x,xx_v,ierr)
1998: xx_v(3) = a
1999: call VecRestoreArrayF90(x,xx_v,ierr)
2000: .ve
2002: Level: beginner
2004: .seealso: VecGetArrayF90(), VecGetArray(), VecRestoreArray(), UsingFortran, VecRestoreArrayReadF90()
2006: M*/
2008: /*MC
2009: VecDestroyVecsF90 - Frees a block of vectors obtained with VecDuplicateVecsF90().
2011: Synopsis:
2012: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
2014: Collective on Vec
2016: Input Parameters:
2017: + n - the number of vectors previously obtained
2018: - x - pointer to array of vector pointers
2020: Output Parameter:
2021: . ierr - error code
2023: Notes:
2024: Not yet supported for all F90 compilers
2026: Level: beginner
2028: .seealso: VecDestroyVecs(), VecDuplicateVecsF90()
2030: M*/
2032: /*MC
2033: VecGetArrayF90 - Accesses a vector array from Fortran90. For default PETSc
2034: vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
2035: this routine is implementation dependent. You MUST call VecRestoreArrayF90()
2036: when you no longer need access to the array.
2038: Synopsis:
2039: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2041: Logically Collective on Vec
2043: Input Parameter:
2044: . x - vector
2046: Output Parameters:
2047: + xx_v - the Fortran90 pointer to the array
2048: - ierr - error code
2050: Example of Usage:
2051: .vb
2052: #include <petsc/finclude/petscvec.h>
2053: use petscvec
2055: PetscScalar, pointer :: xx_v(:)
2056: ....
2057: call VecGetArrayF90(x,xx_v,ierr)
2058: xx_v(3) = a
2059: call VecRestoreArrayF90(x,xx_v,ierr)
2060: .ve
2062: If you ONLY intend to read entries from the array and not change any entries you should use VecGetArrayReadF90().
2064: Level: beginner
2066: .seealso: VecRestoreArrayF90(), VecGetArray(), VecRestoreArray(), VecGetArrayReadF90(), UsingFortran
2068: M*/
2070: /*MC
2071: VecGetArrayReadF90 - Accesses a read only array from Fortran90. For default PETSc
2072: vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
2073: this routine is implementation dependent. You MUST call VecRestoreArrayReadF90()
2074: when you no longer need access to the array.
2076: Synopsis:
2077: VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2079: Logically Collective on Vec
2081: Input Parameter:
2082: . x - vector
2084: Output Parameters:
2085: + xx_v - the Fortran90 pointer to the array
2086: - ierr - error code
2088: Example of Usage:
2089: .vb
2090: #include <petsc/finclude/petscvec.h>
2091: use petscvec
2093: PetscScalar, pointer :: xx_v(:)
2094: ....
2095: call VecGetArrayReadF90(x,xx_v,ierr)
2096: a = xx_v(3)
2097: call VecRestoreArrayReadF90(x,xx_v,ierr)
2098: .ve
2100: If you intend to write entries into the array you must use VecGetArrayF90().
2102: Level: beginner
2104: .seealso: VecRestoreArrayReadF90(), VecGetArray(), VecRestoreArray(), VecGetArrayRead(), VecRestoreArrayRead(), VecGetArrayF90(), UsingFortran
2106: M*/
2108: /*MC
2109: VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2110: VecGetArrayReadF90().
2112: Synopsis:
2113: VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2115: Logically Collective on Vec
2117: Input Parameters:
2118: + x - vector
2119: - xx_v - the Fortran90 pointer to the array
2121: Output Parameter:
2122: . ierr - error code
2124: Example of Usage:
2125: .vb
2126: #include <petsc/finclude/petscvec.h>
2127: use petscvec
2129: PetscScalar, pointer :: xx_v(:)
2130: ....
2131: call VecGetArrayReadF90(x,xx_v,ierr)
2132: a = xx_v(3)
2133: call VecRestoreArrayReadF90(x,xx_v,ierr)
2134: .ve
2136: Level: beginner
2138: .seealso: VecGetArrayReadF90(), VecGetArray(), VecRestoreArray(), VecGetArrayRead(), VecRestoreArrayRead(),UsingFortran, VecRestoreArrayF90()
2140: M*/
2142: /*@C
2143: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2144: processor's portion of the vector data. You MUST call VecRestoreArray2d()
2145: when you no longer need access to the array.
2147: Logically Collective
2149: Input Parameter:
2150: + x - the vector
2151: . m - first dimension of two dimensional array
2152: . n - second dimension of two dimensional array
2153: . mstart - first index you will use in first coordinate direction (often 0)
2154: - nstart - first index in the second coordinate direction (often 0)
2156: Output Parameter:
2157: . a - location to put pointer to the array
2159: Level: developer
2161: Notes:
2162: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2163: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2164: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2165: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
2167: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2169: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2170: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2171: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2172: @*/
2173: PetscErrorCode VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2174: {
2176: PetscInt i,N;
2177: PetscScalar *aa;
2183: VecGetLocalSize(x,&N);
2184: 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);
2185: VecGetArray(x,&aa);
2187: PetscMalloc1(m,a);
2188: for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2189: *a -= mstart;
2190: return(0);
2191: }
2193: /*@C
2194: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2195: processor's portion of the vector data. You MUST call VecRestoreArray2dWrite()
2196: when you no longer need access to the array.
2198: Logically Collective
2200: Input Parameter:
2201: + x - the vector
2202: . m - first dimension of two dimensional array
2203: . n - second dimension of two dimensional array
2204: . mstart - first index you will use in first coordinate direction (often 0)
2205: - nstart - first index in the second coordinate direction (often 0)
2207: Output Parameter:
2208: . a - location to put pointer to the array
2210: Level: developer
2212: Notes:
2213: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2214: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2215: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2216: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
2218: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2220: Concepts: vector^accessing local values as 2d array
2222: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2223: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2224: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2225: @*/
2226: PetscErrorCode VecGetArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2227: {
2229: PetscInt i,N;
2230: PetscScalar *aa;
2236: VecGetLocalSize(x,&N);
2237: if (m*n != N) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 2d array dimensions %D by %D",N,m,n);
2238: VecGetArrayWrite(x,&aa);
2240: PetscMalloc1(m,a);
2241: for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2242: *a -= mstart;
2243: return(0);
2244: }
2246: /*@C
2247: VecRestoreArray2d - Restores a vector after VecGetArray2d() has been called.
2249: Logically Collective
2251: Input Parameters:
2252: + x - the vector
2253: . m - first dimension of two dimensional array
2254: . n - second dimension of the two dimensional array
2255: . mstart - first index you will use in first coordinate direction (often 0)
2256: . nstart - first index in the second coordinate direction (often 0)
2257: - a - location of pointer to array obtained from VecGetArray2d()
2259: Level: developer
2261: Notes:
2262: For regular PETSc vectors this routine does not involve any copies. For
2263: any special vectors that do not store local vector data in a contiguous
2264: array, this routine will copy the data back into the underlying
2265: vector data structure from the array obtained with VecGetArray().
2267: This routine actually zeros out the a pointer.
2269: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2270: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2271: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2272: @*/
2273: PetscErrorCode VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2274: {
2276: void *dummy;
2282: dummy = (void*)(*a + mstart);
2283: PetscFree(dummy);
2284: VecRestoreArray(x,NULL);
2285: return(0);
2286: }
2288: /*@C
2289: VecRestoreArray2dWrite - Restores a vector after VecGetArray2dWrite() has been called.
2291: Logically Collective
2293: Input Parameters:
2294: + x - the vector
2295: . m - first dimension of two dimensional array
2296: . n - second dimension of the two dimensional array
2297: . mstart - first index you will use in first coordinate direction (often 0)
2298: . nstart - first index in the second coordinate direction (often 0)
2299: - a - location of pointer to array obtained from VecGetArray2d()
2301: Level: developer
2303: Notes:
2304: For regular PETSc vectors this routine does not involve any copies. For
2305: any special vectors that do not store local vector data in a contiguous
2306: array, this routine will copy the data back into the underlying
2307: vector data structure from the array obtained with VecGetArray().
2309: This routine actually zeros out the a pointer.
2311: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2312: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2313: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2314: @*/
2315: PetscErrorCode VecRestoreArray2dWrite(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2316: {
2318: void *dummy;
2324: dummy = (void*)(*a + mstart);
2325: PetscFree(dummy);
2326: VecRestoreArrayWrite(x,NULL);
2327: return(0);
2328: }
2330: /*@C
2331: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2332: processor's portion of the vector data. You MUST call VecRestoreArray1d()
2333: when you no longer need access to the array.
2335: Logically Collective
2337: Input Parameter:
2338: + x - the vector
2339: . m - first dimension of two dimensional array
2340: - mstart - first index you will use in first coordinate direction (often 0)
2342: Output Parameter:
2343: . a - location to put pointer to the array
2345: Level: developer
2347: Notes:
2348: For a vector obtained from DMCreateLocalVector() mstart are likely
2349: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2350: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
2352: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2354: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2355: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2356: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2357: @*/
2358: PetscErrorCode VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2359: {
2361: PetscInt N;
2367: VecGetLocalSize(x,&N);
2368: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2369: VecGetArray(x,a);
2370: *a -= mstart;
2371: return(0);
2372: }
2374: /*@C
2375: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
2376: processor's portion of the vector data. You MUST call VecRestoreArray1dWrite()
2377: when you no longer need access to the array.
2379: Logically Collective
2381: Input Parameter:
2382: + x - the vector
2383: . m - first dimension of two dimensional array
2384: - mstart - first index you will use in first coordinate direction (often 0)
2386: Output Parameter:
2387: . a - location to put pointer to the array
2389: Level: developer
2391: Notes:
2392: For a vector obtained from DMCreateLocalVector() mstart are likely
2393: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2394: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
2396: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2398: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2399: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2400: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2401: @*/
2402: PetscErrorCode VecGetArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2403: {
2405: PetscInt N;
2411: VecGetLocalSize(x,&N);
2412: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2413: VecGetArrayWrite(x,a);
2414: *a -= mstart;
2415: return(0);
2416: }
2418: /*@C
2419: VecRestoreArray1d - Restores a vector after VecGetArray1d() has been called.
2421: Logically Collective
2423: Input Parameters:
2424: + x - the vector
2425: . m - first dimension of two dimensional array
2426: . mstart - first index you will use in first coordinate direction (often 0)
2427: - a - location of pointer to array obtained from VecGetArray21()
2429: Level: developer
2431: Notes:
2432: For regular PETSc vectors this routine does not involve any copies. For
2433: any special vectors that do not store local vector data in a contiguous
2434: array, this routine will copy the data back into the underlying
2435: vector data structure from the array obtained with VecGetArray1d().
2437: This routine actually zeros out the a pointer.
2439: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2440: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2441: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2442: @*/
2443: PetscErrorCode VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2444: {
2450: VecRestoreArray(x,NULL);
2451: return(0);
2452: }
2454: /*@C
2455: VecRestoreArray1dWrite - Restores a vector after VecGetArray1dWrite() has been called.
2457: Logically Collective
2459: Input Parameters:
2460: + x - the vector
2461: . m - first dimension of two dimensional array
2462: . mstart - first index you will use in first coordinate direction (often 0)
2463: - a - location of pointer to array obtained from VecGetArray21()
2465: Level: developer
2467: Notes:
2468: For regular PETSc vectors this routine does not involve any copies. For
2469: any special vectors that do not store local vector data in a contiguous
2470: array, this routine will copy the data back into the underlying
2471: vector data structure from the array obtained with VecGetArray1d().
2473: This routine actually zeros out the a pointer.
2475: Concepts: vector^accessing local values as 1d array
2477: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2478: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2479: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2480: @*/
2481: PetscErrorCode VecRestoreArray1dWrite(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2482: {
2488: VecRestoreArrayWrite(x,NULL);
2489: return(0);
2490: }
2492: /*@C
2493: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
2494: processor's portion of the vector data. You MUST call VecRestoreArray3d()
2495: when you no longer need access to the array.
2497: Logically Collective
2499: Input Parameter:
2500: + x - the vector
2501: . m - first dimension of three dimensional array
2502: . n - second dimension of three dimensional array
2503: . p - third dimension of three dimensional array
2504: . mstart - first index you will use in first coordinate direction (often 0)
2505: . nstart - first index in the second coordinate direction (often 0)
2506: - pstart - first index in the third coordinate direction (often 0)
2508: Output Parameter:
2509: . a - location to put pointer to the array
2511: Level: developer
2513: Notes:
2514: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
2515: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2516: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2517: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
2519: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2521: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2522: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2523: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2524: @*/
2525: PetscErrorCode VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2526: {
2528: PetscInt i,N,j;
2529: PetscScalar *aa,**b;
2535: VecGetLocalSize(x,&N);
2536: 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);
2537: VecGetArray(x,&aa);
2539: PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
2540: b = (PetscScalar**)((*a) + m);
2541: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2542: for (i=0; i<m; i++)
2543: for (j=0; j<n; j++)
2544: b[i*n+j] = aa + i*n*p + j*p - pstart;
2546: *a -= mstart;
2547: return(0);
2548: }
2550: /*@C
2551: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
2552: processor's portion of the vector data. You MUST call VecRestoreArray3dWrite()
2553: when you no longer need access to the array.
2555: Logically Collective
2557: Input Parameter:
2558: + x - the vector
2559: . m - first dimension of three dimensional array
2560: . n - second dimension of three dimensional array
2561: . p - third dimension of three dimensional array
2562: . mstart - first index you will use in first coordinate direction (often 0)
2563: . nstart - first index in the second coordinate direction (often 0)
2564: - pstart - first index in the third coordinate direction (often 0)
2566: Output Parameter:
2567: . a - location to put pointer to the array
2569: Level: developer
2571: Notes:
2572: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
2573: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2574: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2575: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
2577: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2579: Concepts: vector^accessing local values as 3d array
2581: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2582: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2583: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2584: @*/
2585: PetscErrorCode VecGetArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2586: {
2588: PetscInt i,N,j;
2589: PetscScalar *aa,**b;
2595: VecGetLocalSize(x,&N);
2596: 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);
2597: VecGetArrayWrite(x,&aa);
2599: PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
2600: b = (PetscScalar**)((*a) + m);
2601: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2602: for (i=0; i<m; i++)
2603: for (j=0; j<n; j++)
2604: b[i*n+j] = aa + i*n*p + j*p - pstart;
2606: *a -= mstart;
2607: return(0);
2608: }
2610: /*@C
2611: VecRestoreArray3d - Restores a vector after VecGetArray3d() has been called.
2613: Logically Collective
2615: Input Parameters:
2616: + x - the vector
2617: . m - first dimension of three dimensional array
2618: . n - second dimension of the three dimensional array
2619: . p - third dimension of the three dimensional array
2620: . mstart - first index you will use in first coordinate direction (often 0)
2621: . nstart - first index in the second coordinate direction (often 0)
2622: . pstart - first index in the third coordinate direction (often 0)
2623: - a - location of pointer to array obtained from VecGetArray3d()
2625: Level: developer
2627: Notes:
2628: For regular PETSc vectors this routine does not involve any copies. For
2629: any special vectors that do not store local vector data in a contiguous
2630: array, this routine will copy the data back into the underlying
2631: vector data structure from the array obtained with VecGetArray().
2633: This routine actually zeros out the a pointer.
2635: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2636: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2637: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2638: @*/
2639: PetscErrorCode VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2640: {
2642: void *dummy;
2648: dummy = (void*)(*a + mstart);
2649: PetscFree(dummy);
2650: VecRestoreArray(x,NULL);
2651: return(0);
2652: }
2654: /*@C
2655: VecRestoreArray3dWrite - Restores a vector after VecGetArray3dWrite() has been called.
2657: Logically Collective
2659: Input Parameters:
2660: + x - the vector
2661: . m - first dimension of three dimensional array
2662: . n - second dimension of the three dimensional array
2663: . p - third dimension of the three dimensional array
2664: . mstart - first index you will use in first coordinate direction (often 0)
2665: . nstart - first index in the second coordinate direction (often 0)
2666: . pstart - first index in the third coordinate direction (often 0)
2667: - a - location of pointer to array obtained from VecGetArray3d()
2669: Level: developer
2671: Notes:
2672: For regular PETSc vectors this routine does not involve any copies. For
2673: any special vectors that do not store local vector data in a contiguous
2674: array, this routine will copy the data back into the underlying
2675: vector data structure from the array obtained with VecGetArray().
2677: This routine actually zeros out the a pointer.
2679: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2680: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2681: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2682: @*/
2683: PetscErrorCode VecRestoreArray3dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2684: {
2686: void *dummy;
2692: dummy = (void*)(*a + mstart);
2693: PetscFree(dummy);
2694: VecRestoreArrayWrite(x,NULL);
2695: return(0);
2696: }
2698: /*@C
2699: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
2700: processor's portion of the vector data. You MUST call VecRestoreArray4d()
2701: when you no longer need access to the array.
2703: Logically Collective
2705: Input Parameter:
2706: + x - the vector
2707: . m - first dimension of four dimensional array
2708: . n - second dimension of four dimensional array
2709: . p - third dimension of four dimensional array
2710: . q - fourth dimension of four dimensional array
2711: . mstart - first index you will use in first coordinate direction (often 0)
2712: . nstart - first index in the second coordinate direction (often 0)
2713: . pstart - first index in the third coordinate direction (often 0)
2714: - qstart - first index in the fourth coordinate direction (often 0)
2716: Output Parameter:
2717: . a - location to put pointer to the array
2719: Level: beginner
2721: Notes:
2722: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
2723: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2724: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2725: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
2727: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2729: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2730: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2731: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2732: @*/
2733: PetscErrorCode VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2734: {
2736: PetscInt i,N,j,k;
2737: PetscScalar *aa,***b,**c;
2743: VecGetLocalSize(x,&N);
2744: 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);
2745: VecGetArray(x,&aa);
2747: PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
2748: b = (PetscScalar***)((*a) + m);
2749: c = (PetscScalar**)(b + m*n);
2750: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2751: for (i=0; i<m; i++)
2752: for (j=0; j<n; j++)
2753: b[i*n+j] = c + i*n*p + j*p - pstart;
2754: for (i=0; i<m; i++)
2755: for (j=0; j<n; j++)
2756: for (k=0; k<p; k++)
2757: c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
2758: *a -= mstart;
2759: return(0);
2760: }
2762: /*@C
2763: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
2764: processor's portion of the vector data. You MUST call VecRestoreArray4dWrite()
2765: when you no longer need access to the array.
2767: Logically Collective
2769: Input Parameter:
2770: + x - the vector
2771: . m - first dimension of four dimensional array
2772: . n - second dimension of four dimensional array
2773: . p - third dimension of four dimensional array
2774: . q - fourth dimension of four 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)
2777: . pstart - first index in the third coordinate direction (often 0)
2778: - qstart - first index in the fourth coordinate direction (often 0)
2780: Output Parameter:
2781: . a - location to put pointer to the array
2783: Level: beginner
2785: Notes:
2786: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
2787: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2788: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2789: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
2791: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2793: Concepts: vector^accessing local values as 3d array
2795: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2796: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2797: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2798: @*/
2799: PetscErrorCode VecGetArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2800: {
2802: PetscInt i,N,j,k;
2803: PetscScalar *aa,***b,**c;
2809: VecGetLocalSize(x,&N);
2810: 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);
2811: VecGetArrayWrite(x,&aa);
2813: PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
2814: b = (PetscScalar***)((*a) + m);
2815: c = (PetscScalar**)(b + m*n);
2816: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2817: for (i=0; i<m; i++)
2818: for (j=0; j<n; j++)
2819: b[i*n+j] = c + i*n*p + j*p - pstart;
2820: for (i=0; i<m; i++)
2821: for (j=0; j<n; j++)
2822: for (k=0; k<p; k++)
2823: c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
2824: *a -= mstart;
2825: return(0);
2826: }
2828: /*@C
2829: VecRestoreArray4d - Restores a vector after VecGetArray3d() has been called.
2831: Logically Collective
2833: Input Parameters:
2834: + x - the vector
2835: . m - first dimension of four dimensional array
2836: . n - second dimension of the four dimensional array
2837: . p - third dimension of the four dimensional array
2838: . q - fourth dimension of the four dimensional array
2839: . mstart - first index you will use in first coordinate direction (often 0)
2840: . nstart - first index in the second coordinate direction (often 0)
2841: . pstart - first index in the third coordinate direction (often 0)
2842: . qstart - first index in the fourth coordinate direction (often 0)
2843: - a - location of pointer to array obtained from VecGetArray4d()
2845: Level: beginner
2847: Notes:
2848: For regular PETSc vectors this routine does not involve any copies. For
2849: any special vectors that do not store local vector data in a contiguous
2850: array, this routine will copy the data back into the underlying
2851: vector data structure from the array obtained with VecGetArray().
2853: This routine actually zeros out the a pointer.
2855: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2856: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2857: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2858: @*/
2859: PetscErrorCode VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2860: {
2862: void *dummy;
2868: dummy = (void*)(*a + mstart);
2869: PetscFree(dummy);
2870: VecRestoreArray(x,NULL);
2871: return(0);
2872: }
2874: /*@C
2875: VecRestoreArray4dWrite - Restores a vector after VecGetArray3dWrite() has been called.
2877: Logically Collective
2879: Input Parameters:
2880: + x - the vector
2881: . m - first dimension of four dimensional array
2882: . n - second dimension of the four dimensional array
2883: . p - third dimension of the four dimensional array
2884: . q - fourth dimension of the four dimensional array
2885: . mstart - first index you will use in first coordinate direction (often 0)
2886: . nstart - first index in the second coordinate direction (often 0)
2887: . pstart - first index in the third coordinate direction (often 0)
2888: . qstart - first index in the fourth coordinate direction (often 0)
2889: - a - location of pointer to array obtained from VecGetArray4d()
2891: Level: beginner
2893: Notes:
2894: For regular PETSc vectors this routine does not involve any copies. For
2895: any special vectors that do not store local vector data in a contiguous
2896: array, this routine will copy the data back into the underlying
2897: vector data structure from the array obtained with VecGetArray().
2899: This routine actually zeros out the a pointer.
2901: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2902: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2903: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2904: @*/
2905: PetscErrorCode VecRestoreArray4dWrite(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2906: {
2908: void *dummy;
2914: dummy = (void*)(*a + mstart);
2915: PetscFree(dummy);
2916: VecRestoreArrayWrite(x,NULL);
2917: return(0);
2918: }
2920: /*@C
2921: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
2922: processor's portion of the vector data. You MUST call VecRestoreArray2dRead()
2923: when you no longer need access to the array.
2925: Logically Collective
2927: Input Parameter:
2928: + x - the vector
2929: . m - first dimension of two dimensional array
2930: . n - second dimension of two dimensional array
2931: . mstart - first index you will use in first coordinate direction (often 0)
2932: - nstart - first index in the second coordinate direction (often 0)
2934: Output Parameter:
2935: . a - location to put pointer to the array
2937: Level: developer
2939: Notes:
2940: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2941: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2942: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2943: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
2945: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2947: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2948: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2949: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2950: @*/
2951: PetscErrorCode VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2952: {
2953: PetscErrorCode ierr;
2954: PetscInt i,N;
2955: const PetscScalar *aa;
2961: VecGetLocalSize(x,&N);
2962: 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);
2963: VecGetArrayRead(x,&aa);
2965: PetscMalloc1(m,a);
2966: for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
2967: *a -= mstart;
2968: return(0);
2969: }
2971: /*@C
2972: VecRestoreArray2dRead - Restores a vector after VecGetArray2dRead() has been called.
2974: Logically Collective
2976: Input Parameters:
2977: + x - the vector
2978: . m - first dimension of two dimensional array
2979: . n - second dimension of the two dimensional array
2980: . mstart - first index you will use in first coordinate direction (often 0)
2981: . nstart - first index in the second coordinate direction (often 0)
2982: - a - location of pointer to array obtained from VecGetArray2d()
2984: Level: developer
2986: Notes:
2987: For regular PETSc vectors this routine does not involve any copies. For
2988: any special vectors that do not store local vector data in a contiguous
2989: array, this routine will copy the data back into the underlying
2990: vector data structure from the array obtained with VecGetArray().
2992: This routine actually zeros out the a pointer.
2994: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2995: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2996: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2997: @*/
2998: PetscErrorCode VecRestoreArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2999: {
3001: void *dummy;
3007: dummy = (void*)(*a + mstart);
3008: PetscFree(dummy);
3009: VecRestoreArrayRead(x,NULL);
3010: return(0);
3011: }
3013: /*@C
3014: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3015: processor's portion of the vector data. You MUST call VecRestoreArray1dRead()
3016: when you no longer need access to the array.
3018: Logically Collective
3020: Input Parameter:
3021: + x - the vector
3022: . m - first dimension of two dimensional array
3023: - mstart - first index you will use in first coordinate direction (often 0)
3025: Output Parameter:
3026: . a - location to put pointer to the array
3028: Level: developer
3030: Notes:
3031: For a vector obtained from DMCreateLocalVector() mstart are likely
3032: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3033: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
3035: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3037: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3038: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3039: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3040: @*/
3041: PetscErrorCode VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3042: {
3044: PetscInt N;
3050: VecGetLocalSize(x,&N);
3051: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
3052: VecGetArrayRead(x,(const PetscScalar**)a);
3053: *a -= mstart;
3054: return(0);
3055: }
3057: /*@C
3058: VecRestoreArray1dRead - Restores a vector after VecGetArray1dRead() has been called.
3060: Logically Collective
3062: Input Parameters:
3063: + x - the vector
3064: . m - first dimension of two dimensional array
3065: . mstart - first index you will use in first coordinate direction (often 0)
3066: - a - location of pointer to array obtained from VecGetArray21()
3068: Level: developer
3070: Notes:
3071: For regular PETSc vectors this routine does not involve any copies. For
3072: any special vectors that do not store local vector data in a contiguous
3073: array, this routine will copy the data back into the underlying
3074: vector data structure from the array obtained with VecGetArray1dRead().
3076: This routine actually zeros out the a pointer.
3078: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3079: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3080: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
3081: @*/
3082: PetscErrorCode VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
3083: {
3089: VecRestoreArrayRead(x,NULL);
3090: return(0);
3091: }
3094: /*@C
3095: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3096: processor's portion of the vector data. You MUST call VecRestoreArray3dRead()
3097: when you no longer need access to the array.
3099: Logically Collective
3101: Input Parameter:
3102: + x - the vector
3103: . m - first dimension of three dimensional array
3104: . n - second dimension of three dimensional array
3105: . p - third dimension of three dimensional array
3106: . mstart - first index you will use in first coordinate direction (often 0)
3107: . nstart - first index in the second coordinate direction (often 0)
3108: - pstart - first index in the third coordinate direction (often 0)
3110: Output Parameter:
3111: . a - location to put pointer to the array
3113: Level: developer
3115: Notes:
3116: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3117: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3118: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3119: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3dRead().
3121: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3123: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3124: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3125: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3126: @*/
3127: PetscErrorCode VecGetArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3128: {
3129: PetscErrorCode ierr;
3130: PetscInt i,N,j;
3131: const PetscScalar *aa;
3132: PetscScalar **b;
3138: VecGetLocalSize(x,&N);
3139: 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);
3140: VecGetArrayRead(x,&aa);
3142: PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
3143: b = (PetscScalar**)((*a) + m);
3144: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3145: for (i=0; i<m; i++)
3146: for (j=0; j<n; j++)
3147: b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;
3149: *a -= mstart;
3150: return(0);
3151: }
3153: /*@C
3154: VecRestoreArray3dRead - Restores a vector after VecGetArray3dRead() has been called.
3156: Logically Collective
3158: Input Parameters:
3159: + x - the vector
3160: . m - first dimension of three dimensional array
3161: . n - second dimension of the three dimensional array
3162: . p - third dimension of the three dimensional array
3163: . mstart - first index you will use in first coordinate direction (often 0)
3164: . nstart - first index in the second coordinate direction (often 0)
3165: . pstart - first index in the third coordinate direction (often 0)
3166: - a - location of pointer to array obtained from VecGetArray3dRead()
3168: Level: developer
3170: Notes:
3171: For regular PETSc vectors this routine does not involve any copies. For
3172: any special vectors that do not store local vector data in a contiguous
3173: array, this routine will copy the data back into the underlying
3174: vector data structure from the array obtained with VecGetArray().
3176: This routine actually zeros out the a pointer.
3178: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3179: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3180: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3181: @*/
3182: PetscErrorCode VecRestoreArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
3183: {
3185: void *dummy;
3191: dummy = (void*)(*a + mstart);
3192: PetscFree(dummy);
3193: VecRestoreArrayRead(x,NULL);
3194: return(0);
3195: }
3197: /*@C
3198: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3199: processor's portion of the vector data. You MUST call VecRestoreArray4dRead()
3200: when you no longer need access to the array.
3202: Logically Collective
3204: Input Parameter:
3205: + x - the vector
3206: . m - first dimension of four dimensional array
3207: . n - second dimension of four dimensional array
3208: . p - third dimension of four dimensional array
3209: . q - fourth dimension of four dimensional array
3210: . mstart - first index you will use in first coordinate direction (often 0)
3211: . nstart - first index in the second coordinate direction (often 0)
3212: . pstart - first index in the third coordinate direction (often 0)
3213: - qstart - first index in the fourth coordinate direction (often 0)
3215: Output Parameter:
3216: . a - location to put pointer to the array
3218: Level: beginner
3220: Notes:
3221: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
3222: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
3223: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
3224: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
3226: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3228: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3229: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3230: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3231: @*/
3232: PetscErrorCode VecGetArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3233: {
3234: PetscErrorCode ierr;
3235: PetscInt i,N,j,k;
3236: const PetscScalar *aa;
3237: PetscScalar ***b,**c;
3243: VecGetLocalSize(x,&N);
3244: 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);
3245: VecGetArrayRead(x,&aa);
3247: PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
3248: b = (PetscScalar***)((*a) + m);
3249: c = (PetscScalar**)(b + m*n);
3250: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
3251: for (i=0; i<m; i++)
3252: for (j=0; j<n; j++)
3253: b[i*n+j] = c + i*n*p + j*p - pstart;
3254: for (i=0; i<m; i++)
3255: for (j=0; j<n; j++)
3256: for (k=0; k<p; k++)
3257: c[i*n*p+j*p+k] = (PetscScalar*) aa + i*n*p*q + j*p*q + k*q - qstart;
3258: *a -= mstart;
3259: return(0);
3260: }
3262: /*@C
3263: VecRestoreArray4dRead - Restores a vector after VecGetArray3d() has been called.
3265: Logically Collective
3267: Input Parameters:
3268: + x - the vector
3269: . m - first dimension of four dimensional array
3270: . n - second dimension of the four dimensional array
3271: . p - third dimension of the four dimensional array
3272: . q - fourth dimension of the four dimensional array
3273: . mstart - first index you will use in first coordinate direction (often 0)
3274: . nstart - first index in the second coordinate direction (often 0)
3275: . pstart - first index in the third coordinate direction (often 0)
3276: . qstart - first index in the fourth coordinate direction (often 0)
3277: - a - location of pointer to array obtained from VecGetArray4dRead()
3279: Level: beginner
3281: Notes:
3282: For regular PETSc vectors this routine does not involve any copies. For
3283: any special vectors that do not store local vector data in a contiguous
3284: array, this routine will copy the data back into the underlying
3285: vector data structure from the array obtained with VecGetArray().
3287: This routine actually zeros out the a pointer.
3289: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3290: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
3291: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3292: @*/
3293: PetscErrorCode VecRestoreArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
3294: {
3296: void *dummy;
3302: dummy = (void*)(*a + mstart);
3303: PetscFree(dummy);
3304: VecRestoreArrayRead(x,NULL);
3305: return(0);
3306: }
3308: #if defined(PETSC_USE_DEBUG)
3310: /*@
3311: VecLockGet - Gets the current lock status of a vector
3313: Logically Collective on Vec
3315: Input Parameter:
3316: . x - the vector
3318: Output Parameter:
3319: . state - greater than zero indicates the vector is locked for read; less then zero indicates the vector is
3320: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3322: Level: beginner
3324: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop()
3325: @*/
3326: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
3327: {
3330: *state = x->lock;
3331: return(0);
3332: }
3334: /*@
3335: VecLockReadPush - Pushes a read-only lock on a vector to prevent it from writing
3337: Logically Collective on Vec
3339: Input Parameter:
3340: . x - the vector
3342: Notes:
3343: If this is set then calls to VecGetArray() or VecSetValues() or any other routines that change the vectors values will fail.
3345: The call can be nested, i.e., called multiple times on the same vector, but each VecLockReadPush(x) has to have one matching
3346: VecLockReadPop(x), which removes the latest read-only lock.
3348: Level: beginner
3350: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPop(), VecLockGet()
3351: @*/
3352: PetscErrorCode VecLockReadPush(Vec x)
3353: {
3356: 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");
3357: x->lock++;
3358: return(0);
3359: }
3361: /*@
3362: VecLockReadPop - Pops a read-only lock from a vector
3364: Logically Collective on Vec
3366: Input Parameter:
3367: . x - the vector
3369: Level: beginner
3371: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockGet()
3372: @*/
3373: PetscErrorCode VecLockReadPop(Vec x)
3374: {
3377: x->lock--;
3378: if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector has been unlocked from read-only access too many times");
3379: return(0);
3380: }
3382: /*@C
3383: VecLockWriteSet_Private - Lock or unlock a vector for exclusive read/write access
3385: Logically Collective on Vec
3387: Input Parameter:
3388: + x - the vector
3389: - flg - PETSC_TRUE to lock the vector for writing; PETSC_FALSE to unlock it.
3391: Notes:
3392: The function is usefull in split-phase computations, which usually have a begin phase and an end phase.
3393: One can call VecLockWriteSet_Private(x,PETSC_TRUE) in the begin phase to lock a vector for exclusive
3394: access, and call VecLockWriteSet_Private(x,PETSC_FALSE) in the end phase to unlock the vector from exclusive
3395: access. In this way, one is ensured no other operations can access the vector in between. The code may like
3398: VecGetArray(x,&xdata); // begin phase
3399: VecLockWriteSet_Private(v,PETSC_TRUE);
3401: Other operations, which can not acceess x anymore (they can access xdata, of course)
3403: VecRestoreArray(x,&vdata); // end phase
3404: VecLockWriteSet_Private(v,PETSC_FALSE);
3406: The call can not be nested on the same vector, in other words, one can not call VecLockWriteSet_Private(x,PETSC_TRUE)
3407: again before calling VecLockWriteSet_Private(v,PETSC_FALSE).
3409: Level: beginner
3411: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop(), VecLockGet()
3412: @*/
3413: PetscErrorCode VecLockWriteSet_Private(Vec x,PetscBool flg)
3414: {
3417: if (flg) {
3418: 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");
3419: 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");
3420: else x->lock = -1;
3421: } else {
3422: 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");
3423: x->lock = 0;
3424: }
3425: return(0);
3426: }
3428: /*@
3429: VecLockPush - Pushes a read-only lock on a vector to prevent it from writing
3431: Level: deprecated
3433: .seealso: VecLockReadPush()
3434: @*/
3435: PetscErrorCode VecLockPush(Vec x)
3436: {
3439: VecLockReadPush(x);
3440: return(0);
3441: }
3443: /*@
3444: VecLockPop - Pops a read-only lock from a vector
3446: Level: deprecated
3448: .seealso: VecLockReadPop()
3449: @*/
3450: PetscErrorCode VecLockPop(Vec x)
3451: {
3454: VecLockReadPop(x);
3455: return(0);
3456: }
3458: #endif